Review of the Natural Time Analysis Method and Its Applications

Abstract

A new concept of time, termed natural time, was introduced in 2001. This new concept reveals unique dynamic features hidden behind time-series originating from complex systems. In particular, it was shown that the analysis of natural time enables the study of the dynamical evolution of a complex system and identifies when the system enters a critical stage. Hence, natural time plays a key role in predicting impending catastrophic events in general. Several such examples were published in a monograph in 2011, while more recent applications were compiled in the chapters of a new monograph that appeared in 2023. Here, we summarize the application of natural time analysis in various complex systems, and we review the most recent findings of natural time analysis that were not included in the previously published monographs. Specifically, we present examples of data analysis in this new time domain across diverse fields, including condensed-matter physics, geophysics, earthquakes, volcanology, atmospheric sciences, cardiology, engineering, and economics.

1. Introduction

Natural time analysis (NTA) was introduced in 2001 [1,2,3]. It originally focused on the discrimination of low-frequency ($f\le 1$ Hz) electric signals that are statistically significant [4] precursors observed [5] before earthquakes (EQs) in Greece [6,7], China [8], Japan [9,10,11,12,13], and Mexico [14], which are called seismic electric signals (SESs) (cf. for their physical properties and method of observation, see [15,16,17,18,19,20]), from similar looking signals of human-made origin [21,22]. SESs [17] provide information on the magnitude, M, of the impending EQ and the probable epicentral area, as well as an estimation of its occurrence time. In these original publications, the fact that NTA may reveal that SESs exhibited critical dynamics has also been used for the identification of criticality in seismicity within a probable epicentral area, which allowed for the shortening of the occurrence time-window of the impending EQ.

Natural time, which is a new concept of time, is inherently discrete and, hence, drastically different from conventional (clock) time [23,24,25]. NTA unveils new dynamical features in time-series originating from complex systems that may remain hidden if a study is conducted in the framework of conventional time. In general, NTA facilitates the surveillance of an evolving complex system and recognizes when the system approaches a critical stage [24,25,26,27].

NTA, as such, has found useful applications in various fields. The scope of this paper is to review the applications of NTA in various complex systems. Initially, we focus on briefly presenting applications that appeared upon the publication of the first monograph on NTA that appeared [24] in 2011. Subsequently, we include here the most recent NTA applications (see, e.g., Section 3.3.1, Section 3.5, Section 3.6, Section 3.7 and Section 3.8) that were not reviewed in the second monograph [25] published in 2023.

This review is organized as follows: In the next section, we will summarize the methodology, which is the background of NTA, while in Section 3, various applications in diverse systems will be discussed. Specifically, in Section 3.1, applications in the general field of condensed matter physics will be reviewed. These include the electric signals that precede a rupture in Section 3.1.1, the magnetic flux avalanches observed in high-T_c superconductors in Section 3.1.2, the avalanches observed in 3D rice piles in Section 3.1.3, the acoustic emission observed in granular materials in Section 3.1.4, the fluctuations of electrical resistance before a fracture in Section 3.1.5, the electromagnetic emissions before a fracture in LiF in Section 3.1.6, and self-organized critical systems in general in Section 3.1.7. In Section 3.2, applications in geophysics are presented that include fracture-induced electromagnetic emissions before EQs in Section 3.2.1, the subionospheric propagation anomalies before EQs in Section 3.2.2, and the crustal deformation before the 2016 Kumamoto EQs in Section 3.2.4. Section 3.3 reviews NTA applications to EQs, including the most recent applications, in Section 3.3.1, that involve the combination of NTA with nonextensive statistical mechanics (NESM) [28], which improves the estimation of the occurrence time of an impending EQ. Applications to volcanology are reviewed in Section 3.4. NTA applications to atmospheric sciences are the subject of Section 3.5. There, we briefly discuss the very recent introduction [29,30] of the notion of a natural rank for studying extreme natural phenomena related to rain and clouds. We also discuss applications on ozone-hole dynamics over Antartica in Section 3.5.1 and the precursory signals of major El Niño events in Section 3.5.2. Section 3.6 reviews the applications of NTA to cardiology, while Section 3.7 and Section 3.8 discuss applications to engineering and economics, respectively. Finally, in Section 4, we compile the most important conclusions.

2. Methodology: Natural Time Analysis

In a time-series comprising N events (cf. a temporal point pattern [31], as suggested in [26]), like the occurrence of the six EQs shown in Figure 1, the natural time, ${\chi }_k$, corresponding to the k-th event is defined as ${\chi }_k=k/N$ [1,2,3]. This way, we discard the interoccurrence times between successive events but retain their energy, $Q_k$, together with their order, k. In NTA, the pair $({\chi }_k,p_k)$, where $p_k=Q_k/{\sum }_{n=1}^N{Q}_n$ is the normalized energy, is attributed to each event, k. In order to proceed, we consider the normalized power spectrum $\Pi \left(\omega \right)\equiv {|\Phi \left(\omega \right)|}^2$, where $\Phi \left(\omega \right)={\sum }_{k=1}^N{p}_{k}exp\left(i\omega {\chi }_k\right)$. In that $p_k$ can be considered (see, e.g., [25]) as probabilities for the occurrence of each event, $\Phi \left(\omega \right)$ with $\omega \in \mathcal{R}$ is the characteristic function for $p_k$ [32]. As is well known, the study of $\Phi \left(\omega \right)$ as $\omega \to 0$ reveals important properties of the distribution of $p_k$, such as, for example, the fact that the moments of this distribution equal the derivatives $d^m\Phi \left(\omega \right)/d{\omega }^m$ ($m>0$) at $\omega \to 0$. To this end, a quantity, ${\kappa }_1$, was obtained from the Maclaurin series

