## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Natural time analysis of critical phenomena

Contributed by Seiya Uyeda, May 21, 2011 (sent for review May 10, 2011)

## Abstract

A quantity exists by which one can identify the approach of a dynamical system to the state of criticality, which is hard to identify otherwise. This quantity is the variance of natural time *χ*, where and *p*_{k} is the normalized energy released during the *k*th event of which the natural time is defined as *χ*_{k} = *k*/*N* and *N* stands for the total number of events. Then we show that *κ*_{1} becomes equal to 0.070 at the critical state for a variety of dynamical systems. This holds for criticality models such as 2D Ising and the Bak–Tang–Wiesenfeld sandpile, which is the standard example of self-organized criticality. This condition of *κ*_{1} = 0.070 holds for experimental results of critical phenomena such as growth of rice piles, seismic electric signals, and the subsequent seismicity before the associated main shock.

- short-term earthquake prediction
- dynamic exponent
- fractional Gaussian noise
- fractional Brownian motion
- Burridge–Knopoff “train” model

It has been shown that some unique dynamic features hidden behind can be derived from the time series of complex systems, if we analyze them in terms of natural time *χ* (1–3). For a time series comprising *N* events, we define an index for the occurrence of the *k*th event by *χ*_{k} = *k*/*N*, which we term natural time. In doing so, we ignore the time intervals between consecutive events, but preserve their order and energy *Q*_{k}. We, then, study the evolution of the pair (*χ*_{k}, *Q*_{k}) by using the normalized power spectrum [1]defined by , where *ω* stands for the angular natural frequency and [2]is the normalized energy for the *k*th event. In the time-series analysis using natural time, the behavior of Π(*ω*) at *ω* close to zero is studied for capturing the dynamic evolution, because all the moments of the distribution of the *p*_{k} can be estimated from Φ(*ω*) at *ω* → 0 (see ref. 4, p. 499). For this purpose, a quantity *κ*_{1} is defined from the Taylor expansion Π(*ω*) = 1 - *κ*_{1}*ω*^{2} + *κ*_{2}*ω*^{4} + ⋯, where . We found that this quantity, the variance of natural time *χ*_{k}, is a key parameter for the distribution of energy within the natural time interval (0,1]. Note that *χ*_{k} is “rescaled” as natural time changes to *χ*_{k} = *k*/(*N* + 1) together with rescaling upon the occurrence of any additional event. It has been demonstrated that this analysis enables recognition of the complex dynamic system under study entering the critical stage (1–3). This occurs when the variance *κ*_{1} converges to 0.070. Originally the condition *κ*_{1} = 0.070 for the approach to criticality was theoretically derived for the seismic electric signals (SES) (1, 2), which are transient low frequency (≤ 1 Hz) electric signals that have been repeatedly observed before earthquakes (3, 5, 6). The experimental data showed that *κ*_{1} obtained from SES activities in Greece and Japan attain the value 0.070 (1–3, 7–10). The emission of SES was attributed to a phase transition of second order. It was also shown empirically that the same condition *κ*_{1} = 0.070 holds for other time series, including turbulence (8) and seismicity preceding main shocks (3, 7–11). Moreover, it has been found empirically that main shocks occur, in terms of the conventional time, a few days up to one week after the condition *κ*_{1} = 0.070 was attained for the seismicity subsequent to SES activity (1, 3, 7–10) supporting the concept that seismicity is a critical phenomenon (e.g., refs. 12 and 13 and references therein). Despite these numerous successes (e.g., see refs. 3, 14), however, the condition *κ*_{1} = 0.070 for criticality has remained a major stumbling block for wider acceptance, because the validity of this condition has not been theoretically demonstrated for the cases other than the SES activities and the Burridge–Knopoff “train” model for earthquakes (15). In order to remedy this situation, in this paper, we will try to identify the origin of the validity of the *κ*_{1} = 0.070 condition for a wider range of critical systems.

## Explanation of *κ*_{1} = 0.070 for Critical Systems

