text:
# Detection of invisible and crucial events: from seismic fluctuations to the war against terrorism

## I Introduction

The main aim of this paper is to prove the efficiency of a new method for the detection of crucial events that might have useful applications to the war against terrorism. This has to do with the search for rare but significant events, a theme of research that has been made of extreme importance by the tragedy of September 11. This method is applied here to defining the statistics of seismic main-shocks, as done in an earlier publication earlier . However, the emphasis here is more on the conceptual issues behind the interesting results obtained in Ref. earlier  than on their geophysical significance. In fact, the discussion of these conceptual issues aims at supporting the conjecture that the method has a wider range of validity. We shall help the reader to understand this general discussion with a dynamic model, originally proposed in Ref. memorybeyondmemory . We point out that this model was proposed for purposes different from geophysical applications. However, it is a case where the crucial events to detect are under our control, thereby making it possible for us to check the accuracy of the method of detection of invisible and crucial events that we propose here for a more general purpose, including the war against terrorism. Furthermore, for this model an analytical treatment has been recently found memory2 , supporting the claims that we make in this paper for the accuracy of the method of detection. For the reader’s convenience, the results on the seismic fluctuations are suitably reviewed, and discussed in the light of the more general perspective of this paper. We also review the model for seismic fluctuations, proposed in the earlier work of Ref. earlier . This model shares with the model of Ref.memorybeyondmemory  the property that the crucial events are imbedded in a sea of secondary events, but it allows us to reveal with accuracy the statistics of the crucial events for different mathematical reasons.

## II Crucial events, memory and predictability

The analysis shown in action on the seismic fluctuations, should serve the more general purpose of detecting the statistical properties of *crucial* events that are *invisible*. By crucial events we mean events causing other events, which would be predictable if the time occurrence of their causes were known. By invisible and crucial events we mean crucial events embedded in a sea of many other events, either caused by the invisible events or by environmental fluctuations. These secondary events play a camouflage action that makes it difficult to detect the crucial events with accuracy. We discuss first the property that the crucial events must have, regardless of whether they are visible or not. We want also to address the delicate issue of the extent to which the crucial events are predictable and the extent to which they are random, this being a problem that has caused much confusion in the past. Let us consider the following dynamic model. A particle moves in the interval $`I[0,1]`$. Its trajectory $`x(t)`$ is governed by the equation:

$$\frac{dx}{dt}=\alpha x^z.$$
(1)
The parameter $`\alpha `$ is a positive number, which has to be kept much smaller than $`1`$, if the integration time step is $`1`$. This fundamental equation serves the purpose of generating non-Poisson statistics. In fact, as we shall see, $`z=1`$ generates Poisson statistics, while the wide dominion of non-Poisson statistics is given by $`z>1`$. When the particle reaches the border $`x=1`$ it is injected back to a new initial condition $`x(0)>0`$, with uniform probability. Now, let us imagine that the times of sojourn within the interval are recorded via direct observation. For tutorial purposes, we assume these random events to be visible. We shall get the time series

$$\{\tau \}\tau (1),\tau (2),\mathrm{}$$
(2)
Let us consider the initial condition $`x(0)`$. The time spent by the particle, moving from this initial condition, within the interval $`I`$, before reaching the border is

$$\tau =\frac{1}{\alpha }\left[\frac{1}{1z}\frac{x(0)^{1z}}{1z}\right],$$
(3)
as predicted by the time integration of Eq. (1). The connection between the waiting-time distribution and the injection process into the initial condition is given by

$$\psi (\tau )d\tau =p(x(0))dx(0).$$
(4)
In the case of a uniform injection, $`p(x(0))=1`$, inserting Eq. (3) into Eq. (4) yields after some algebra

$$\psi (\tau )=(\mu 1)\frac{T^{\mu 1}}{(T+\tau )^\mu },$$
(5)
where

$$\mu \frac{z}{z1}$$
(6)
and

$$T=\frac{1}{\alpha (z1)}.$$
(7)
It is worth devoting some comments to this result. Let us imagine that we convert the time serie $`\{\tau \}`$ into a diffusion process with the same prescription as that adopted in the earlier work on the seismic fluctations in Southern California. This means that the random walkers makes a jump ahead by a given quantity, equal to $`1`$, for instance, at the times $`t(1)=\tau _1`$, $`t(2)=\tau _1+\tau _2`$, $`t(3)=\tau _1+\tau _2+\tau _3`$, and so on. An ensemble of random walkers obeying the same prescription (see Ref. giacomo , for details on how to derive this ensemble from the single sequence $`\{\tau \}`$) undergoes a diffusion process that in the specific case $`2<\mu <3`$ yields a diffusion process. The probability distribution function, $`p(x,t)`$, for this process, in the time asymptotic limit is expected to obey the scaling condition

$$p(x,t)=\frac{1}{t^\delta }F(\frac{xwt}{t^\delta }),$$
(8)
where $`w`$ is the mean velocity produced by the walking rule adopted and $`\delta `$ is the scaling index given by

$$\delta =\frac{1}{\mu 1}.$$
(9)
$`F(y)`$ is an asymmetric function of $`y`$, whose detailed analytical form is discussed in Ref. giacomo . Ref. giacomo  discusses other walking prescriptions, and physical conditions different from $`2<\mu <3`$ as well. For simplicity, in this paper we discuss only the earlier walking prescription and the case where $`\mu <3`$, so as to create a strong departure from ordinary statistical mechanics, namely from the condition where the second moment of $`\psi (\tau )`$ is finite. We also set the condition $`\mu >2`$ which keeps the system far from the condition of a diverging first moment. Why do we assign to $`\psi (\tau )`$ this condition of strong departure from ordinary statistical mechanics? We shall answer this important question after making the reader familiar with the intriguing issue of the memory emerging from the breakdown of the Poisson statistics. This is a poorly understood property, in spite of the fact that 32 years ago Bedeaux, Lindenberg and Shuler katja  wrote a clarifying paper on this subject. We have seen that our method of analysis rests on turning a time series into a diffusion process. If we imagine the one-dimensional axis on which this diffusion process is realized as a chain of infinite discrete sites, we can denote the state at time $`t`$ of the diffusing system through the vector $`𝐩(t)`$, with $`p_i(t)`$ denoting the probability for the random walker to be at the $`i`$-th site at time $`t`$. Thus, it is legitimate to ask the important question of whether the knowledge of $`𝐩(t)`$ allows us to determine $`𝐩(t^{})`$ with $`t^{}>t`$. In other words, the question is the following: does the information $`𝐩(t)`$ at a given time $`t`$ allow us to predict the state of the system at a later time? We want to point out that the question refers to a set of random walkers, not to a single random walker, whose walk, at a time scale larger than the first moment of the waiting time distribution $`\psi (\tau )`$, is certainly unpredictable.

A careful reading of Ref. katja  allows us to answer this question with this apparently striking statement: this is possible only in the Poisson case gerardo . In the non-Poisson case an infinitely extended memory emerges. The Poisson case is the only one where the state $`𝐩(t)`$, with $`t>0`$, determines the time evolution of the system of interest. In all the other cases, the future time evolution of the system also depends on $`𝐩(t^{\prime \prime })`$, with $`t^{\prime \prime }<t`$. In other words, the system time evolution retains memory of the initial condition $`𝐩(0)`$ forever. It is evident that the concept of crucial event implies a departure from the Poisson condition. In fact, the probability of occurrence of a main-shock is expected to have memory, this corresponds to the fact that the geophysical processes, responsible for the main-shocks, do not generate random fluctuations, but long-range correlation.

In the case of the terrorist network, we conjecture that the crucial events, having either ideological or religious origin, and so historical roots, are driven by non-Markovian master equations: This implies that the irrelevant degrees of freedom, playing the role of a thermal bath, are not equivalent to white noise, as it would be in the case of Poisson statistics katja . The choice of the condition of $`2<\mu <3`$ makes the resulting diffusion process depart dramatically from the Gaussian state, thereby assigning to the crucial events, either main-shocks or main secret events of the terrorist network, a condition of striking departure from ordinary statistical mechanics. In other words, we conjecture that the crucial events, which, although invisible, influence cascades of secondary events, are located in a basin of attraction of anomalous rather than normal statistics, an assumption that fits the traditional wisdom of the researchers in the field of complexity. It has to be pointed out that from a formal point of view a condition of infinite memory is realized by $`\mu <\mathrm{}`$, without necessarily implying $`\mu <3`$. However, the condition $`\mu >3`$ would not make the crucial event generate a visibly anomalous diffusion, and an even more sensitive procedure should be planned, to discover the existence of this kind of crucial events. Thus, the condition $`\mu <3`$, which, as earlier pointed out, seems to be a plausible property of complex systems, corresponds to a case where the procedure illustrated in this paper, is already adequate, in the present form, to reveal their existence. Before ending this Section, we must clarify a problem that is a frequent source of confusion. The infinite memory associated with non-Poisson statistics might be mistaken as a way to make prediction.

We have to point out that the infinite memory is a concept referring to probabilities, or to a set of walkers. The concept of Gibbs ensemble, although fundamental for statistical mechanics, is based on the assumption that many identical copies of the system are available to us. Actually, we study only single systems. For instance, the predictability of earthquakes, implies that knowing that a crucial event occurred at time $`t=t_1`$, we can predict that the next will take place at a time $`t=t_2>t_1`$. The time distance $`t_2t_1`$ cannot be predicted, if the events under study are crucial in the sense earlier defined. However, if Eq. (1) were a reliable model for the process under study, one might conjecture that a specific observation of the geophysical motion is equivalent to informing us about the new initial condition, after the back injection taking place at $`t=t_1`$. The instant of the back injection is the genuinely random event. The laminar motion ensuing this random event is deterministic and, consequently, compatible with predictability, at least in principle. Why do we leave room for randomness, in the moment of selection of the new initial condition? This is equivalent to associating crucial events to randomness, and a thorough discussion of this issue beyond the limits of this paper. If we adopt the usual view that randomness is an expression of our ignorance about the infinitely many and irrelevant degrees of freedom in a system, this choice is equivalent to a drastic simplification of the problem under study. Going beyond that would be equivalent to predict the occurrence time of main-shocks, in the case of seismic fluctuations, and of terrorist actions, in the case of the war to terrorism. For the time being, our purpose is much less ambitious.

## III memory beyond memory

Notice that the title of the paper of Ref. memorybeyondmemory , *memory beyond memory*, is probably incomprehensible to all those who do not know the fundamental work of Ref. katja . On the basis of the results of Ref. katja  we can explain the meaning of this title. The majority of events under observation in Ref. memorybeyondmemory  are not crucial events. The crucial events, which are rare, are imbedded in a sea of secondary events, also called pseudo-events. These secondary events are influenced by the crucial events and play a camouflage role that makes the really crucial events invisible. As a consequence of being secondary, the pseudo-events have memory of the crucial events influencing them.