$$\Pi \left(\omega \right)=1-{\kappa }_1{\omega }^2+{\kappa }_2{\omega }^4+\dots ,$$

(1)

where

$${\kappa }_1=\langle {\chi }^2\rangle -{\langle \chi \rangle }^2=\sum _{k=1}^N{p}_k{\left({\chi }_k\right)}^2-{\left(\sum _{k=1}^N{p}_k{\chi }_k\right)}^2.$$

(2)

It has been shown that

$${\kappa }_1\approx 0.070$$

(3)

at criticality for various dynamical systems [1,24,25,33,34], e.g., see Figure 1 of [34]. In NTA, ${\kappa }_1$ can be used to identify the approach to a critical point. The NTA of EQ catalogs has shown [35] that ${\kappa }_1$ can be used as an order parameter (OP) for seismicity. The value diminishes abruptly when a strong EQ takes place. The statistics of ${\kappa }_1$ fluctuations share properties that can be observed in well-known critical systems [35,36,37]. In particular, the scaled distribution of ${\kappa }_1$ collapses on that of other critical systems exhibiting, for several orders of magnitude, a characteristic exponential tail; e.g., see [38].

It has been proven that, with NTA, we can monitor the ${\kappa }_1$ fluctuations within an EQ catalog, which constitutes a great advantage [24,25]. Let us assume an excerpt of an EQ catalog with W consecutive events. We want to evaluate all possible ${\kappa }_1$ values resulting from subexcerpts of consecutive 6 to W EQs (cf. initially, Sarlis et al. [39] used 40 EQs as the maximum subexcerpt length). To this end, we apply Equation (2) for all subexcerpts of $N=$ 6, 7, and up to W consecutive EQs using as $Q_k$ the quantity ${10}^{1.5{\left(M_W\right)}_k}$, where ${\left(M_W\right)}_k$ is the moment magnitude [40] corresponding to the k-th EQ, which is proportional to the energy emitted during this EQ [41]. We calculate, for this set of ${\kappa }_1$ values, the mean value and the standard deviation, $\mu \left({\kappa }_1\right)$ and $\sigma \left({\kappa }_1\right)$, respectively. Subsequently, we can estimate the variability in ${\kappa }_1$:

$${\beta }_W={\displaystyle \frac{\sigma \left({\kappa }_1\right)}{\mu \left({\kappa }_1\right)}}.$$

(4)

The variability, ${\beta }_W$, quantifies the fluctuations, ${\kappa }_1$, within this excerpt of the EQ catalog. The definition of Equation (4) is reminiscent [42] of the Ginzburg criterion [43], the usefulness of which has been discussed by Holliday et al. [44] within the context of EQs. Finally, by sliding the window of W consecutive EQs, event by event, in the EQ catalog, we can pursue how ${\beta }_W$ varies in conventional time. Each ${\beta }_W$ obtained is associated with the occurrence time of the EQ that followed the last one in the window studied.

The entropy, S, in natural time [1,22], is given by

$$S=\langle \chi ln\chi \rangle -\langle \chi \rangle ln\langle \chi \rangle =\sum _{k=1}^N{p}_k{\displaystyle \frac{k}{N}}ln\left({\displaystyle \frac{k}{N}}\right)-\left(\sum _{k=1}^N{p}_k{\displaystyle \frac{k}{N}}\right)ln\left(\sum _{k=1}^N{p}_k{\displaystyle \frac{k}{N}}\right).$$

(5)

When it comes to “uniform” distribution [21], $Q_k$ are positive random variables that are independent and identically distributed (IID), and S becomes [45] $S_u\equiv {\displaystyle \frac{ln2}{2}}-{\displaystyle \frac{1}{4}}\approx 0.0966$. The application [24,25,46] of time reversal $\widehat{T}$, i.e., $\widehat{T}{p}_k=p_{N-k+1}$, in Equation (5) leads to

$$S_{-}=\sum _{k=1}^N{p}_{N-k+1}{\displaystyle \frac{k}{N}}ln\left({\displaystyle \frac{k}{N}}\right)-\left(\sum _{k=1}^N{p}_{N-k+1}{\displaystyle \frac{k}{N}}\right)ln\left(\sum _{k=1}^N{p}_{N-k+1}{\displaystyle \frac{k}{N}}\right).$$

(6)

The quantity $S_{-}$ is different from S. This leads to an entropy change

$$\Delta S\equiv S-S_{-}$$

(7)

upon time reversal. Hence, S is time-reversal-asymmetric [24,46,47]. Furthermore, S is dynamic entropy; it is not static like Shannon entropy [45,47,48], and it exhibits [46] positivity, concavity, and experimental stability [49,50].

The entropy change, $\Delta S$, is important in order to determine the approach of the system to a dynamic phase transition [25,26,27]. Using Equations (5) and (6) for $N=W$, the quantities S and $S_{-}$ are calculated in the natural time domain for a window length of W consecutive events. When this window slides, event by event, through the whole time-series of events, the corresponding time-series of S and $S_{-}$ can be deduced. Then, Equation (7) is used for the estimation of $\Delta S$. This results in a time-series of consecutive $\Delta S_W$ values where each $\Delta S_W$ is associated with the occurrence time of the event that followed the last one of the window studied.

Evaluating the standard deviation $\sigma (\Delta S_W)$ of the latter time-series, we can estimate [24,51] the complexity measure

$${\Lambda }_W={\displaystyle \frac{\sigma (\Delta S_W)}{\sigma (\Delta S_{100})}}.$$

(8)