We deal with time series of signals emitted from complex dynamical systems. When the system is in thermodynamic equilibrium, it should produce stationary time series with *p*_{k} independent of *χ*_{k}. The situation is different when the system is not in equilibrium. When the system approaches the critical state, clusters of the new phase are formed by enhanced fluctuation and their size increases as does the correlation length (16–18). But this happens not instantly because long-range correlations develop gradually leading to the dynamic phase transition of the second order (17). Thus, the time series emitted in such a nonequilibrium process will be nonstationary and *p*_{k} will be not any more independent of *χ*_{k}.

Using , which is the distribution corresponding to *p*_{k}, the normalized power spectrum of Eq. **1** can be rewritten as [3]

Taylor expansion of Eq. **3** around *ω* → 0 leads to the value [4]

We are interested in *p*(*χ*) of a dynamic system approaching criticality, which characterizes the way energy is released during the evolution of the dynamic transition. The newly formed phase may in general be coupled with an existing external field and the interaction energy is expected to be proportional to the linear dimension of the newly forming phase and hence to the correlation length *ξ* (for example, once “charge” is conserved, in the new phase we may only have charge separation leading to dipole moment). Thus, we expect *p*(*χ*) ∝ *ξ*. Because of the critical slowing down when approaching dynamic transition, the time-dependent correlation length *ξ* becomes as expressed by *ξ* ∼ *t*^{1/z}, where *z* is the dynamic critical exponent. Here we assume this relation holds also for the natural time domain as *ξ* ∝ *χ*^{1/z}, which leads to [5]where *N*_{c} is a normalization constant to make . In fact, Eq. **5** is plausible from the definition of *p*_{k}; i.e., it represents the normalized energy emitted during the *k*th event, and the energy at criticality has a power law distribution.

The above would be applied to the case of earthquakes because the state just before a big earthquake may be characterized by a long chain of dislocations or faults just like the long chain of aligned spins in the Ising model in the critical state. Substituting Eq. **5** into Eq. **4**, we obtain [6]

Substituting the value of the dynamic critical exponent *z* for various universality classes of critical systems (19), we can obtain the values of *κ*_{1} depicted in Fig. 1. Note that for most universality classes, *z* varies in the region from *z* = 2 to *z* = 2.4, and thus (see Fig. 1) the value of *κ*_{1} obtained by Eq. **6** is in the range of 0.068 to 0.071, including the 2D Ising model, which is qualitatively similar to the process of SES emission [early and most recent Monte Carlo calculations showed *z* = 2.165 (see ref. 20) and *z* = 2.154 (see ref. 21) leading through Eq. **6** to *κ*_{1} ≈ 0.070]. These results seem to justify the substitution of *t* by *χ* because the time *t* used for the computation of *z* in Monte Carlo steps (MCS) is the internal clock of the system, which can be considered as equivalent to the natural time.

## The Case of a 2D Ising Model

We now show numerically that in a 2D Ising system quenched from a high temperature to a temperature close to (but below) the critical temperature the value of *κ*_{1} approaches 0.070. The calculations are carried out as follows: A 2D Ising system of linear size *L*, with periodic boundary conditions, is prepared in a high temperature state and then quenched to a temperature just below *T*_{c}. Considering that the Hamiltonian for the interaction between two spins is given by , where *s*_{i} = ± 1 and *J* > 0 stands for the coupling constant between *s*_{i} and *s*_{j}, the evolution of the magnetization per spin is simulated by the standard Metropolis algorithm (22) and studied as a function of the number *k* of MCS. The number *k* is set to zero when the system is quenched and increases by 1 at each MCS following the standard Metropolis algorithm. For the purpose of the present simulation, *k* runs from *k* = 1 to 10^{4} MCS. Fig. 2*A* depicts the ensemble average of |*M*_{k}|, which corresponds to the correlation length *ξ*, obtained from 10^{3} replicas for various sizes *L* = 100, 200, 400, and 1,000. It is observed in Fig. 2*A* that, due to the well-known phenomenon of critical slowing down (22), systems of larger linear size need larger number of MCS to finally reach the equilibrium magnetization. We now present in Fig. 2*B* a log–log plot of the values shown in Fig. 2*A*. This clearly reveals that, practically independent of *L*, the dynamics of is a power law of the form with *z* very close to *z* = 2.165, which is the value given in ref. 20 for the dynamic exponent for the 2D Ising model (see the cyan straight line in Fig. 2*B*). This dynamic model was then analyzed in natural time by setting *Q*_{k} = |*M*_{k}|. Fig. 2*C* depicts the results for *κ*_{1} as a function of the number *k* of MCS that followed the instantaneous quench. It is clear that *κ*_{1} ≈ 0.070.