It is worth remarking that in the case of seismic events the crucial events are the main shocks and the secondary events are the Omori swarms of aftershocks. We make an important conjecture: the case of the terrorist network rests on the picture of the passive supporters of terrorist activities. These supporters generate events that, although secondary, are dependent on the main terrorist events, of which they bear memory. This is the memory of the second kind, the memory of first kind being, as pointed out in the earlier Section, the memory corresponding to the non-Poisson statistics of the crucial events. The terrorists trigger events, with either religious or ideological memory, and these crucial events influence secondary events, the action of passive supporters, which are characterized by *memory beyond memory*. Although, the term *memory beyond memory* was originally coined for the completely different purpose of helping the search for the key physiological processes behind heart-beating memorybeyondmemory , we find it to be especially adequate to describe a procedure of statistical analysis aiming at helping the Intelligence Community in the war against terrorism. For all these reasons, it is convenient to review a model that was originally proposed to illustrate the origin of memory of the second type memorybeyondmemory . The model consists of two particles. The first particle is the visible one. We maintain the same dynamic rule as that established by Eq. (1), though now we change the role of this particle from the generator of crucial events to the generator of pseudo-events. This means that we keep using the visible particle to generate events, and consequently the time series to study, with a different back injection rule, though, for the purpose of generating events that are not random. To do that, following Ref. memorybeyondmemory , we introduce a second particle ruled by an equation of the same kind as Eq. (1). However, while Eq. (1) refers to events that we can monitor, now the second particle refers to hidden events. Therefore we refer to this particle as the invisible particle. In conclusion, we describe this model by means of the following set of equations:

$`{\displaystyle \frac{dx_{vis}}{dt}}=\alpha x_{vis}^z,`$ (10)
$`{\displaystyle \frac{dx_{invis}}{dt}}=\beta x_{invis}^\zeta .`$ (11)
We assume that the dynamics of the visible particle are much faster than the dynamics of the invisible particle. Thus, the visible particle gets to the border and is injected back many times before the occurrence of the leading, or crucial, event. The crucial event occurs when the invisible particle reaches the border and is injected back randomly to a new initial condition in the interval $`I`$. Before the occurrence of this crucial event, the visible particle has been injected back following a very simple deterministic prescription. In the earlier work of Ref. memorybeyondmemory  to check the efficiency of our method of analysis we have made the assumption that the initial condition is always the same, and it is changed randomly only when the invisible particle is injected back.

At this stage, we wonder if it is possible to distinguish the crucial events from the surrounding pseudo events; in particular we wonder if a statistical method of analysis exists that detects the waiting time distribution of the crucial events. The answer is positive, and can be found in the paper of Ref.memorybeyondmemory . First of all, we have to convert the time series into a diffusion process. According to the prescriptions of Ref. giacomo , we evaluate the Shannon entropy of this diffusion process. This is why this technique of analysis is called Diffusion Entropy (DE) method. As is well known, the distance $`x`$ travelled by the walkers is related to time by the relation $`xt^\delta `$, where $`\delta `$ is termed diffusion index. If the diffusion process is not the sum of uncorrelated fluctuations, the scaling parameter $`\delta `$ departs from the prescription of ordinary statistical mechanics, namely, from $`\delta =0.5`$. The DE method is an efficient way to determine $`\delta `$ giacomo .

In the case where the visible events are not correlated, and the walking rule of Section II is adopted, the scaling $`\delta `$, determined by means of the DE method, and the power index $`\mu `$ of the waiting time distribution $`\psi (\tau )`$ of Eq. (5), are related by means of Eq. (9). The violation of this crucial condition suggests that the events under observation are not genuine events, but rather pseudo-events, bearing, as a consequence of that, memory of the second type. Actually, the way to proceed is as follows. We evaluate numerically the waiting time distribution $`\psi _{exp}(\tau )`$, by running the two-walkers model. The observation of visible events determines the waiting time distribution

$$\psi _{exp}(\tau )=(\mu ^{}1)\frac{T^{\mu ^{}1}}{(T+\tau )^\mu ^{}}.$$
(12)
We do not address here the interesting problem of establishing $`\mu ^{}`$ as a function of the parameter of the two-walkers model. This is not crucial for the discussion of this paper. Let us limit ourselves to noticing that $`\mu ^{}>\mu `$. The numerical results of Ref. memorybeyondmemory  show that $`\delta `$ does not have anything to do with $`1/(\mu ^{}1)`$. These numerical results rather prove the attractive fact that Eq. (9) applies, but with $`\mu `$ denoting the power coefficient of $`\psi (\tau ^{[m]})`$, and $`\tau ^{[m]}`$ the time distance between two consecutive crucial events (the subfix $`m`$ here stands for main events, in analogy with the definition used in earlier ). In other words, the scaling coefficient $`\delta `$, detected by means of the DE method, reveals an important statistical property of crucial and invisible events. Let us summarize the procedure that we propose to detect the statistical properties of invisible and crucial events. First of all, we adopt an experimental view, and we derive from the real sequence under study the waiting time distribution $`\psi _{exp}(\tau )`$, referring to the time distance between two consecutive events. If we find that the waiting time distribution is not exponential, we have a first indication of complexity. If the distribution is an inverse power law, we record the power law index, $`\mu ^{}`$. Then we use the DE method to measure the scaling parameter $`\delta `$. The condition $`\delta =1/(\mu ^{}1)`$ is a plausible indication that we are observing a time sequence of significant events. If $`\delta `$ significantly departs from $`1/(\mu ^{}1)`$, there are good reasons to believe that

$$\mu =1+\frac{1}{\delta }$$
(13)
is a reliable indicator of the complexity of invisible and crucial events. It is important to note that this important conclusion is supported by the analytical treatment of Ref. memory2 . The authors of Ref. memory2  shows that in the long-time limit the memory of the pseudo-events is lost, and the process under study becomes equivalent to a Lévy flight, corresponding to the power index $`\mu `$ of the crucial events.

## IV The Omori’s law as a source of pseudo-events

In this and in the next Section we review the work of Ref. earlier  for the main purpose of proving that the results of this paper are a realization of the method for the search of invisible and crucial eventes illustrated in Section III. In Fig. 1 we report the sketch of the typical earthquakes frequency *vs* time in the catalog that we shall consider in the next Section. By $`\tau _i=t_{i+1}t_i`$ we indicate the time interval between an earthquake and the next. Each peak of frequency (cluster) in figure includes the time location of a main-shock. The time interval between one peak and the next is reported in figure and is denoted by the symbol $`\tau _i^{[m]}`$, where the superscript $`m`$ stands for main-shock, since the main-shocks are the main event in the case of seismic fluctuations. According to the definition of crucial events given in Section I, we must make the assumption that two different $`\tau ^{[m]}`$’s are not correlated, i.e. that the correlation function is:

$$\tau _i^{[m]}\tau _j^{[m]}=\left(\tau ^{[m]}\right)^2\delta _{i,j}.$$
(14)
Note that with the symbol $`\tau ^{[m]}`$ we denotes distances between two consecutive crucial events. Thus the corresponding waiting time distribution is equivalent to that distribution of Eq. (5).

The experimental determination of this distribution would imply the adoption of a way to identify the main-shocks. Although the geologists might suggest reliable criteria for their identification, for instance through the magnitude, with the use of our method we can determine their statistical properties without identifying them. Thus, the symbols $`\tau ^{[m]}`$ denote distances between consecutive events that are assumed to be invisible. One of the models adopted to describe the time distribution of earthquakes is the Generalized Poisson (GP) model shlien ; gardner ; gasperini ; console ; godano . Basically the GP model assumes that the earthquakes are grouped into temporal clusters of events and these clusters are not long-range correlated: in fact the clusters are distributed at random in time and therefore the time intervals between one cluster and the next one follow a Poisson distribution. On the other hand, the intra-cluster earthquakes are correlated in time as it is expressed by the Omori’s law omori ; utsu , an empirical law stating that the main-shock, i.e. the highest magnitude earthquake of the cluster, occurring at time $`t_0`$ is followed by a swarm of correlated earthquakes (after-shocks) whose number (or frequency) $`n(t)`$ decays in time as a power law, $`n(t)(tt_0)^p`$, with the exponent $`p`$ being very close to $`1`$. The Omori’s law implies bak  that the distribution of the time intervals between one earthquake and the next, denoted by $`\tau `$, is a power law $`\psi (\tau )\tau ^p`$. This property has been recently studied by the authors of Ref. bak  by means of a unified scaling law for $`\psi _{L,M}(\tau )`$, the probability of having a time interval $`\tau `$ between two seismic events with a magnitude larger than $`M`$ and occurring within a spatial distance $`L`$. This has the effect of taking into account also space and extending the correlation within a finite time range $`\tau ^{}`$, beyond which the authors of Ref. bak  recover Poisson statistics. Let us discuss the GP model in the light of the general remarks of Sections II and III. The Poisson assumption about the distribution of time distances between one main-shock and the next is equivalent to assigning no memory to the geophysical process responsible for the main-shock. This conflicts with our definition of crucial events and with our conviction that the crucial events, as unpredictable as the time duration of a laminar region is, cannot be determined by erratic bath fluctuations. The GP model, if supported by the statistical analysis of data, would imply that our definition is not correct, and that crucial events can be generated also from within ordinary statistical mechanics. This would conflict also with the tenets of complexity, which seem to connect cooperation and inverse power law relaxation. In Section V we shall prove that the GP model must be dismissed. In fact, using the method of statistical analysis reviewed in this paper, it is shown earlier  that the asymptotic scaling generated by the GP scaling would be $`\delta =0.5`$, which in fact corresponds to the prescription of Section II, when $`\mu >3`$. It is worth recalling that the Poisson condition sets the exponential decay of $`\psi (\tau )`$, and thus $`\mu =\mathrm{}`$. The statistical analysis of real data, discussed in Section V, will prove that $`\delta =0.94`$, that the GP model is incorrect, and that we cannot rule out the possibility that the main-shocks fulfill our definition of crucial events.

## V Data and results