Importantly, a different choice of the window length ($W=100$) in the denominator does change the values obtained, but the qualitative results and their physical interpretation concerning the time evolution of ${\Lambda }_W$ remain unaltered [52]. As is evident from Equation (8), ${\Lambda }_W$ measures how the statistics of the $\Delta S_W$ time-series changes when increasing the scale from $W=100$ to another scale, W. The latter has major significance when we study a dynamically evolving complex system; see, e.g., Chapter 3 of [24], and see also [51,53,54,55,56,57]. Notably, based on the entropy, S, various additional complexity measures that have been reviewed in [58] can be defined.

3. Applications

This section is devoted to the discussion of the applications of NTA in various complex systems. These applications are arranged according to the corresponding descipline they serve.

3.1. Condensed-Matter Physics

3.1.1. Electric Signals That Precede Rupture

In the 1980s, a physical model termed the “pressure stimulated polarization currents” (PSPC) model [59] (see also [16,60]) was suggested, and it revealed that transient electric signals precede EQs [17]; see also [5,15,16,17,61]. This is schematically shown in Figure 2. In the ionic constituents of rocks in the Earth, electric dipoles always exist [59] due to lattice defects; see, e.g., [62]. These electric dipoles initially have random orientations (see Figure 2c) in the future focal region of an EQ. For example, there begins a gradual stress, $\sigma$, increase as a result of the excess stress accumulation. This is what we call stage A; see Figure 2a. The electric dipoles exhibit a cooperative orientation (see Figure 2e) when this gradually increasing stress reaches a critical value (${\sigma }_{cr}$) (since ”cooperativity” is a hallmark of criticality [63]). This results in the release of a transient electric current termed an SES (see Figure 2b) of current density j. We call this stage B. As noticed by Uyeda et al. [64], the PSPC model does not require any sudden stress variation like microfracturing [65], and hence, it is unique among all electromagnetic EQ precursor-generation models [17]. SES constitutes a critical phenomenon in view of the “cooperativity” described above. Intense SES preceding major EQs are accompanied by detectable variations [66] in the Earth’s magnetic field, predominantly involving the vertical component [17,67].

A multitude of SESs emitted consecutively within a short time is called SES activity [15]. A study of the SES properties may reveal the epicentral location and magnitude of the forthcoming EQ; see, e.g., [4]. The EQ prediction method based on SES is broadly called [68,69,70,71]. The VAN method is named according to the initials of the scientists Varotsos, Alexopoulos, and Nomicos who invented it. After long experimentation [17] during the 1980s and 1990s, a network of SES measuring stations was established in Greece; see Figure 3. It is hereafter called the VAN network.

The NTA of SES activities and other similar-looking signals of human-made origin that are usually recorded in the VAN network has revealed [1,2,3,21,22,72] that only SESs are critical. In particular, it was shown that, apart from satisfying Equation (3), SES activities also satisfy

$$S\le S_u$$

(9)

and

$$S_{-}\le S_u,$$

(10)

see, e.g., Table 4.6 at pp. 227 and 228 of [24].

3.1.2. Time-Series of Magnetic-Flux Avalanches Observed in High-T_c Superconductors

One of the first applications of NTA in condensed-matter physics is that [33] in time-series of magnetic-flux avalanches observed in high-T_c superconductors. The critical state in superconductors was suggested (see, e.g., [73]) to be a self-organized critical (SOC) [74,75] system. The similarity between a SOC, i.e., sand piles, and superconductors was first pointed out by de Gennes [76]. If type II superconductors are exposed to a slowly ramped-up, external magnetic field, magnetic vortices will enter the sample from its edges, but they may become pinned by crystallographic defects such as dislocations [59]. In this way, a flux gradient can be formed that is hardly stable, as is the slope of slowly built sandpiles. As a result, it may happen that a minor increase in the applied magnetic field can lead to significant flux rearrangements inside the sample; the latter are called flux avalanches [77,78,79].

The NTA of the magnetic flux, $\Delta {\Phi }_k$, avalanches observed in a thin film of YBa₂Cu₃O_7−x was based on detailed experiments made by Aegerter et al. [80] (see their Figure 2). The magnetic flux avalanche time-series in a typical such experiment were analyzed by Sarlis et al. [33] by setting $Q_k=\Delta {\Phi }_k$, see their Figure 1. This figure reveals that, after 140 avalanches, the ${\kappa }_1$ value is 0.070(5), thus verifying criticality, as mentioned in Section 2; see Equation (3). As far as the value of S is concerned, this results in $S\approx 0.085$, i.e., smaller than $S_u$.

Moreover, the above results were compared [33] with those obtained from the NTA of a generalized stochastic SOC model introduced by Carbonne and Stanley [81]. The probability density functions (PDFs) of ${\kappa }_1$ and S are shown in Figure 3 of [33], and they reveal that, for appropriate SOC-model parameters, and for $N=140$ avalanches, they reach their maximum values of ${\kappa }_1=0.070\left(10\right)$ and $S=0.080\left(10\right)$, respectively. Such values, which also point to criticality (see Equation (3)), are compatible with the ones obtained from the NTA of the penetration of magnetic flux into thin films of YBa₂Cu₃O_7−x.

3.1.3. Time-Series of Avalanches Observed in 3D Rice Piles

NTA has been applied to well-controlled experimental results obtained from three-dimensional (3D) piles of rice by Aegerter et al. [82,83]. In these works, Aegerter et al. studied the evolution of a 3D rice pile that was initially far from the critical state as it gradually approached SOC. They showed [83] that concepts [84] based on extremal dynamics can describe their experimental findings. Additionally, Aegerter et al. studied how avalanche sizes and the corresponding PDFs evolve, which constitutes a further test of extremal dynamics.