## The Case of Self-Organized Criticality.

Natural time analysis has been applied to the experimental dataset of a self-organized criticality (SOC) system such as rice pile (23) as well as to the time series obtained numerically from a SOC model based on the Burridge–Knopoff train model for earthquakes (15). In both cases it has been shown that *κ*_{1} converges to 0.070 at the onset of the critical stage. Here, we present the theoretical results obtained from the natural time analysis of time series of the avalanches in the archetypal system that exhibits SOC, e.g., the Bak–Tang–Wiesenfeld (BTW) sandpile model (24). The BTW model is a multiparticle dynamical system wherein the dynamics cannot be reduced to a few degrees of freedom (24, 25). After some initial transient period, the system settles down to a steady state described by power law distributions as in the case of the second-order phase transitions.

Let us consider the BTW model on a *D*-dimensional hypercubic lattice of linear size *L* in which the number of sand particles at each lattice site is given by the integer variables *z*_{i}≥0. We perturb the system by adding a sand particle at a site *i* that means *z*_{i} → *z*_{i} + 1. When *z*_{i} equals the value 2*D* and the site becomes unstable, the site relaxes (topples). At that time, its *z*_{i} value decreases by 2*D*, and the number of sand particles of its 2*D* nearest neighbors (nn) increases by one: [7][8]If the neighboring sites become unstable, an avalanche may proceed. This avalanche stops when all sites are stable again. An avalanche is characterized by its size *s* (the total number of topplings). According to the basic hypothesis of BTW (24), in the SOC state the probability distribution of the avalanche sizes exhibits power law behavior: [9]where *τ* is the size exponent.

In order to proceed to numerical simulations, we study a deterministic version of the BTW sandpile model (25), where the random site seeding is replaced by seeding at the central site. Wiesenfeld et al. (25) showed that the system for *D* = 2 also evolves into a SOC state. We found that the natural time analysis of the series of avalanches with initial condition *z*_{i} = 0 leads to the *κ*_{1} values plotted in Fig. 3 for *D* = 2 to 7.

The *κ*_{1} values for various *D* plotted in Fig. 3 fluctuate around the following values: *κ*_{1} = 0.056, 0.064, 0.069, 0.071, 0.073, and 0.075 for *D* = 2 to 7, respectively. Interestingly, these values are given by Eq. **6** for *z* = *D*/2. This result can be understood on the following grounds: Because an avalanche occurs every time when 2*D* sand particles are fed into the central site, the number of avalanches is equal to that of particles fed *n* divided by 2*D*. Natural time increases by 1/*N* when an avalanche occurs; therefore we have [10]where *N* is the total number of avalanches and the brackets [.] denote the integer part.

According to Dhar (26), formulas **7** and **8** lead to the expected number of toppling *G*_{ij} at site *j* upon adding a particle at site *i*: [11]where *r*_{ij} is the distance between the sites *i* and *j*. Because we deal with a centrally fed sandpile, the total expected number of topplings is found by integrating formula **11** in the hypersphere of radius *l* of the sandpile: [12]

With regard to *l*, recent mathematical studies (27) have shown that the linear dimension of the formed sandpile grows as [13]

Inserting Eq. **10** and formula **13** into formula **12**, we obtain , which explains why the observed *κ*_{1} values are compatible with those obtained from Eq. **6** with *z* = *D*/2.