The catalog we have studied covers the period 1976-2002 in the region of Southern California spanning $`20^0`$ N -$`45^0`$ N latitude and 100<sup>0</sup> W 125<sup>0</sup> W longitude scsn . This region is crossed by the most seismogenetic part of the San Andrea fault, which accommodates by displacement the primarily strike-slip motion between the North America and the Pacific plates, producing velocities up to $`47mm/yr`$ turcotte . The total number of recorded earthquakes in the catalog is $`383687`$ and includes the June 28 1992 Landers earthquakes (M = 7.3), the January 17 1994 Northridge earhquake (M = 6.7), and the October 16 1999 Hector Mine earthquake (M = 7.1). Geophysical observations point out that these large earthquakes have triggered a widespread increase of seismic events at remote distances in space and in time hill ; parsons . The coupling of the sources of stress change (i.e. large earthquakes occurrence) and seismicity triggering mechanisms is a primary target of geophysical investigations, and, as shown below, is revealed by the DE analysis. In Fig. 2 we report the results of the DE method. The analysis was performed by setting $`\xi (t)=1`$ when an earthquake occurs at time $`t`$ (independently of whether it is a main-shock or an after-shock), and $`\xi (t)=0`$ if no earthquake happens. By means of the full circles we denote the entropy $`S(t)`$ as a function of time when all the seismic events of the catalog are considered (independently of their magnitude $`M`$). After a short transient, the function $`S(t)`$ is characterized by a linear dependence on $`\mathrm{ln}t`$. A fit in the linear region gives a value of the scaling parameter $`\delta =0.94\pm 0.01`$ at $`95\%`$ of confidence level. We next consider only the earthquakes with magnitude larger than a fixed value $`\overline{M}=2,3,4`$. We see that, regardless of the value of the threshold $`\overline{M}`$ adopted, the function $`S(t)`$ is characterized by the same long-time behavior with the same slope. This indicates that we are observing a property of the time location of large earthquakes. This leads us to conclude that the time intervals between two large events fit the distribution of Eq.(5), with the value of $`\mu `$ related to $`\delta `$ through Eq.(9), $`\mu =2.06\pm 0.01`$. Our conclusion is also supported by other two numerical analyses based on different prescriptions to construct the diffusion process. The former rests on assuming $`\xi (t)`$ equal to the magnitude $`M`$ of the earthquake, at each time when an earthquake occurs. The latter sets with equal probability either $`\xi (t)=1`$ or $`\xi (t)=1`$ when an earthquake occurs giacomo . Both the methods give the same exponent $`\mu =2.06\pm 0.01`$. In conclusion, the statistical analysis of real data rules out the GP model, which would conflict with our definition of crucial events. Thus, there is still room for the main-shocks to fit our definition of crucial events. The authors of Ref. earlier  prove that, under the stationary condition, they do. Are these crucial events also invisible? This is a question of fundamental importance for the war against terrorism. The answer is that the main-shocks are not quite invisible. For professional geologists it is possible to identify all of them. However, our method of analysis works, regardless of whether the crucial events are invisible or not. In fact, Fig. 2 (a) shows that the asymptotic properties of the entropy indicator are independent of the threshold $`M`$ adopted. Since the magnitude is a distinctive property of the main-shocks, we conclude that the results of our analysis do not imply that the main-shocks are identified. This is the reason of our conviction that our method can be successfully used even when the crucial events are quite invisible.

## VI Generalized Poisson (GP) and Long-range (LR) model

We now illustrate how the DE method works on two artificial earthquakes time series: the first generated by means of the GP model, and the second generated by a new model, the Long-Range (LR) model, that we propose as a better model to reproduce the properties of the catolog considered. In the LR model the earthquakes are grouped into temporal clusters, and, as in the GP model, the number of earthquakes in a cluster follows the Pareto law, i.e. a power law distribution with exponent equal to 2.5 console . The events within the same cluster are distributed according to the Omori’s law: the interval $`\tau `$ follows a power law with exponent $`p=1`$. However, in the LR model the time distance $`\tau ^{[m]}`$ between one cluster and the next follows a power law with exponent $`\mu =2.06`$, rather than a Poisson prescription as in the GP. Notice that this value of $`\mu `$ is close to the border between stationary and non-stationary condition giacomo . The two sequences have the same time length. We choose the number of clusters in order to have the same total number of earthquakes as in the real data note1 . The result of the DE on the artificial sequences is reported in Fig. 2(b). The GP model is characterized by a long-time behavior that, as expected, fits very well the prescription of ordinary statistical mechanics, with $`\delta =0.5`$. The LR model yields the quite different scaling $`\delta =0.94`$. It is also clear that the LR model yields a behavior qualitatively similar to that produced by the real data of Fig. 2 (a) as well as the same scaling parameter $`\delta =0.94`$, while the GP fail reproducing both properties. Note that from the reasons why the DE method reveals the genuine statistical properties of the crucial events in this case, are different from those justifying the efficiency of the method in the case of the model of Section III. In this case, the true reason seems to be that for large distances between one main-shock and next the pseudo-events tend to concentrate immediately after the last main-shock with a so low density immediately before the occurrence of the next as to create a condition where the number of pseudo-events is essentially independent of the length of the laminar region. In the case of the model of Section III, on the contrary, the number of pseudo-events is proportional to the length of the laminar region.

## VII Conclusions

It is the time for us to balance the detection of crucial and invisible events in general, with possible applications to the war against terrorism. Our definition of crucial events implies a strong departure from Poisson statistics. In the recent literature there is a general agreement about the fact that complex networks, including the terrorist network, are scale-free systems. The authors of another paper of these Proceedings linguistics  show that there exists a connection between the scale-free condition and non-Poisson statistics. This seems to support the definition of crucial event adopted in this paper. On the other hand, there exists an interesting connection with the Small Words theory illustrated by Latora and Marchiori vito , a paper of these Proceedings, explicitly devoted to the war against terrorism. This paper, in turn, rests on a perspective that is related to the sociological picture illustrated Elliott and Kiel euel  for the same purpose. It is worth mentioning that a local version of the DE method can be applied to the timely detection of toxicants massi . In conclusion, the present paper belongs to a set of contributions to these Proceedings linguistics ; vito ; euel ; massi , which might bear benefits to a program of research to combat terrorism. As to the detection of invisible and crucial events for specific purpose of the war to terrorism, the procedure to follow depends on the data to analyze, and on the forms, under which they will be made available to the investigators. Nevertheless, with the present paper, we are convinced that at least the first few steps of the search for crucial and invisible events, take a clear shape. It seems to be evident what the first step of this procedure will be the evaluation of $`\psi _{exp}(\tau )`$, and of the corresponding power index, denoted by the symbol $`\mu ^{}`$ in this paper. The second step will be the evaluation of $`\delta `$, by means of the DE method, and the comparison of $`\delta `$ with $`1/(\mu ^{}1)`$. If the two values do not coincide, and the difference is larger than the statistical error, this has to be thought of as a plausible indication that there are invisible and crucial events involved. Then, we shall have to decide whether or not recourse can be done to the simple prescription of Eq.(13) to establish the degree of complexity of the invisible crucial events. This will require further research work determined by the specific nature of the data to analyze.

PG acknowledges support from ARO, through Grant DAAD19-02-0037.


format: text

text:
# Observation of a Narrow Meson Decaying to 𝑫⁺_𝒔𝝅^𝟎𝜸 at a Mass of 2.458 𝐆𝐞𝐕/𝒄^𝟐

thanks: Deceased

The BABAR Collaboration

## Abstract

A narrow state, which we label $`D_{sJ}(2458)^+`$, with mass $`2458.0\pm 1.0(\text{stat.})\pm 1.0(\text{syst.})`$ $`\mathrm{MeV}/c^2`$, is observed in the inclusive $`D_s^+\pi ^0\gamma `$ mass distribution in 91 $`\mathrm{fb}^1`$ of $`e^+e^{}`$ annihilation data recorded by the BABAR detector at the PEP-II asymmetric-energy $`e^+e^{}`$ storage ring. The observed width is consistent with the experimental resolution. The data favor decay through $`D_s^{}(2112)^+\pi ^0`$ rather than through $`D_{sJ}^{}(2317)^+\gamma `$. An analysis of $`D_s^+\pi ^0`$ data accounting for the influence of the $`D_{sJ}(2458)^+`$ produces a $`D_{sJ}^{}(2317)^+`$ mass of $`2317.3\pm 0.4(\text{stat.})\pm 0.8(\text{syst.})`$ $`\mathrm{MeV}/c^2`$.

preprint: BABAR-PUB-03/030preprint: SLAC-PUB-10215

Interest in the spectrum of charmed mesons has been heightened by the discovery by this collaboration Aubert et al. (2003) of a narrow state, produced in $`e^+e^{}c\overline{c}`$ collisions at the PEP-II collider, decaying to $`D_s^+\pi ^0`$ cpf , with mass 2317 $`\mathrm{MeV}/c^2`$, approximately 41 $`\mathrm{MeV}/c^2`$ below the $`DK`$ mass threshold. This state, the $`D_{sJ}^{}(2317)^+`$, has been confirmed by CLEO Besson et al. (2003) and Belle Mikami et al. (2003); Krokovny et al. (2003). Along with the $`D_{sJ}^{}(2317)^+`$, we noted Aubert et al. (2003) the presence of a narrow peak in the $`D_s^+\pi ^0\gamma `$ mass distribution near 2.46 $`\mathrm{GeV}/c^2`$. Because this signal is near the kinematic overlap of the $`D_{sJ}^{}(2317)^+\gamma `$ and $`D_s^{}(2112)^+\pi ^0`$ systems, special attention is required to remove associated backgrounds and to distinguish between the two possible decay modes. Such an analysis is the subject of this paper.

This state near 2.46 $`\mathrm{GeV}/c^2`$ has been seen by CLEO Besson et al. (2003) and Belle Mikami et al. (2003) in the inclusive $`D_s^+\pi ^0\gamma `$ mass spectrum and by Belle Krokovny et al. (2003) in exclusively reconstructed $`B`$ decays.

To investigate the $`D_s^+\pi ^0\gamma `$ spectrum, we study $`D_s^+`$ candidates from $`e^+e^{}c\overline{c}`$ (at a center-of-mass energy near 10.6 GeV) that decay to $`K^{}K^+\pi ^+`$. Particle identification is used to provide clean samples of charged $`K`$ and $`\pi `$ candidates, which are combined using a geometric fit to a common vertex. Backgrounds are suppressed by selecting decays to $`\overline{K}{}_{}{}^{0}K_{}^{+}`$ and $`\varphi \pi ^+`$. A description of this sample and additional details can be found elsewhere Aubert et al. (2003). Events with $`1.954<m(K^{}K^+\pi ^+)<1.981`$ $`\mathrm{GeV}/c^2`$ are taken as $`D_s^+`$ candidates.

A candidate $`\pi ^0`$ is formed by constraining a pair of photons each with energy greater than 100 MeV to emanate from the intersection of the $`D_s^+`$ trajectory with the beam envelope, performing a one-constraint fit to the $`\pi ^0`$ mass, and requiring a fit probability greater than 5%. A given event may yield several acceptable $`\pi ^0`$ candidates. We retain only those candidates for which neither photon belongs to another otherwise acceptable $`\pi ^0`$.