Sarlis et al. [33] employed NTA to the time-series of the size $\Delta V$ of avalanches, which was measured directly in one experiment by [82]. By considering $Q_k=\Delta V_k$, Sarlis et al. [33] found the results shown in their Figure 4. These results show that, actually, at later times (when the system starts to approach the critical state), the ${\kappa }_1$ value scatters in the region around 0.070(10), pointing to criticality, according to Equation (3). Moreover, it was found that $S\approx 0.070\left(10\right)<S_u$.

In order to further theoretically supplement these experimental results, Sarlis et al. [33] performed NTA in the avalanche time-series obtained from a simplification of the Burridge and Knopoff (BK) EQ model [85] suggested by de Sousa Vieira [86]. The latter model has been extensively studied [87] as it approached SOC. The results in Figure 5 of [33] show that, when the system approached criticality, ${\kappa }_1\approx 0.070$ with $S\le S_u$, which was compatible with the values obtained from the NTA of 3D rice piles.

3.1.4. Acoustic Emission Observed in Granular Materials

Tsuji and Katsuragi [88] performed the following detailed experiment: Glass beads of various grain diameters were poured into a cylindrical Plexiglas container. A steel sphere was then penetrated into the granular bed, and spheres of various radii were used. The penetration speed was kept fixed at three different values, and the acoustic emissions (AEs) due to the penetration were monitored using a fixed sensor buried in the granular bed.

Tsuji and Katsuragi [88] performed an NTA of their experimental results by considering the squared maximum amplitude of each AE event as $Q_k$. They examined the validity of Equation (3) and found that ${\kappa }_1$ fluctuates around 0.070 in some experimental setups, indicating criticality; for example (see their Figure 9a,b). However, this tendency was not universal.

Tsuji and Katsuragi [88] concluded that NTA showed that ${\kappa }_1$ is, in general, distributed at around 0.070–0.083, while Equation (3) can be established in the brittle-like regime of the granular material.

3.1.5. Fluctuations of Electrical Resistance Before Fracture

NTA was applied [89] to the study of electric resistance changes before the fracture of a cement mortar and Luserna stone samples (cf. the recorded AE were also analyzed in natural time). Assuming an undamaged specimen before the loading protocol starts, Niccolini et al. [89] managed to relate the electrical resistance of the damaged sample during the test to an active cross-sectional area that was reduced due to freshly formed microcracks. They further estimated the energy of each cracking step (event), which is necessary in NTA, by means of fracture mechanics. Equation (3), together with the relations (9) and (10), were checked in order to determine criticality. The NTA indicated that electrical resistance criticality systematically preceded that of AE for all specimens [89].

3.1.6. Electromagnetic Emissions Before Fracture in LiF

Lithium fluoride (LiF) is an insulator that possesses the largest reported bandgap [90], and it is expected [91] to remain transparent under stress exceeding 1000 GPa. Potirakis and Mastrogiannis [92] recorded electromagnetic emissions simultaneously with AE during fracture experiments in LiF crystals that were irradiated and non-irradiated by gamma rays.

An NTA of the electromagnetic emissions was carried out [92] by assuming various thresholds and estimating ${\kappa }_1$, S, and $S_{-}$ in a manner similar to that suggested by Varotsos et al. [93]. Criticality was identified by means of Equation (3) and relations (9) and (10) in the electromagnetic emissions before fracture for various thresholds. Potirakis and Mastrogiannis [92] also performed NTA for the recorded AE, finding critical behavior that was expected to arise during the fracture process and was observed before that of the electromagnetic emissions.

3.1.7. Self-Organized Critical Systems

As already mentioned in the previous sections, Section 3.1.2 and Section 3.1.3, the SOC models suggested by Carbonne and Stanley [81] and de Sousa Vieira [86] were studied via NTA that may have revealed their approaches to criticality on the basis of Equation (3). However, the archetypal example of SOC is the growing of a sand pile [74,75]. The centrally fed Bak–Tang–Wiesenfeld (BTW) sand-pile model [74] was studied in natural time by Varotsos et al. [34]. The results compiled in their Figure 3 reveal that, after initial fluctuations, and as the system gradually approached SOC, ${\kappa }_1\approx 0.070$ satisfied Equation (3) within 10% when the underlying dimension, D, of the sand pile was $D\in [3,7]$. It is noteworthy that the NTA results found in the centrally fed BTW model were recently experimentally observed [94] in a study of acoustic emission before the fracture of notched fiber-reinforced concrete specimens that were subjected to three-point bending.

Another model that constitutes an approximation of the BK model [85] and that simulates EQ behavior is the Olami—Feder–Christensen (OFC) model [95], which concerns cellular automaton. The OFC model introduced dissipation into the family of SOC systems, although its criticality was debated [96,97]. The OFC model’s SOC behavior can be destroyed through some small changes to the model’s rules [98,99,100,101], while it cannot account for certain aspects of the spatiotemporal clustering of seismicity [102]. Nonetheless, the OFC model is considered [103] closer to reality than other models concerning EQ predictability [104] or the Omori law [105,106]. Additionally, for appropriate dissipation values, the OFC model may lead to an avalanche-size PDF that is compatible [107] with the law of Gutenberg–Richter (GR). Hence, the OFC model is considered [108] to be exemplary for a supposedly SOC system describing seismicity. We note, however, that whether EQs are appropriately described or not through SOC, or whether different kinds of mechanisms are necessary (see, e.g., [109,110]), remains unsolved [56,105,106,111,112,113].

The NTA of the OFC model was presented in detail in Section 8.3 of [24]. There, it was shown that, as the system approaches SOC, the criticality condition of Equation (3) is observed; see Figures 8.5 and 8.6 of [24]. Moreover, before strong avalanches, a finite entropy change, $\Delta S$, in natural time was found upon the time reversal [24,114]. This unveils a breaking of time symmetry, thus revealing the predictability of the OFC model; see, e.g., [115]. In particular, Figure 8.12 of [24] shows that $\Delta S$ exhibits a clear minimum before large avalanches.