Our results in Fig. 3 indicate that *κ*_{1} ≈ 0.07 within 10% for *D*≥3 but not for *D* = 2. Note that Ktitarev et al. (28) showed that avalanches for *D* = 2 deviate from power law behavior.

## Fractional Gaussian Noises and Fractional Brownian Motions

It has been shown (7) that when the self-similarity index deduced from the detrended fluctuation analysis (DFA exponent *α*) of the time series for fractional Brownian motions and fractional Gaussian noises approaches unity, reflecting the infinitely long-range temporal correlations, the quantity *κ*_{1} approaches 0.070. It may be added here that the presence of long-range temporal correlations in SES activities has been established because they also lead to the values of *α* close to unity (29–31).

## Conclusion

Based on the concept of natural time, an explanation has been proposed for the experimental fact that becomes equal to 0.070 when a variety of dynamical systems enter the critical stage.

The case of the Ising model was studied here because it is widely known and in addition it is qualitatively similar to the generation mechanism of SES (5, 6, 32). The only difference is that the factor that brings about the critical state is the temperature in the case of the Ising model, whereas it is the stress in the focal region in the case of SES.

Results exhibiting similar behavior were presented for other critical systems including SOC on which unpredictability of earthquakes has been erroneously claimed. The fact that *κ*_{1} becomes equal to 0.070 for the seismicity before the main shock can be used for earthquake prediction purposes. Actually, the occurrence time of a main shock is specified in advance by analyzing in natural time the seismicity subsequent to the initiation of the SES activity (1, 3, 7–10, 15). This analysis identifies the time when the seismicity approaches the critical state, i.e., when the condition *κ*_{1} = 0.070 is attained. The main shock was found empirically to follow usually within a few days up to one week. This has been successfully applied to several major earthquakes in Greece, including the strongest one (*M*_{w}6.9) during the last 28 years (14).

## Acknowledgments

We thank D. Turcotte and H. Kanamori for valuable comments on an earlier version and reviewers, in particular H. Ezawa, on the present version. This research was partially supported by the Ministry of Education, Science, Sports and Culture, Grant-in-Aid for Young Scientists (B), 21710180, 2009–2011 (to M.K.), and Scientific Research (C), 20510171, 2008–2010, and 23510218, 2011 (to S.U. and M.K.), and Observation and Research Program for Prediction of Earthquakes and Volcanic Eruptions, 2010–2011 (to S.U. and M.K.).

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: suyeda{at}st.rim.or.jp.

Author contributions: P.V. and N.V.S. designed research; P.V., N.V.S., E.S.S., S.U., and M.K. performed research; P.V., N.V.S., and E.S.S. contributed new reagents/analytic tools; N.V.S. and E.S.S. analyzed data; and P.V., N.V.S., E.S.S., S.U., and M.K. wrote the paper.

The authors declare no conflict of interest.

## References

- ↵
- Varotsos PA,
- Sarlis NV,
- Skordas ES

- ↵
- ↵
- Varotsos P,
- Sarlis N,
- Skordas ES

- ↵
- Feller W

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Keilis-Borok VI,
- Soloviev AA

- ↵
- ↵
- Uyeda S,
- Kamogawa M

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Landau DP,
- Binder K

- ↵
- ↵
- ↵
- ↵
- ↵
- Fey A,
- Levine L,
- Peres Y

- ↵
- ↵
- ↵
- ↵
- ↵
- Varotsos P,
- Alexopoulos K

- ↵

## Citation Manager Formats

## Sign up for Article Alerts

## Jump to section

## You May Also be Interested in

_{2}conditions, by the late 21st century, increasing temperatures could lead to reduced snowpack, drier summers, and increased fire risk, independent of changes in winter precipitation.

### More Articles of This Classification

### Physical Sciences

### Related Content

- No related articles found.