Each $`D_s^+`$ candidate is combined with all combinations of accompanying $`\pi ^0`$ candidates with momentum greater than 300 $`\mathrm{MeV}/c`$ and photon candidates of energy greater than 100 MeV. To suppress background, photons that belong to any $`\pi ^0`$ candidate are excluded and we require the momentum, $`p^{}`$, of each $`D_s^+\pi ^0\gamma `$ combination in the $`e^+e^{}`$ center-of-mass frame to be greater than 3.5 $`\mathrm{GeV}/c`$. The last requirement also removes any $`D_s^+\pi ^0\gamma `$ combination from $`B`$ decay.

The $`D_s^+\pi ^0\gamma `$ invariant mass distribution is shown in Fig. 1a. A clear enhancement is observed near 2.46 $`\mathrm{GeV}/c^2`$. The background underneath this peak is from several sources, which can be described in terms of mass differences defined as

$`\mathrm{\Delta }m_\gamma `$ $``$ $`m(D_s^+\gamma )m(D_s^+)`$ (1)
$`\mathrm{\Delta }m_{\pi ^0}`$ $``$ $`m(D_s^+\gamma \pi ^0)m(D_s^+\gamma ).`$ (2)
A scatter plot of the data is shown in Fig. 1b. Particular background patterns are visible: $`D_s^{}(2112)^+D_s^+\gamma `$ decay combined with an unassociated $`\pi ^0`$, which appears as a horizontal band, and $`D_{sJ}^{}(2317)^+D_s^+\pi ^0`$ decay combined with an unassociated $`\gamma `$, which appears as a band that is almost vertical.

To demonstrate the existence of a signal above these backgrounds, the upper histogram of Fig. 1c shows $`D_s^+\pi ^0\gamma `$ combinations in the $`D_s^{}(2112)^+`$ signal region, and the gray histogram, scaled to the area of the signal region, corresponds to the two $`D_s^{}(2112)^+`$ sidebands. We conclude that a signal for a state decaying to $`D_s^+\pi ^0\gamma `$ exists over a background resulting from $`D_{sJ}^{}(2317)^+`$ and an unassociated $`\gamma `$. This background peaks at a mass slightly higher than that of the signal. A Gaussian fit to the subtracted mass distribution (Fig. 1d) indicates a narrow signal at $`\mathrm{\Delta }m_{\pi ^0}=346.2\pm 0.9`$ $`\mathrm{MeV}/c^2`$ (statistical error only).

The state corresponding to this signal, which we label $`D_{sJ}(2458)^+`$, may decay to $`D_s^+\pi ^0\gamma `$ through $`D_s^{}(2112)^+\pi ^0`$ or $`D_{sJ}^{}(2317)^+\gamma `$. To disentangle these modes and reliably extract the parameters of the signal, we apply an unbinned maximum likelihood fit simultaneously to the $`D_s^+\pi ^0\gamma `$, $`D_s^+\pi ^0`$, and $`D_s^+\gamma `$ invariant masses of all $`D_s^+\pi ^0\gamma `$ combinations using the channel likelihood method Condon and Cowell (1974). This fit describes the probability density function of the two $`D_{sJ}(2458)^+`$ decay channels as the product of a Gaussian shape in the $`D_s^+\pi ^0\gamma `$ mass distribution and a Gaussian shape projected into the $`D_s^+\pi ^0`$ or $`D_s^+\gamma `$ mass axes, as appropriate. Because the daughter resonances are narrow, interference between the two $`D_{sJ}(2458)^+`$ decay modes cannot be resolved, and so is ignored.

Sources of background in the $`D_s^+\pi ^0\gamma `$ spectrum included in the fit are purely combinatorial background ($`D_s^+`$ meson combined with an unassociated $`\pi ^0`$ and $`\gamma `$), $`D_s^{}(2112)^+D_s^+\gamma `$ decay combined with an unassociated $`\pi ^0`$, and $`D_{sJ}^{}(2317)^+D_s^+\pi ^0`$ decay combined with an unassociated $`\gamma `$. The fit also includes a contribution from $`D_{sJ}(2458)^+D_s^{}(2112)^+\pi ^0`$ decay but with an unassociated $`\gamma `$ replacing the $`\gamma `$ from $`D_s^{}(2112)^+`$ decay. The fit determines the relative size of the background and signal contributions, the mass and width of the $`D_{sJ}(2458)^+`$, and the $`D_{sJ}^{}(2317)^+`$ mass.

The likelihood fit is validated using Monte Carlo (MC) simulation. This simulation includes $`e^+e^{}c\overline{c}`$ events and all known charm states and decays, including the $`D_{sJ}^{}(2317)^+`$ and the signal under study. The generated events were processed by a detailed detector simulation Agostinelli et al. (2003) and subjected to the same reconstruction and event-selection procedure as the data.

As shown in Fig. 2a, the fit provides a good description of the $`D_s^+\pi ^0\gamma `$ mass distribution observed in the data. The $`D_{sJ}(2458)^+`$ signal for a particular decay mode can be isolated by calculating a weight for each $`D_s^+\pi ^0\gamma `$ combination proportional to the relative likelihood contributed by the decay mode of interest. Distributions of events so weighted can be compared to the likelihood function to validate the fit. This is shown in Figs. 2b and 2c. A $`\chi ^2`$ probability calculation gives 22%, 74% and 11% for fig. 2a, b and c respectively. The resulting yield of correctly reconstructed $`D_{sJ}(2458)^+D_s^{}(2112)^+\pi ^0`$ ($`D_{sJ}(2458)^+D_{sJ}^{}(2317)^+\gamma `$) decays is $`195\pm 26`$ ($`0\pm 23`$), consistent with the fit shown in Fig. 1d. Excluding the $`D_{sJ}(2458)^+`$ from the likelihood fit decreases the log likelihood by approximately 57, corresponding to a significance of more than $`10`$ standard deviations. The fit yields a $`D_{sJ}(2458)^+`$ mass of $`2458.0\pm 1.0`$ $`\mathrm{MeV}/c^2`$ with an rms width of $`8.5\pm 1.0`$ $`\mathrm{MeV}/c^2`$.

The likelihood fit uses the shapes of the $`D_s^+\pi ^0`$ and $`D_s^+\gamma `$ mass distributions to distinguish between the two possible decay modes, $`D_{sJ}(2458)^+D_s^{}(2112)^+\pi ^0`$ and $`D_{sJ}(2458)^+D_{sJ}^{}(2317)^+\gamma `$. These shapes are influenced by the kinematic constraints of $`D_{sJ}(2458)^+`$ decay shown in Fig. 3a. Figs. 3b–3c show the sideband-subtracted $`D_s^+\pi ^0`$ and $`D_s^+\gamma `$ mass projections compared with MC simulations of the two hypotheses (scaled to match the data yield). The $`D_{sJ}(2458)^+D_s^{}(2112)^+\pi ^0`$ decay mode (solid histograms) produces a narrow $`D_s^+\gamma `$ mass distribution and a wide $`D_s^+\pi ^0`$ mass distribution. In contrast, the $`D_{sJ}(2458)^+D_{sJ}^{}(2317)^+\gamma `$ decay mode (dashed histograms) produces a wide $`D_s^+\gamma `$ mass distribution and a narrow $`D_s^+\pi ^0`$ mass distribution. Figures 3b and 3c show that the $`D_{sJ}(2458)^+D_s^{}(2112)^+\pi ^0`$ hypothesis is in better agreement with the data.

Our previous measurement Aubert et al. (2003) of the $`D_{sJ}^{}(2317)^+`$ mass using the decay $`D_{sJ}^{}(2317)^+D_s^+\pi ^0`$ did not explicitly consider background from $`D_{sJ}(2458)^+D_s^{}(2112)^+\pi ^0`$ decay. This background peaks in the $`D_s^+\pi ^0`$ mass spectrum just below the $`D_{sJ}^{}(2317)^+`$ mass. Shown in Fig. 2d is the $`D_s^+\pi ^0`$ invariant mass distribution for a sample of $`D_s^+`$ candidates combined with all $`\pi ^0`$ candidates, with $`p^{}>3.5`$ $`\mathrm{GeV}/c`$. Superimposed on this distribution is a binned fit that includes the contribution from the $`D_{sJ}(2458)^+`$ as estimated from MC simulation and a quadratic background function. The result is a $`D_{sJ}^{}(2317)^+`$ yield of $`1022\pm 50`$ events, a mass of $`2317.3\pm 0.4`$ $`\mathrm{MeV}/c^2`$, and measured rms width $`7.3\pm 0.2`$ $`\mathrm{MeV}/c^2`$. These results are an improvement over our earlier measurement Aubert et al. (2003).

We divide the sources of systematic uncertainty in the $`D_{sJ}(2458)^+`$ and $`D_{sJ}^{}(2317)^+`$ mass values and production rates into three categories. The first category is associated with the fit procedure. Likelihood fits to MC samples that include samples of $`D_{sJ}(2458)^+D_s^{}(2112)^+\pi ^0`$ and $`D_{sJ}(2458)^+D_{sJ}^{}(2317)^+\gamma `$ decays correctly reproduce the given sample sizes within statistical errors. The average values of the fit results obtained using statistically distinct MC samples corresponding to the measurements in the data are used to place limits on any fit bias.

We obtain the background distribution used in the likelihood from a random selection of $`D_s^+`$, $`\pi ^0`$, and $`\gamma `$ candidates taken from the MC $`D_s^+\pi ^0\gamma `$ sample. To test our sensitivity to this distribution, various selection requirements are altered within reasonable bounds to provide alternate background samples for use in the fit. The resulting changes in yield and mass are used as the second category of systematic uncertainty.

Reconstruction of the decay sequences is the third source of systematic uncertainty. To evaluate the reliability of the MC determination of $`\pi ^0`$ efficiency and momentum calibration, we use control samples of $`K_S\pi ^0\pi ^0`$ and $`\tau \pi ^0X`$. On this basis, we assign a systematic uncertainty of $`\pm 5`$% in $`\pi ^0`$ reconstruction efficiency and a relative $`\pm 1`$% in $`\pi ^0`$ momentum bias. Similar studies for $`\gamma `$ reconstruction reveal a systematic uncertainty of $`\pm 3`$% in $`\gamma `$ reconstruction efficiency and $`\pm 1`$% in energy bias. Uncertainties in the $`D_s^+`$ and $`D_s^{}(2112)^+`$ masses, taken from world averages Hagiwara et al. (2002), also contribute to the systematic uncertainty.