3.2. Geophysics

Apart from the case of SESs presented in Section 3.1.1, NTA has been applied to the identification of criticality in a variety of other EQ precursors.

3.2.1. Fracture-Induced Electromagnetic Emissions Before EQs

Potirakis et al. [116] presented fracture-induced (or fracto-) electromagnetic emissions that were recorded before strong earthquakes in Greece; see also [117,118]. NTA was applied [116,117,118,119,120] by employing a threshold, as in the case of long-duration SES activities [93] and Equation (3); moreover, conditions (9) and (10) were found to be valid for a wide range of threshold values.

3.2.2. Subionospheric Propagation Anomalies Before EQs

In order to detect subionospheric, very low-frequency (VLF) propagation anomalies preceding EQs, custom-designed VLF/LF receivers were installed [121] at various measuring stations. NTA was applied after a variety of phenomena (like magnetic storms, solar flares, typhoons, tsunamis, and volcanic eruptions; see, e.g., [122,123,124]) that may affect measurement quality were removed. The criteria of Equation (3) and relations (9) and (10) for criticality were examined. This analysis revealed [116] that a lower ionosphere reached criticality (probably associated with the Kumamoto EQs) 2 to 10 days prior to the 15 April 2016 $M_{W}7.0$ EQ; see also Section 4 of [121]. An NTA of the VLF propagation anomalies prior to the 30 October 2020 Aegean Sea $M_W$7.0 EQ also revealed [125] criticality that was observed within the last two weeks before the main shock.

3.2.3. Ultra Low-Frequency (ULF) Magnetic Field Variations Before EQs

For this application of NTA, geomagnetic measurements from stations in China, Japan, and Greece were used [116,118,126,127,128,129]. An NTA of ULF magnetic field variations revealed criticality—by means of Equation (3) and the relations (9) and (10)—before the 11 March 2011 Tohoku $M_W$9.0 EQ [127], the 12 May 2008 Sichuan $M_W$7.9 EQ [128], the 15 April 2016 Kumamoto $M_W$7.0 EQ [129], the 12 June 2017 Lesvos $M_W$6.3 EQ [118], and the 12 April 2013 Kobe $M_{W}5.8$ EQ [126].

3.2.4. Deformation Before the 2016 Kumamoto EQs

Yang et al. [130] studied the Global Navigation Satellite System (GNSS) deformation data of five stations prior to the 15 April 2016 Kumamoto $M_W$7.0 EQ. In order to support the atmospheric gravity wave (AGW) hypothesis of lithosphere–atmosphere–ionosphere coupling (LAIC), they attempted to recover the AGW components from these data. An NTA of the GNSS deformations before the 15 April 2016 Kumamoto $M_{W}7.0$ EQ revealed [130] that, in late March, three stations displayed criticality, i.e., Kumamoto with an epicentral distance of r = 10 km, Aso (with r = 38 km), and Sagara (with r = 57 km). In early April, criticality was reached at all five stations studied [130].

3.3. Earthquakes

The first applications of NTA to the physics of EQs were published in the early 2000s [1,3], and they focused on the estimation of the occurrence time of a strong EQ after SES detection. These studies focused on the two strong 1995 EQs at Grevena-Kozani $M_{W}6.6$ and Eratini-Egio $M_{W}6.5$, the 1997 Strofades $M_{W}6.6$ EQ in the Ionian Sea (see, e.g., [131]), and the 2001 Aegean Sea $M_{W}6.5$ EQ.

The continuation of this NTA of seismicity later revealed [35] that ${\kappa }_1$ may serve as an OP for EQs (which can be considered a critical phenomenon; see, e.g., [44] and references therein) that results in a universal curve both regionally [35] and globally [132,133], the details of which may reveal important properties of the EQ preparation process [134,135,136,137].

Interestingly, this OP, i.e., the quantity ${\kappa }_1$ can be used in the identification of correlations between EQ magnitudes [132,133,134,138,139,140]. More importantly, ${\kappa }_1$ exhibits characteristic fluctuations (quantified by ${\beta }_W$; see Section 2) before main shocks [39,42,56,135,141,142,143,144,145,146,147,148,149,150,151,152,153,154] that enable the estimation of the occurrence time and epicenter location of strong EQs based solely on EQ catalogs; for more details, see [25]. More recently, Varotsos et al. [155] demonstrated that the combination of ${\kappa }_1$, ${\beta }_W$, and ${\Lambda }_W$ reveals characteristic fluctuations in seismicity before strong EQs that are simultaneous with variations in the electric and magnetic fields of the Earth that precede [9,11,12,66,67,156] strong EQs.

In general, through an NTA of the seismicity after an SES activity in the area likely to suffer a strong EQ, we can determine the occurrence time of the forthcoming main shock with good accuracy, practically at around a few days to one week [1,3,10,13,35,138,157,158,159,160,161,162,163]. The NTA of seismicity has also been succesfully applied [119,120,164,165,166] to the identification of criticality around the epicenter of a future EQ. Other applications of natural time in seismicity include its providing the basis for predicting aftershock magnitudes [167,168] and providing assistance for a better description of seismic sequences [169] with encouraging results.

The most recent method for seismic risk estimation [153,163,170,171,172,173,174,175,176,177,178,179,180,181], which is termed EQ nowcasting, was introduced by Rundle et al. [182] based on the concept of natural time, as already mentioned in the introduction. This method has found several applications, including the study of induced seismicity [183,184], the temporal clustering of global EQs [185], epicentral information for a future EQ [152,153,154], etc. For example, the recent application of nowcasting to volcanic eruptions [186] is briefly described in the following section, Section 3.4.

3.3.1. Newer Primary Results

Very recently, the entropy change $\Delta S$ under time reversal, as well as its fluctuations in the NTA of seismicity, led [52,55,57,187] to additional information concerning the occurrence time of a strong EQ that may shorten [26,27,188] the expected time window. The latter advancement was based on the combination of NTA with NESM [28] (see Section 1), pioneered by Tsallis [189,190], who introduced the entropy $S_q$, and specifically in the determination of the $S_q$ of a fault by Posadas and Sotolongo-Costa [191]. A comparative study of $\Delta S_i(\equiv \Delta S_W$ for $W=i$) and $S_q$, calculated within the same window of i consecutive EQs after the criticality condition ${\kappa }_1=0.070$ has been fulfilled, reveals simultaneous variations before the strong EQ occurrence; e.g., see Figure 4. This figure presents, on the same time scale, the evolution of $\Delta S_i$ of NTA and the Tsallis entropy, $S_q$, for the seismicity reported by the Turkish Disaster and Emergency Management Authority (AFAD) within a 3° × 3° area around the 6 February 2023 Kahramanmaraş–Gaziantep $M_W$7.8 Pazarcik EQ epicenter. The red arrows indicate the simultaneous transient changes in $\Delta S_i$ and $S_q$ after the validity of the criticality condition ${\kappa }_1=0.070$. They start and end with the occurrence of two EQs of magnitudes 4.4 and 4.6, respectively. These two EQs occurred on 29 January 2023 at 16:12:39 UTC and on 3 February 2023 at 11:05:08 UTC with epicenters at 35.88 °N 35.78 °E and 37.21 °N 36.40 °E, respectively. The $M_W$7.8 Pazarcik EQ occurred almost 38 h later, on 6 February 2023 at 01:17:35 UTC; it was followed, within 9 h, by a $M_W$7.7 EQ with an epicenter almost 90 km north of that of the $M_W$7.8 EQ [192].

Apart from the aforementioned application to Turkey [188], this new method has been, so far, applied to retrospective studies of EQs in Japan [26] (2011 M9 Tohoku EQ), Mexico [27] (2017 M8.2 Chiapas and M7.1 Mexican flat-slab EQs), and California [27] (2019 Ridgecrest EQ). This is certainly a limitation, but we hope that, in the near future, we will have the opportunity to apply this method to real-time EQ prediction.

3.4. Volcanology

The 2018 Kīlauea volcano eruptive sequence (e.g., see [193]) was studied by means of nowcasting by Fildes et al. [186]. From mid-May to early August, 62 collapse events of seismic energy equal to that of a magnitude $M\ge 4.7$ EQ occurred. Fildes et al. [186] applied nowcasting to three time windows (17 May–2 August 2018, 29 May–2 August 2018, and 14 June–2 August 2018) to match the time periods when changes in the eruption sequence were identified by other means. The method suggested by Luginbuhl et al. [183] was applied, and the strongest agreement between the actual and nowcasted collapse events was achieved for the last time window, 14 June to 2 August 2018. This can be understood in the sense that, once a stable stress state is reached during an eruption, then successful nowcasting becomes possible.

3.5. Atmospheric Sciences

Very recently, Varotsos et al. [29], through an approach based on NTA, introduced the method of a natural rank for the use of remotely sensed geometric data on rain and clouds in studying extreme events. The same method was applied to highlight the relation between rain, clouds, and cosmic rays [30] (cf. for applications in nowcasting extreme cosmic ray events, see also [194]). These results are promising, and we hope this new method will soon be applied in real time. We also note that recent applications of natural time in the Earth-atmosphere system include the nowcasting of extreme low daily minimum air temperatures in the Arctic [195], the development of a novel nowcasting tool for the estimation of the average waiting time between extreme values of climate parameters [196], and the nowcasting of air-pollution episodes in modern megacities [197].

We now briefly summarize two additional NTA applications in atmospheric physics.

3.5.1. Ozone-Hole Dynamics over Antarctica

Varotsos and Tzanis [198] employed NTA in order to study the dynamics of the ozone hole over Antarctica using the maximum daily ozone hole area (MD-OHA). They considered each event to last approximately one year, and the value of the yearly maximum MD-OHA was considered to be $Q_k$. The values of S,$S_{-}$, and $\Delta S$ were studied [198] for various windows with lengths of 3 to 15 years, sliding year by year from 1977 to 2009; see their Figure 2. Interestingly, precursory changes were identified [198] that are consistent with the approach of the system to a dynamic phase-transition critical point approximately two years before September 2002, when the Antarctica ozone hole split [199].

3.5.2. Precursory Signals of Major El Niño Events

Varotsos et al. [200] studied the El Niño southern oscillation (ENSO) from January 1876 to November 2011 by means of an NTA. For this purpose, they used the Southern Oscillation Index (SOI) values, calculated by employing Troup’s formula [201], and these values are dimensionless. The SOI is related to the difference in the monthly surface air pressure between Tahiti and Darwin, and it measures the strength of the ENSO phenomenon [202]. Under the assumption that each event lasts approximately one month, the minimum monthly SOI values were considered in an NTA. For this purpose, the minimum monthly SOI value observed since 1876 (SOI_min = −38.8, measured in May 1896) was removed, and $Q_k$ = SOI_k − SOI_min, which are positive, were analyzed in natural time. Varotsos et al. [200] calculated the time-series of the entropy change $\Delta S_W$ under time reversal for moving window lengths of W months, and those authors examined their ability to predict weak or strong El Niño events one month in advance. They investigated various values of W and concluded, by means of the ROC [203], that the highest predictive skill would be calculated for $\Delta S_{20}$ to $\Delta S_{24}$. They also noticed that the quasi-biennial oscillation (QBO) in the zonal wind of the tropical stratosphere may provide a mechanism suggesting a time window of approximately 2 years (20 to 24 months). This link was further studied via an NTA in [204].