The resulting total systematic uncertainty in the $`D_{sJ}(2458)^+`$ \[$`D_{sJ}^{}(2317)^+`$\] mass is $`\pm 1.0`$ \[$`\pm 0.8`$$`\mathrm{MeV}/c^2`$.

Using the yields from our fit and correcting for efficiency, we estimate the relative production rate

$$R=\frac{\sigma (D_{sJ}(2458)^+)(D_{sJ}(2458)^+D_s^{}(2112)^+\pi ^0)}{\sigma (D_{sJ}^{}(2317)^+)(D_{sJ}^{}(2317)^+D_s^+\pi ^0)}$$
(3)
to be $`0.25\pm 0.03(\text{stat.})\pm 0.03(\text{syst.})`$, requiring $`p^{}>3.5`$ $`\mathrm{GeV}/c`$ for both states. We also estimate, at 95% C.L.,

$$\frac{(D_{sJ}(2458)^+D_{sJ}^{}(2317)^+\gamma )}{(D_{sJ}(2458)^+D_s^{}(2112)^+\pi ^0)}<0.22.$$
(4)
The observed rms width of the $`D_{sJ}(2458)^+`$ is consistent with detector resolution, as determined by Monte Carlo studies. We conclude that the intrinsic width of the $`D_{sJ}(2458)^+`$ is small ($`\mathrm{\Gamma }10`$ $`\mathrm{MeV}/c^2`$).

The mass of the $`D_{sJ}(2458)^+`$ lies above $`DK`$ and below $`D^{}K`$ thresholds. The narrow width and the isospin-violating decay to $`D_s^{}(2112)^+\pi ^0`$ indicate that decay to $`DK`$ is forbidden and suggest an unnatural spin-parity assignment for the state. Belle has observed the decay $`D_{sJ}(2458)^+D_s^+\gamma `$ in production from both $`c\overline{c}`$ continuum Mikami et al. (2003) and $`B`$ decay Krokovny et al. (2003). Such a decay rules out $`J=0`$ and favors a $`1^+`$ interpretation. Decay distributions studied by Belle further support $`J=1`$ for $`D_{sJ}(2458)^+`$ and also $`J^P=0^+`$ for $`D_{sJ}^{}(2317)^+`$. The apparent absence of the decay $`D_{sJ}(2458)^+D_{sJ}^{}(2317)^+\gamma `$ may indicate that the electromagnetic decay mechanism cannot compete with $`D_s^{}(2112)^+\pi ^0`$, which may be a strong, but isospin-violating, process resulting from $`\eta `$-$`\pi ^0`$ mixing, as discussed by Cho and Wise Cho and Wise (1994).

Our measurement of the $`D_{sJ}(2458)^+`$ mass ($`2458.0\pm 1.4`$ $`\mathrm{MeV}/c^2`$, with combined statistical and systematic uncertainties) agrees with that obtained by Belle ($`2456.5\pm 1.7`$ $`\mathrm{MeV}/c^2`$Mikami et al. (2003), but is two standard deviations smaller than that obtained by CLEO ($`2463.1\pm 2.1`$ $`\mathrm{MeV}/c^2`$Besson et al. (2003). We obtain a relative yield ($`R=0.25\pm 0.04`$) which agrees with that of Belle ($`0.26\pm 0.08`$). Both values are somewhat smaller than that reported by CLEO ($`0.44\pm 0.13`$). Our reanalysis of the $`D_{sJ}^{}(2317)^+D_s^+\pi ^0`$ sample to account for background from the $`D_{sJ}(2458)^+`$ gives a mass of $`2317.3\pm 0.4(\text{stat.})\pm 0.8(\text{syst.})`$ $`\mathrm{MeV}/c^2`$, which remains consistent with results from CLEO Besson et al. (2003) and Belle Mikami et al. (2003).

In summary, in 91 $`\mathrm{fb}^1`$ of data collected from the BABAR experiment, we have observed a narrow state that decays to $`D_s^+\pi ^0\gamma `$ with a mass of $`2458.0\pm 1.0(\text{stat.})\pm 1.0(\text{syst.})`$ $`\mathrm{MeV}/c^2`$. The only significant $`D_s^+\pi ^0\gamma `$ decay mode we observe is through $`D_s^{}(2112)^+\pi ^0`$. We measure a mass and yield relative to the $`D_{sJ}^{}(2317)^+`$ similar to those measured by Belle though smaller than those reported by CLEO. The observed width is compatible with our mass resolution. After including the influence of this state, our new measurement of the $`D_{sJ}^{}(2317)^+`$ mass is $`2317.3\pm 0.4(\text{stat.})\pm 0.8(\text{syst.})`$ $`\mathrm{MeV}/c^2`$.

###### Acknowledgements.

We are grateful for the extraordinary contributions of our PEP-II colleagues in achieving the excellent luminosity and machine conditions that have made this work possible. The success of this project also relies critically on the expertise and dedication of the computing organizations that support BABAR. The collaborating institutions wish to thank SLAC for its support and the kind hospitality extended to them. This work is supported by the US Department of Energy and National Science Foundation, the Natural Sciences and Engineering Research Council (Canada), Institute of High Energy Physics (China), the Commissariat à l’Energie Atomique and Institut National de Physique Nucléaire et de Physique des Particules (France), the Bundesministerium für Bildung und Forschung and Deutsche Forschungsgemeinschaft (Germany), the Istituto Nazionale di Fisica Nucleare (Italy), the Foundation for Fundamental Research on Matter (The Netherlands), the Research Council of Norway, the Ministry of Science and Technology of the Russian Federation, and the Particle Physics and Astronomy Research Council (United Kingdom). Individuals have received support from the A. P. Sloan Foundation, the Research Corporation, and the Alexander von Humboldt Foundation.


text:
# 1 Introduction

## 1 Introduction

Computer simulation of opinion dynamics (consensus models) (Axelrod 1997; Deffuant 2000; Deffuant 2002; Weisbuch 2002; Hegselmann 2002; Hegselmann 2004; Krause 1997; Sznajd-Weron 2000; Stauffer 2000; Stauffer 2002; Galam 1990; Galam 1997; Stauffer 2003) is an important part of sociophysics (Weidlich 2000; Moss de Oliveira 1999; Schweitzer 2003). One checks if, starting from a random distribution of opinions (Monte Carlo method), one ends up with a consensus or a diversity of final opinions. The simulated people (”agents”) are located on lattices, on scale-free networks (Albert 2002; Barabási 2002), or form a purely topological structure where everybody can be connected with everybody. For the particular case of the consensus model of Deffuant et al (Deffuant 2000; Deffuant 2002; Weisbuch 2002), it was shown that on a Barabási-Albert (BA) network (Stauffer 2004) (see also Weisbuch 2004) the number $`S`$ of different surviving opinions (if no complete consensus was achieved) was an extensive quantity, i.e. it varied proportional to the number $`N`$ of agents, while it is intensive (independent of $`N`$ for large $`N`$) when everybody can be connected to everybody (Ben-Naim 2003). The literature on Barabási-Albert networks contains many comparisons with reality, e.g. for the computer networks of the Internet.

The motivation of the present work is two-fold: We want to have an unambiguous criterion whether two opinions agree or disagree, and thus use discrete instead of continuous variables for the opinions, section 2. Then we want to make the model more realistic by introducing noise representing events outside the opinion dynamics of Deffuant et al, by using in section 3 a more realistic network (Davidsen 2002; Holme 2002; Szabó 2003) with a higher clustering coefficient that the BA network, by allowing for advertising through mass media, and by taking into account more than one layer in order to implement an age structure; the last two effects are dealt with in section 4. An appendix gives the basic Fortran program.

## 2 Basic model and noise

Instead of allowing for the opinions any real number between 0 and 1, we take them as discrete numbers $`q=1,2,\mathrm{}Q`$, as in the Sznajd model (Sznajd-Weron 2000; Stauffer 2000; Stauffer 2002). Now it is well defined if two opinions differ or agree, while for real numbers it depends on the accuracy of the simulation. At first, only people differing by $`\pm 1`$ in their opinion can convince each other (bounded confidence (Hegselmann 2002; Hegselmann 2004; Krause 1997, Deffuant 2000; Deffuant 2002; Weisbuch 2002)); thus $`1/Q`$ corresponds to the confidence interval of the previous models. If two agents with opinions differing by one unit talk to each other, randomly one of them takes the opinion of the other (Axelrod 1997). We put agents on a directed Barabási-Albert (Albert 2002; Barabási 2002; Stauffer 2004) network, starting with $`m=3`$ agents connected with each other and with themselves; thereafter, $`N`$ agents are added, each of which selects $`m`$ pre-existing agents to be connected with. These randomly selected old agents are not regarded as connected to the new agent, i.e. the connections are directed. In this BA network the number of agents having $`k`$ neighbours is known to be proportional to $`1/k^3`$. The size of the network is the number $`N`$ of agents on it, i.e. the population.

First we construct the network, then we start the opinion dynamics from a random distribution of opinions. For each iteration we go through all agents in the order in which they were added to the network, and each selects randomly one of the $`m`$ agents it had chosen before to be connected with. The simulation stops if no agent changed opinion during one iteration. (About the same results are obtained from random instead of regular updating, provided we stop if for ten consecutive iterations no opinion changed.)

Figure 1 shows that for large $`N`$ the number $`S`$ of surviving final opinions roughly equals $`Q`$ for not too small $`Q`$; for $`Q=2`$, on the other hand, nearly always a complete consensus $`S=1`$ was found. (We averaged over 1000 samples except for $`N=10^5`$ when only 100 samples were used.) If, however, $`Q`$ grows to values closer to $`N`$, then the finite size of the network is felt: $`S`$ is lower than $`Q`$ and approaches $`N+m`$, that means everybody keeps its own opinion and the simulation stops soon. A finite-size scaling formula

$$S=(Q1)f(Q/N);f(x0)=1,f(x\mathrm{})=1/x$$

fits reasonably the same data, Fig.2, except for small $`Q`$. Thus, the number $`S`$ of final opinions is an extensive quantity if $`Q`$ is varied as $`N`$, and it is an intensive quantity if $`Q`$ is kept constant when $`N\mathrm{}`$. In this sense the new results are in between the intensive $`S`$ of (Deffuant 2000; Deffuant 2002; Weisbuch 2002; Ben-Naim 2003) and the extensive $`S`$ of (Stauffer 2004). (Also in contrast to (Stauffer 2004), the histograms for the number of people sharing the same opinion show a single peak. As in (Stauffer 2004) people may share the same opinion even if they are disconnected.) Fig.3 shows that for a fixed ratio $`N/Q=100`$ the cluster numbers are extensive.

Roughly the same results, Fig.4, are obtained if noise is introduced to simulate outside information in the opinion dynamics. Thus at each iteration, every agent after the above dynamics shifts the opinion by +1 with probability 1/4, by $`1`$ with probability 1/4, and keeps it unchanged with probability 1/2. (However, the opinion cannot leave the interval from 1 to $`Q`$. The simulation stops if without this noise no opinion would have changed.) Thus the model is robust against this noise. Of course, now no complete consensus is found, not even at $`Q=2`$. (If we follow (Deffuant 2000; Deffuant 2002; Weisbuch 2002) and allow everybody to interact with everybody, the results with noise are nearly the same as in Fig.4.)

The present model gets closer to the original Deffuant model if we introduce another free parameter $`L`$ such that two people convince each other if their opinions do not differ by more than $`L`$ units; $`L=1`$ then is our previous discrete model. If $`Q`$ and $`L`$ both go to infinity at constant ratio $`d=L/Q`$, then this ratio is the $`d`$ of Deffuant et al. The parameter $`\mu `$ of (Deffuant 2000; Deffuant 2002; Weisbuch 2002) was taken as $`0.1^{1/2}`$. Fig. 5 shows the variation of the number of surviving opinions with $`L`$, at fixed $`Q=1000`$ and various $`N`$. No noise is used here since noise prevents convergence except for very small $`L`$. For $`Q=N=5L`$ the number of opinion clusters seems to vary proportional to $`N`$ for large $`N`$.

## 3 Triads

Although the Barabási-Albert (Albert 2002; Barabási 2002; Stauffer 2004) network has successfully explained the scale-free nature of many networks, a striking discrepancy between it and real networks is that the value of the clustering coefficient - which is the probability that two nearest neighbours of the same node are also mutual neighbours - predicted by the theoretical model decays very fast with the network size and for large systems is typically several orders of magnitude lower than found empirically (it vanishes in the thermodynamic limit, Fig.6a). In social networks (Davidsen 2002; Holme 2002; Szabó 2003), for instance, the clustering coefficient distribution $`C(k)`$ exhibits a power-law behaviour, $`C(k)k^\gamma `$, where $`k`$ is number of neighbours (degree or connectivity) of the node and $`\gamma 1`$ (everyone in the network knows each other).

Very recently, by adding a triad formation step on the Barabási-Albert prescription, this problem has been surmounted and scale-free models with high clustering coefficient have been investigated (Davidsen 2002; Holme 2002; Szabó 2003). The algorithm is defined as follows: First, each new node performs a preferential attachment step, i.e, it is attached randoml to one of the existing nodes with a probability proportional to its degree; then follows a triad formation step with a probability $`p_{tf}`$: the new node selects at random a node in the neighbourhood of the one linked to in the previous preferential attachment step. If all neighbours are already connected to it, then a preferential attachment step is performed (“friends of friends get friends”). In this model, the original Barabási-Albert network corresponds to the case of $`p_{tf}=0`$. It is expected that a nonzero $`P_{tf}`$ gives a finite nonzero clustering coefficient as $`N`$ is increased, and with the clustering coefficient going to zero when $`P_{tf}=0`$ (the BA scale-free network model).

Simulations of the discretized version of the Deffuant model on this network produces similar results (see Fig. 6) to those obtained using a Barabási-Albert network (Fig. 1). The same behaviour was found for any value of the probability $`p_{\mathrm{tf}}`$ to perform a triad formation step. Furthermore, the clustering coefficient (Fig. 6b) agrees with the one predicted in (Davidsen 2002; Holme 2002; Szabó 2003), with a nearly linear increase with probability $`p_{\mathrm{tf}}`$. (Other values of $`m`$ give qualitatively the same behaviour.)

For completeness we end with some results from the directed model of (Stauffer 2004) on the size distribution of opinion clusters, those for the undirected case are similar. Fig.7 shows one peak for small sizes, with cluster numbers increasing proportional to the system size, and another peak for large sizes of the order of the network size, with small size-independent cluster numbers.

## 4 Several layers

Following Schulze (Schulze 2004), we now put $`A`$ copies of the same directed Barabási-Albert network on top of each other, with $`N`$ agents on every layer. Each layer corresponds to a certain age cohort, with babies on the bottom (layer $`A`$), the oldest old on the top (layer 1), and intermediate ages in the $`A2`$ layers in between. Each person in every layer dies with probability $`p`$ at each iteration. In case of death, the younger people on that position in the network move one layer up, keeping their opinion, and the lowest layer is occupied by a newly born baby getting the opinion of the parent. Initially, all people of different ages on the same position in the network share the same opinion, but later each layer draws its random numbers independently for random sequential updating. Thus it is crucial that the same network appears in $`A`$ copies on top of each other, in order to have a unique identification of ancestors and offspring. We refrain from comparing the model with university teachers (full professors, associate professors, …, down to teaching assistants) waiting for the superior to leave the job and make it available to younger ones.

Our simulations averaged over 1000 samples and four combinations of the number $`N`$ of agents and the number $`Q`$ of possible opinions: $`(N,Q)`$ = (10,10), (100,10), (1000,10) and (10,1000). We used mortalities $`p=0.01`$ and 0.5, maximum age $`A`$ = 2, 3, 5, 10, and lengths $`L`$ of the confidence interval = 1,2, … 6 for $`Q=10`$, and $`L=10,20,\mathrm{}60`$ for $`Q=100`$. Figure 8 shows that only $`L/Q`$ is really important: If it is 0.6 or higher, all samples lead to a full consensus with only one opinion surviving in all layers together. For small $`L`$ the number $`S`$ of surviving opinions is seen to be slightly below 10, the usual maximum number of opinions. Even for $`L=1`$ and $`Q=100`$ or 1000 (not shown) $`S`$ remains of order ten if $`N=10`$; we have to increase $`N`$ to get larger $`S`$ for larger $`Q`$. Also with $`N=10000`$ and $`Q=10`$ the results (not shown) look like in Fig.8a. Thus the basic result is similar to the monolayer results in Sec. 2: For small $`L/Q`$ the number $`S`$ of surviving opinions is of the order of $`N`$ or $`Q`$, whatever is smaller. The number $`A`$ of layers and the mortality $`p`$ hardly influence the results.

In Fig.8a we let the simulations run, until during one iteration none of the $`A\times N`$ agents changed opinion. Figure 8b shows that the results are nearly the same if we follow (Schulze 2004) and stop the simulation when the baby layer remained unchanged during ten consecutive iterations; however, for $`L/Q0.6`$ instead only one surviving opinion we find on average between one and two.

Finally, we introduce advertising (Schulze 2003; Sznajd-Weron 2003) in favour of opinion $`q=1`$: with 3 percent probability, every agent at every iteration had its opinion reduced by one unit. Now $`S`$ is reduced and depends stronger on the various parameters, Fig.9. Again, for a confidence interval $`L/Q0.6`$ a complete consensus is found. Fig.10 shows for the same runs the success fraction: With what probability is the final opinion within one layer a consensus in the advertised opinion? For $`L/Q0.5`$ nearly complete success is seen.

## 5 Conclusions

By discretizing the opinions, the simulations of the Deffuant model could be simplified and made less ambiguous. Two limits are quite trivial: With many people and few opinions, nearly all opinions have some followers, and the number of final opinion clusters nearly agrees with the total number of opinions. In the opposite limit of many opinions for few people, nearly every person forms a separate opinion cluster. For the transition between these two limits, a simple scaling law is observed for the discretized opinions. At a fixed ratio of the number of people to the number of opinions, the number of final opinion clusters is extensive. Noise and a more realistic network with stronger clustering (Davidsen 2002; Holme 2002; Szabó 2003) do not change the results much in the discretized model. An ageing model with several layers representing different age groups gave results not much different from those of one single layer, also if advertising is included.

J. Hołyst suggested to include noise, and J. Kertész to use the network of (Davidsen 2002; Holme 2002; Szabó 2003). A.O. Sousa thanks a grant from Alexander von Humboldt Foundation.

## 6 Appendix

To facilitate others to continue this research, to state details unambiguously and to allow checks for possible errors, the main program used for section 2 is reproduced here. An electronic version is available from stauffer@thp.uni-koeln.de as deffuant14.f.

Loop 35 goes over various values of the confidence interval idis, loop 17 over nrun different samples. ibm is a random odd integer with 64 bits. All computations with derri... involve the Derrida-Flyvbjerg (Derrida 1986) parameter omitted from the text for simplicity and may be ignored by the reader. (The results for this parameter were similar to Weisbuch 2004.)

After initialization, loop 7 connects the initial core of the BA network, and loop 2 adds to it max sites i each of which builds m directed bonds to neighbours neighb(i,new), new=1,2,...m, selected randomly according to the BA rule with the help of the Kertész list of length L. After this network is built up, loop 5 initializes the opinions is randomly, and loop 9 makes the maxt iterations of the Deffuant process. In this process for each site j in loop 10 randomly a neighb(j,i) is selected, and if these two agree already or are too far from each other in their opinion, noting is done: goto 12. Otherwise the counter ichange is increased by one, thus giving the number of opinion pairs which have changed during this iteration. If both opinions differ only by one unit, one of them is selected randomly. If the opinion difference is greater than one (but not exceeding the confidence interval idis) then both opinions change by the same amount idiff according to the usual Deffuant rule. The lines adding noise to the opinion dynamics are commented out in this version. If loop 10 ends going through all sites j without having changed any opinion, the iterations stop: goto 11.

Now the analysis starts: nhist(i) counts how often opinion i is found in the final set of opinions is: loops 28 and 29. Loop 27 gives icount as the number of different surviving opinions, later averaged over many samples using ict to give the average number of different surviving opinions, i.e. the crucial quantity of this study. number is used for the binned size distribution ns of opinion clusters, and loops 33, 34, 31 can be ignored as mentioned above.

```
      parameter(nsites=10  ,m=3,iseed=4711,maxt=1000000,iq=1000
     1  ,nrun=1000,max=nsites+m, length=1+2*m*nsites+m*m)
      integer*8 ibm,mult
      dimension list(length),neighb(max,m),is(max),nhist(0:iq),
     1 number(max), ns(31),irand(0:3)
      data irand/0,0,-1,1/
      w=sqrt(0.1)
      print 100, iq, nsites, m, iseed, maxt, nrun, w
 100  format(’# directed, more confidence’, 6i9,f6.3)
      factor=(0.25d0/2147483648.0d0)/2147483648.0d0
      facto2=factor*2*iq
      mult=13**7
      mult=mult*13**6
      ibm=2*iseed-1
      do 35 idis=100,900,100
      derrisum=0.0
      do 16 i=1,31
 16     ns(i)=0
      ict=0
      do 17 irun=1,nrun
        do 29 n=1,max
 29       number(n)=0
        do 7 i=1,m
          do 7 nn=1,m
            neighb(i,nn)=nn
 7          list((i-1)*m+nn)=nn
        L=m*m
c       All m initial sites are connected with each other and themselves
        do 1 i=m+1,max
          do 2 new=1,m
 4          ibm=ibm*16807
            j=1+(ibm*factor+0.5)*L
            if(j.le.0.or.j.gt.L) goto 4
            j=list(j)
            list(L+new)=j
            list(L+m+new)=i
 2          neighb(i,new)=j
 1        L=L+2*m
c     end of network and neighbourhood construction, start of opinion change

      n=max
      do 5 i=1,n
        ibm=ibm*16807
 5      is(i)=1+iabs(ibm)*facto2
c     print *, is
      do 9 iter=1,maxt
       ichange=0
       do 10 j=1,n
 6      ibm=ibm*16807
        i=1+(ibm*factor+0.5)*m
        if(i.le.0.or.i.gt.m) goto 6
        i=neighb(j,i)
        if(is(i).eq.is(j) .or. iabs(is(i)-is(j)).gt.idis) goto 12
        ichange=ichange+1
        if(iabs(is(i)-is(j)).eq.1) then
          ibm=ibm*16807
          if(ibm.lt.0) then
            is(i)=is(j)
          else
            is(j)=is(i)
          end if
        else
          idiff=isign(ifix(0.5+w*iabs((is(i)-is(j)))),is(i)-is(j))
          is(j)=is(j)+idiff
          is(i)=is(i)-idiff
        endif
c10     print *, iter, is(j)
 12     continue
c       ibm=ibm*mult
c       index=ishft(ibm,-62)
c       noise
c       is(j)=min0(iq,max0(irand(index)+is(j),1))
 10     continue
c      if(iter.eq.(iter/1000 )*1000 ) print *, iter,ichange
       if(ichange.eq.0)  goto 11
 9    continue
      print *, ’not converged’
 11   continue
      do 28 i=0,iq
 28     nhist(i)=0
      do 25 i=1,n
        j=is(i)
 25     nhist(j)=nhist(j)+1
c     print *, iter, nhist
      icount=0
      do 27 i=1,iq
        if(nhist(i).gt.0) icount=icount+1
        if(icount.gt.0) number(icount)=number(icount)+nhist(i)
 27   continue
      ict=ict+icount
c     print *, irun,icount,iter
      icount=0
      do 33 i=1,max
 33     icount=icount+number(i)
      derrida=0.0
      fact=1.0/icount**2
      do 34 i=1,max
 34     derrida=derrida+fact*number(i)**2
      do 31 i=1,max
       if(number(i).eq.0) goto 31
       ibin=1+alog(float(number(i)))/0.69315
       ns(ibin)=ns(ibin)+1
 31   continue
c     print *, irun,icount,iter,iq,derrida
      derrisum=derrisum+derrida
 17   continue
      derrida=derrisum/nrun
c     do 32 i=1,31
c32   if(ns(i).gt.0) print *, 2**(i-1), ns(i)
      call flush(6)
 35   print *, idis, ict*1.0/nrun, derrida
      stop
      end

```

References

ALBERT, R. and Barabási, A.L. (2002), “Statistical mechanics of Complex networks”, Rev. Mod. Phys. 74, pp. 47-97.

AXELROD, R. (1997), “The Dissemination of Culture: A Model with Local Convergence and Global Polarization”, J. Conflict Resolut. 41, pp. 203-226; The Complexity of Cooperation: Agent-Based Models of Competition and Collaboration, Princeton University Press, Princeton NJ.

BARABÁSI, A.L. (2002), Linked: The New Science of Networks, Perseus Books Group, Cambridge MA.

BEN-NAIM, E., Krapivsky, P. and Redner, S. (2003), “Bifurcations and Patterns in Compromise Processes”, Physica D 183, pp. 190-204.

DAVIDSEN, J., Ebel, H. and Bornholdt, S. (2002), “Emergence of a small world from local interactions: Modeling acquaintance networks”, Phys.Rev.Letters 88, pp. 128701-1-128701-4.

DEFFUANT, G., Amblard, F., Weisbuch G. and Faure, T. (2002), “How can extremism prevail? A study based on the relative agreement interaction model”, Journal of Artificial Societies and Social Simulation 5 (4), paper 1 (http://jasss.soc.surrey.ac.uk/5/4/1.html).

DEFFUANT, G., Neau, D. Amblard, F. and Weisbuch, G. (2000), “Mixing beliefs among interacting agents. Advances in Complex Systems”, Adv. Complex Syst. 3, pp. 87-98.

DERRIDA, B. and Flyvbjerg, H. (1986), “Multivalley structure in Kauffman’s model: analogy with spin glasses”, J. Phys. A 19, pp. L1003-L1008.

GALAM, S. (1990), “Social paradoxes of majority rule voting and renormalisation group”, J. Stat. Phys. 61, pp. 943-951 (1990).

GALAM, S. (1997), “Rational Group Decision Making: a random field Ising model at T=0”, Physica A 238, pp. 66-80.

HEGSELMANN, R. and Krause, U. (2002), “Opinion Dynamics and Bounded Confidence Models, Analysis and Simulation”, Journal of Artificial Societies and Social Simulation 5 (3), paper 2 (http://jasss.soc.surrey.ac.uk/5/3/2.html).

HEGSELMANN, R. and Krause, U. (2004), “ Opinion Dynamics Driven by Various Ways of Averaging”, Physica A (in press) = http://pe.uni-bayreuth.de/?coid=18.

HOLME, P. and Kim, B.J. (2002), “Growing scale-free networks with tunable clustering”, Phys. Rev. E 65, pp. 026107-1-026107-4.

KRAUSE, U. (1997), “Soziale Dynamiken mit vielen Interakteuren. Eine Problemskizze”. In Krause, U. and Stöckler, M. (Eds.), Modellierung und Simulation von Dynamiken mit vielen interagierenden Akteuren, pp. 37-51, Bremen University, Bremen.

MOSS DE OLIVEIRA, S., de Oliveira, P.M.C. and Stauffer, D. (1999): Evolution, Money, War and Computers, Teubner, Stuttgart and Leipzig.

SCHULZE, C. (2003), “Advertising in the Sznajd Marketing Model”, Int. J. Mod. Phys. C 14, pp. 95-98.

SCHULZE, C. (2004), “Advertising, consensus, and ageing in multilayer Sznajd model”, Int. J. Mod. Phys. C 15, No. 4 (in press) = cond-mat/0312342.

SCHWEITZER, F. (2003), Brownian Agents and Active Particles, Springer, Berlin.

STAUFFER, D. (2002), “Monte Carlo simulations of Sznajd models”, Journal of Artificial Societies and Social Simulation 5 (1), paper 4 (http://jasss.soc.surrey.ac.uk/5/1/contents.html).

STAUFFER, D. (2003), “Sociophysics Simulations”, Computing in Science and Engineering 5 (3), pp. 71-75; “How to Convince Others? Monte Carlo Simulations of the Sznajd Model” , in AIP Conference Proceedings, The Monte Carlo Method on the Physical Sciences: Celebrating the 50th Anniversary of the Metropolis Algorithm, edited by Gubernatis, J.E., 690 (1), pp. 147-155 = cond-mat/0307133.

STAUFFER, D., Sousa, A.O and Moss de Oliveira, S. (2000), “Generalization to Square Lattice of Sznajd Sociophysics Model”, Int. J. Mod. Phys. C 11, pp. 1239-1245.

STAUFFER, D. and Meyer-Ortmanns, H. (2004), “Simulation of Consensus Model of Deffuant et al on a Barabasi-Albert Network”, Int. J. Mod. Phys. C 15 (2) (in press) = cond-mat/0308231.

SZABÓ, G., Alava, M. and Kertész, J. (2003), “Structural transitions in scale-free networks”, Phys. Rev. E 67, pp. 056102-1-056102-1.

SZNAJD-WERON, K. and Sznajd, J. (2000), “Opinion Evolution in Closed Community”, Int. J. Mod. Phys. C 11, 1157.

SZNAJD-WERON, K. and Weron, R.(2003), “How effective is advertising in duopoly markets?”, Physica A 324, pp. 437-444.

WEIDLICH, W. (2000), Sociodynamics; A Systematic Approach to Mathematical Modelling in the Social Sciences, Harwood Academic Publishers, Amsterdam.

WEISBUCH, G., Deffuant, G., Amblard, F. and Nadal, J.-P. (2002), “Meet, Discuss, and Segregate!”, Complexity 7, pp. 55-63.

WEISBUCH, G. (2004), “Bounded confidence and social networks”, Eur. Phys. J. B, Special Issue: Application of Complex Networks in Biological Information and Physical Systems (in press) = cond-mat/0311279 .


text:
# Large scale inhomogeneity of inertial particles in turbulent flow

## Abstract

Preferential concentration of inertial particles in turbulent flow is studied by high resolution direct numerical simulations of two-dimensional turbulence. The formation of network-like regions of high particle density, characterized by a length scale which depends on the Stokes number of inertial particles, is observed. At smaller scales, the size of empty regions appears to be distributed according to a scaling law.

The transport of inertial particles in fluids displays properties typical of compressible motion even in incompressible flows. This is a consequence of the difference of density between particles and fluid. The most peculiar effect is the spontaneous generation of inhomogeneity out of an initially homogeneous distribution. The clustering of inertial particles has important physical applications, from rain generation FFS02 , to pollutant distribution and combustion DWM96 . Starting from the first examples in laminar flow Maxey87 , it is now demonstrated both numerically SE91 ; WM93 ; HC01  and experimentally FKE94 ; ACHL02  that also in fully developed turbulence there is a tendency of inertial particles to form small scale clusters. The parameter characterizing the effect of inertia is the Stokes number $`\mathrm{St}`$, defined as the ratio between the particle viscous response time $`\tau _\mathrm{s}`$ and a characteristic time of the flow $`\tau _\mathrm{v}`$. In the limit $`\mathrm{St}0`$ inertial particles recover the motion of fluid particles and no clusterization is expected. In the opposite limit $`\mathrm{St}\mathrm{}`$ particles become less and less influenced by the velocity field. The most interesting situation is observed for intermediate values of $`\mathrm{St}`$ where strong clusterization is observed FKE94 ; HC01 .

In the case of a smooth velocity field the Eulerian characteristic time $`\tau _\mathrm{v}`$ is a well defined quantity as it can be identified with the inverse Lyapunov exponent $`\tau _\mathrm{v}=\lambda _1^1`$ of fluid trajectories. In this case some general theoretical predictions are possible EKR96 ; BFF01  such as the exponential growth of high order concentration moments. Detailed numerical simulations in a chaotic random flow have shown maximal clusterization (measured in terms of the dimension of the Lagrangian attractor) for a value $`\mathrm{St}0.1`$ Bec03 .

In the case of turbulent flow, where the velocity field is not smooth, a simple scaling argument suggests that maximal compressibility effects are produced by the smallest, dissipative scales BFF01 . Nevertheless, for sufficiently large values of $`\mathrm{St}`$, the particle response time introduces a characteristic scale in the inertial range which breaks the scale invariance of the velocity field and produces, as we will see, large scales inhomogeneity in particle distribution.

The motion of a spherical particle in an incompressible flow, when the size of the particle is so small that the surrounding flow can be approximated by a Stokes flow, is governed by the set of equations MR87

$$\{\begin{array}{c}\dot{𝒙}=𝒗\hfill \\ \\ \dot{𝒗}=\frac{1}{\tau _\mathrm{s}}\left[𝒗𝒖(𝒙(t),t)\right]+\beta \frac{\mathrm{d}}{\mathrm{d}t}𝒖(𝒙(t),t)\hfill \end{array}$$
(1)
where $`𝒗`$ represents the Lagrangian velocity of the particle, $`\beta =3\rho _0/(\rho _0+2\rho )`$ and $`\rho `$ and $`\rho _0`$ are the density of particle and fluid respectively. In (1) $`𝒖(𝒙,t)`$ represents the velocity field whose evolution is given by Navier-Stokes equations

$$\frac{𝒖}{t}+𝒖\mathbf{}𝒖=\mathbf{}p+\nu \mathrm{\Delta }𝒖+𝒇$$
(2)
In what follows, we will consider the limit heavy particles such that $`\beta 0`$. In this limit it is easy to show that the Lagrangian velocity possesses a compressible part BFF01 : expanding (1) to first order in $`\tau _\mathrm{s}`$ and using $`\mathbf{}𝒖=0`$, one obtains, from (2)

$$\mathbf{}𝒗\tau _\mathrm{s}\mathbf{}\left(𝒖\mathbf{}𝒖\right)0$$
(3)
From (3) it is possible to give a dimensional estimation of the relative importance of the compressible part for a turbulent velocity field with scaling exponent $`h`$, $`\delta _{\mathrm{}}uU(\mathrm{}/L)^h`$, $`L`$ and $`U`$ being a characteristic large scale and velocity. The scaling exponent for the compressible component of $`𝒗`$ is $`\delta _{\mathrm{}}vU(\mathrm{}/L)^{2h1}`$ and thus the relative compressibility scales as $`(\mathrm{}/L)^{h1}`$, i.e. reaches the maximum value at the viscous scale BFF01 . Nevertheless, we will see that the presence of inertial range of scales in the turbulent flow generates large scale structures in the particle distribution at large St.

In this letter we address the problem of transport of heavy inertial particles in fully developed two-dimensional turbulence in the inverse cascade regime. High resolution direct numerical simulations of two-dimensional Navier-Stokes equation with white in time, random forcing $`𝒇`$ localized at small scales $`\mathrm{}_\mathrm{f}`$ have been performed. As customary, a friction term $`\alpha 𝒖`$ is added to (2) in order to extract energy from the system at the friction scale $`\mathrm{}_{\mathrm{fr}}\epsilon ^{1/2}\alpha ^{3/2}`$. The intermediate scales $`\mathrm{}_\mathrm{f}\mathrm{}\mathrm{}_{\mathrm{fr}}`$ define the inertial range in which Kolmogorov scaling $`\delta _{\mathrm{}}uU(\mathrm{}/L)^{1/3}`$ is observed BCV00 . Lagrangian tracers are placed at random with initial zero velocity and integrated according to (1) with a given $`\tau _\mathrm{s}`$. After a scratch run long several $`\tau _\mathrm{s}`$, Lagrangian statistics is accumulated for typically some tens of $`\tau _\mathrm{s}`$. Stokes time is made dimensionless by rescaling with the Lagrangian Lyapunov exponent of fluid particles, $`\mathrm{St}\lambda _1\tau _\mathrm{s}`$. Fig. 1 shows typical distributions of inertial tracers in stationary conditions at different values of $`\mathrm{St}`$, obtained starting from the same initial random distribution. One observes in both cases strong inhomogeneous distributions with empty “holes”, in the second case on much larger scales.

Maximum compressibility effects are expected at small scales and can be described by the Lyapunov spectrum for inertial particles. We recall that for a generic dynamical system the sum of the Lyapunov exponents gives the exponential rate of expansion (or contraction) of the hypervolume in phase space. In our case, from (1) we have $`_{i=1,4}\lambda _i=2/\tau _\mathrm{s}`$, thus volumes are contracted at a constant rate. Let us observe that when $`\tau _\mathrm{s}\mathrm{}`$, phase space contraction rate vanishes, and we thus expect less clusterization. As a consequence of the structure of (1) we find that two Lyapunov exponents are close to $`1/\tau _\mathrm{s}`$, representing the rate of adjustment of Lagrangian velocity to the Eulerian one. The first Lyapunov exponent is found positive, as the trajectories are chaotic and the second, negative, determines the dimension of the attractor according to the definition of Lyapunov dimension ott

$$d_\mathrm{L}=K+\frac{\underset{i=1}{\overset{K}{}}\lambda _i}{|\lambda _{K+1}|}$$
(4)
where $`K`$ is defined as the largest integer such that $`_{i=1}^K\lambda _i0`$. In the inset of Fig. 2 we show the dependence of the first Lyapunov exponent on Stokes number. As $`\mathrm{St}`$ increases, $`\lambda _1(\mathrm{St})`$ decreases monotonically from the neutral value $`\lambda _1(0)0.72`$. The behavior of Lyapunov dimension is also shown in Fig. 2. In the limit $`\mathrm{St}0`$, particles become neutral and thus one recovers the homogeneous distribution with $`d_\mathrm{L}=2`$. The presence of a minimum around $`\mathrm{St}0.1`$ was already discussed in the case of smooth flow Bec03  and indicates a value for which compressibility effects are maximum. We remark that the curve $`d_L(St)`$ of Fig. 2 is almost identical to the one obtained in synthetic smooth flow Bec03 , and is thus insensitive to the presence of the hierarchy of scales typical of a turbulent flow.

In the turbulent scenario, instead of becoming more homogeneous, for increasing $`\mathrm{St}`$ the inertial particle distribution develops structures on larger scales, as is evident in Fig. 1. The most evident dishomogeneities are related to the presence of empty regions (“holes”) on different scales. We have thus studied the hole statistics at varying Stokes number.

We have performed a coarse graining of the system by dividing it into small boxes forming the sites of a square lattice and counting the number of particles contained in each small box. From this coarse grained density we have computed the probability density function of holes, defined as connected regions of empty boxes.

Probability density functions (pdf) of hole areas are shown in Fig. 3 for different Stokes numbers. The hole distributions follow a power law with an exponent $`1.8\pm 0.2`$ up to an exponential cutoff at a scale $`A_{\mathrm{St}}`$ which moves to larger sizes with $`\mathrm{St}`$, as shown in the inset of Fig. 3. At variance with the smooth flow case, in which particle density recovers homogeneity for $`\mathrm{St}>0.1`$ Bec03 , in the turbulent case inhomogeneities are pushed to larger scales when increasing $`\mathrm{St}`$.

The hole area pdf is independent on the number of particles used in the simulations. This is a non trivial property, reflecting the fact that inertial particles cluster on network-like structures where a clear-cut distinction between empty regions and particle-rich regions can be observed. We have also verified that the choice of the small coarse graining scale does not modify the hole area pdf at larger scales.

These results support the following physical picture. According to (1), heavy particles filter out velocity fluctuations on timescales shorter than $`\tau _S`$. We may thus expect that inhomogeneity does appear at a scale $`\mathrm{}_S`$ whose characteristic time is of order of $`\tau _S`$. For a velocity field with a roughness exponent $`h1`$, this would lead to $`A_{\mathrm{St}}\mathrm{St}^{2/(1h)}`$, i.e. to the increase of hole areas with $`\mathrm{St}`$, in qualitative agreement with Fig. 3. Clusterization of inertial particles in turbulence appears to be a selfsimilar process up the cutoff scale $`A_{\mathrm{St}}`$. Indeed the pdfs of hole areas collapse when areas are rescaled with the cutoff area, as shown in Fig. 4.

A natural measure of the degree of clusterization of inertial particles which has been considered in the literature SE91 ; HC01 ; FKE94  is the deviation of the particle density distribution from Poisson law, $`D(R)=(\sigma _R\sigma _R^{(P)})/\lambda _R`$, where $`\lambda _R`$ and $`\sigma _R`$ represent the mean and standard deviation of particle distribution coarse grained at scale $`R`$ and $`\sigma _R^{(P)}`$ is the corresponding standard deviation for a Poisson distribution. The maximum of the deviation $`D(R)`$, which depends on $`St`$ (see Fig. 5) is then used for defining the typical scale of clusters, $`R_{St}`$. We have found that this definition of $`D(R)`$ is strongly sensitive to particle statistics, nevertheless, the increase of $`R_{St}`$ on $`St`$ (inset of Fig. 5) seems to be a robust result.

The presence of structures in the inertial particle distribution is often attributed SE91  to the fact that heavy particles are expelled from vortical regions. Although structures are related to the presence of many active scales in the turbulent flow, one should recall that in 2D turbulence, as a consequence of the direct vorticity cascade, vorticity is concentrated at the small scales BCV00  and no large scale coherent structures appear. Holes emerge as a result of the delayed dynamics (1), which filters the scales of the underlying turbulent flow characterized by times of the order of the Stokes time $`\tau _\mathrm{s}`$.

As discussed before, clustering occurs also in synthetic flows where a hierarchy of time scales is absent, just as a consequence of the dissipative character of the motion Bec03 . However, it appears from our simulations that to fully understand the geometry of inertial particle distribution in a turbulent flow the presence of structures characterized by a large set of time scales cannot be ignored.

We conclude that the geometry of inertial particle clusters in developed turbulence is controlled both by the dissipative effective dynamics of the particle motion at small scales, and by the tendency of inertial particles to filter the active scales characterized by times of the order of the characteristic relaxation time of the particles. A full understanding of the geometry of particle clusters in developed turbulence is particularly relevant for several applications, such as coalescence processes or chemical reactions. The reaction rate of two chemical species is a function of their concentration. When particles of different species are transported by the same turbulent flow, their local concentrations are not independent and the presence of large scale correlations in the particle distribution can in principle influence the reaction velocity.

This work was supported by MIUR-Cofin2001 contract 2001023848. We acknowledge the allocation of computer resources from INFM “Progetto Calcolo Parallelo”.