The method presented above was also employed [205] in predicting the strength of an ongoing El Niño event, the 2015–2016 El Niño event that has been reported to possibly become “one of the strongest on record” [206]. NTA suggested that this El Niño would become a “moderate to strong” or even “strong” one but not “one of the strongest on record”. Indeed, it was revealed that, upon the completion of the whole 2015–2016 El Niño event, the SOI values reported (see, e.g., http://www.bom.gov.au/climate/history/enso/ (accessed on 4 September 2023)), are compatible with a “moderate to strong” or even “strong” event, thus confirming what has been suggested well in advance by Varotsos et al. [205]. Very recently, Varotsos et al. [207] successfully applied this method, as well as the modified NTA method for nowcasting anomalies, to the 2023–2024 El Niño event, and they identified beforehand (August 2023) that this event would not be one of the strongest.

3.6. Cardiology

Heart dynamics have been studied in natural time using electrocardiograms (ECGs), photoplethysmography (PPG), or both. The results obtained are briefly summarized below.

3.6.1. NTA of ECGs

The turning points in ECGs labeled with the letters P, Q, R, S, and T are shown in Figure 5. In the NTA of ECGs, the various intervals that can be defined within a heartbeat, like QRS and QT, as well as interbeat intervals like RR (see Figure 5) or NN (which are [208] intervals between adjacent QRS complexes resulting from sinus-node depolarizations), were considered [45,47,48,51,209,210,211,212] $Q_k$.

In these studies, the entropy, S, in the NTA, as well as its fluctuations, $\delta S$ (defined as the standard deviation of the S time-series), were used to show the non-Markovianity of ECGs [45], as well as the distinction between healthy (H) humans and humans who experience sudden cardiac death (SCD) [45,48] based on either QT, RR, or QRS intervals. Varotsos et al. [48] introduced appropriate complexity measures to NTA using the RR, QRS, and QT intervals to achieve a classification of individuals with SCD, patients with heart disease, and truly H individuals; see also Chapter 9.3 of [24]. When this focus is solely on RR intervals, the importance of the entropy change, $\Delta S$, in natural time, as well as that of the corresponding fluctuations, was first identified by Varotsos et al. [47]. The results obtained [47] showed that $\Delta S_{13}$ may identify, in SCD individuals, the occurrence time of impending cardiac arrest. Moreover, when complexity measures based on $\Delta S$ are used, the SCD risk can be identified; see Figure 3 of [47]. An NTA of heart rate variability (HRV) was applied to various HRV models [209], as well as the study of HRV correlations during relaxing visualization [210]. The complexity measure ${\Lambda }_W$, mentioned in Section 2, was applied [51] to long-lasting—with a duration of several hours to 24 h—ECG recordings for the distinction of three classes: H, individuals with SCD, and patients suffering from congestive heart failure (CHF). The same complexity measure applied to ECGs was also used for comparisons in recent studies [211,212] that employed PPG and analyzed the data obtained in natural time.

3.6.2. NTA of Heart Dynamics Through PPG

The PPG technique has significantly simplified [211,213,214,215,216] the monitoring of HRV. A typical PPG signal received from the index finger of a young male individual is shown in Figure 6. The maximum of the PPG signal during a heartbeat is labeled as P, and an NTA is conducted under the assumption of $Q_k$ as the PP interval; i.e., the time period between adjacent P measures.

The results obtained [211] concerning the separation of H from SCD via 20 min PPG recordings are encouraging, and when Support Vector Machines (SVMs) [217,218,219,220,221,222] were used for the separation of these two classes, the sensitivity of detecting CHF may have reached 97.7%. Additionally, an appropriate PPG device for the NTA of HRV was constructed, and it enabled [212] remote sensing, thus simplifying follow-up studies. This method was incorporated [223] in an innovative e-health cloud-based system that is non-invasive and can be used easily in any setting. Employing modern methods of machine learning, Tatsis et al. [223] suggested classifiers for the discrimination of CHF patients from healthy controls, even from samples taken in a 10-min time window. At the current stage, larger studies are needed to further validate this novel e-health cloud-based system before its use in everyday clinical practice.

3.7. Engineering

The structural integrity and remaining life of constructions and/or critical structural elements upon the application of mechanical loads close to fracture are very important in both science and engineering [224]. AE is one of the most useful methods for the detection of pre-failure indicators [225]. Vallianatos et al. [226] employed an NTA for the study of AE before fracture in laboratory experiments on Etna basalt. The study revealed correlations in the magnitude of evolution of AE. Furthermore, this study revealed that the scaled distribution of the order parameter ${\kappa }_1$ showed [226] a characteristic feature that was similar in cases of seismicity and other equilibrium or non-equilibrium critical systems, including self-organized critical systems [35,132,133,135]; see also Section 3.1.7 and Section 3.3. This result pointed to the existence of a universal fracture behavior from the laboratory scale to the global scale.

Since fracture can be thought of as a phase change, an NTA of AE signals recorded when various materials are loaded to fracture can reveal information concerning the underlying dynamical system.

To this end, an NTA of AE was studied in fracture experiments for a variety of materials, like LiF [92], marble [227,228,229,229], steel [230], cement mortar [37,224], fiber-reinforced concrete specimens [94], rods of Luserna stone [89], and structures of technological interest [231].

The identification of criticality has been conducted by means of Equation (3) and the relations (9) and (10), revealing that an NTA of AE before fracture is potentially very significant to technological applications; see, e.g., [228,230].

Within this context, Loukidis et al. [37,224,232] examined the presence of criticality during the mechanical loading of cement mortar and Dionysos marble specimens until fracture. It was found [37,224,232] that ${\kappa }_1$ of NTA is useful in the identification of the last stages preceding fracture.

Kourkoulis et al. [228], among others, showed that, in notched marble specimens, criticality before fracture was first observed via an AE sensor that lay close to the crown of the notch where the fracture had started.

Niccolini et al. [230] in particular studied, through AE, an in-service double-girder bridge showing an asymmetric damage pattern. The NTA was able to identify criticality in the damage zone, while it revealed that the other monitored portions of the girder crane remained non-critical.

Also, recently, a study by Friedrich et al. [231] focused on using time-series methods to analyze AE data from two experiments with two totally different materials and structures: (i) a fiber-reinforced polymer plate undergoing a three-point bending test and (ii) a spaghetti bridge model. The NTA showed that the convergence of the order parameters and the entropies could identify the approach of a structure to a critical stage. This was possible either using the AE energy or by counting the ruptured samples. In both examples, it was possible to clearly perceive the correlation between the critical interval that was determined via the NTA from which the imminent increments in the acoustic emission activity followed.

Very recently, Triantis et al. [233] studied the evolution of both AE and electric activity for three kinds of concrete specimens under three-point bending until fracture. Specifically, centrally notched, beam-shaped, unreinforced concrete specimens, specimens reinforced with plastic, and specimens reinforced with metallic fibers were subjected to the same experimental protocol and loaded until fracture. The comparative analysis of AE and electric activities revealed that both sensing techniques, when analyzed in natural time, provided clear, in-time signals, indicating the entrance of the specimens into criticality. A further study of these signals, which are related to the damage mechanisms involved, may be very useful in engineering projects, such as structural health monitoring, EQ engineering, etc. [233].

3.8. Economics

NTA has the potential to be applied to problems related to economics. Mintzelas and Kiriakopoulos [234] introduced NTA to markets for price prediction and algorithmic trading, finding interesting results.

As an additional example, we present here an NTA of the Standard and Poor’s 500, or simply the S&P 500. As is well known, this is a stock market index that tracks the stock performance of the 500 largest companies included in the United States stock exchanges. Figure 7 shows the daily opening (open) values of the S&P 500, together with the daily volume according to Yahoo Finance [235] versus conventional time during the period of 1950 to 2012.

When setting as $Q_k$ the daily opening value of the S&P 500 and starting the NTA either on 1 January 1950, 1 January 1960, 1 January 1970, 1 January 1980, 1 January 1990, or 1 January 1995, we obtain the ${\kappa }_1$ values shown in Figure 7. Interestingly, we observe that, when the NTA was started on 1 January 1950, the criticality condition of Equation (3) was satisfied on 1 September 1987 (almost fifty days before the “Black Monday” crash on 19 October 1987; see the Reuters article written by Caroline Valetkevitch [236]). When starting the NTA of the S&P 500 on 1 January 1960, 1 January 1970, and 1 January 1980, Equation (3) was satisfied on 8 April 1998, 2 June 1997, and 19 August 1998, respectively. These dates mark the first significant increase in the S&P 500, and they are almost three years before the 24 March 2000 S&P 500 maximum during the dot-com bubble. A further inspection of Figure 7 reveals that, when the NTA started on 1 January 1990 and 1 January 1995, Equation (3) became valid on 19 March 2002 and 6 December 2001, respectively, which are several months before the S&P 500 minimum on 10 October 2002. Finally, when the S&P 500 was analyzed since 1 January 1995, we found that, again, ${\kappa }_1=0.070$ on 21 July 2005, i.e., almost two years before the 2007 S&P 500 maximum. We observed that this preliminary NTA of the S&P 500 revealed that Equation (3) identifies critical stages since it is precursory either to extreme events, like Black Monday, or to significant inflection (tipping) points in S&P 500 dynamics.

Certainly, the potential of NTA applications in economics has not been fully exploited, and many more such applications are expected in the future.

4. Summary and Conclusions

Natural time is based on the premise that the most important properties of a complex system can be revealed when we focus solely on the emitted energy of the events taking place in the complex system and on their order of occurrence. Using the normalized energy of these events and elements of probability theory, we were able to identify two key quantities: (i) the variance in natural time, ${\kappa }_1$, and (ii) the entropy, S, in natural time; see Section 2. Natural time analysis, introduced in 2001 [1,2,3] and is based on this novel concept of time called natural time, has found a multitude of applications in diverse, complex systems that range from geophysics (Section 3.2) to heart dynamics (Section 3.6) and from atmospheric sciences (Section 3.5) to economics (Section 3.8).

It enables (Section 3.3) a reliable estimation of the occurrence time (see, e.g., Figure 4) and epicentral area (cf. the EQ magnitude is estimated from the SES amplitude; Section 3.1) of impending major EQs, and it is currently considered the basis for the most recent method in seismic risk estimation by Turcotte, Rundle, and coworkers, termed earthquake nowcasting. A real-world example can be mentioned, the successful prediction of three EQs in Greece [6,7] based on NTA and the SES activities recorded via the VAN telemetric network shown in Figure 3. In the same sense, predicting the risk of natural hazards, the identification beforehand that the 2015–2016 and the 2023–2024 El Niño events were not to become some of the strongest was successful, as discussed in Section 3.5.2. Another real-world example is the medical studies presented in Section 3.6 concerning the discrimination of healthy individuals from those at high risk of sudden cardiac death by means of ECGs (Figure 5) and/or PPG (Figure 6). The application of NTA to engineering (Section 3.7) also includes real-life applications to various building materials, double-girder bridges, and structural health monitoring. Finally, applications to economics, like that in Figure 7, have the potential to become important real-world implementations of NTA.