Order Parameter and Entropy of Seismicity in Natural Time before Major Earthquakes: Recent Results

Abstract

A lot of work in geosciences has been completed during the last decade on the analysis in the new concept of time, termed natural time, introduced in 2001. The main advances are presented, including, among others, the following: First, the direct experimental verification of the interconnection between a Seismic Electric Signals (SES) activity and seismicity, i.e., the order parameter fluctuations of seismicity exhibit a clearly detectable minimum when an SES activity starts. These two phenomena are also linked closely in space. Second, the identification of the epicentral area and the occurrence time of an impending major earthquake (EQ) by means of the order parameter of seismicity and the entropy change of seismicity under time reversal as well as the extrema of their fluctuations. An indicative example is the M9 Tohoku EQ in Japan on 11 March 2011. Third, to answer the crucial question—when a magnitude 7 class EQ occurs—whether it is a foreshock or a mainshock. This can be answered by means of the key quantities already mentioned, i.e., the order parameter of seismicity and the entropy change of seismicity under time reversal along with their fluctuations. The explanation of the experimental findings identified before major EQs is given in a unified way on the basis of a physical model already proposed in the 1980s.

1. Introduction

There is a widespread belief that it is not space but time that in the end poses the greatest challenge in science. A monograph published in 2011 [1] reviewed the new view on time, termed natural time $\chi$ (from the Greek word “$\chi \rho$ó$\nu o\varsigma$” which means “time”), which was introduced by the authors in 2001 [2,3,4]. In this new view, time is not continuous, thus being in sharp contrast with the hitherto used conventional time t which is modeled as the one-dimensional continuum $\mathcal{R}$ of the real numbers.

The results deduced until 2011 on the basis of this new domain revealed that novel dynamical features hidden behind time series in complex systems can emerge upon analyzing them in natural time, but cannot do so when the analysis is carried out within the frame of conventional time. It was also shown that the analysis in natural time enables the study of the dynamical evolution of a complex system and identifies when the system enters a critical stage. This reflects that natural time plays a key role in predicting impending catastrophic events in general. In several chapters of this monograph [1] examples of data analysis in natural time were presented (that have been published mainly in Physical Review and Physical Review Letters) in diverse fields including Earth Sciences (Geophysics, Seismology), Physics (Condensed Matter, Statistical Physics, Physics of Complex Systems), Biology and Cardiology.

In [1] (Chapter 2), it was demonstrated that natural time is optimal for enhancing the signals in time–frequency space when employing the Wigner function and measuring its localization property. In other words, natural time analysis conforms to the desire to reduce uncertainty and extract signal information as much as possible. The normalized power spectrum $\Pi \left(\omega \right)$ introduced in natural time whose Taylor expansion leads, at low natural (cyclic) frequencies $\omega$ ($\omega \to 0$), to the expression $\Pi \left(\omega \right)\approx 1-{\kappa }_1{\omega }^2$ where the coefficient ${\kappa }_1$ is just the variance of natural time, i.e., ${\kappa }_1=\langle {\chi }^2\rangle -{\langle \chi \rangle }^2$. This quantity is useful in identifying the approach to a critical point. In the natural time analysis of the seismicity, a careful inspection reveals (see [5]; see also Chapter 6 of [1]) that the quantity ${\kappa }_1$ may be considered as an order parameter of seismicity. This allows the determination of the constant b in the Gutenberg–Richter law for earthquakes (EQs), $N(\ge M)={10}^{a-bM}$, by applying the Maximum Entropy Principle [6,7]. It leads to $b\approx 1$, which agrees with real seismic data.

In [1] (Chapter 3, see also [8,9]), the entropy S in natural time was defined and additionally it was shown that the entropy $S_{-}$ deduced from the natural time analysis of the time series obtained upon time reversal (i.e., reversing the time arrow) is in general different from S, thus, the entropy in natural time does satisfy the condition to be “causal” [9,10]. In addition, the physical meaning of the change $\Delta S\equiv S-S_{-}$ of the entropy in natural time under time reversal was discussed, which is of profound importance when studying [11] the evolution of complex systems. Furthermore, complexity measures were introduced that quantify the fluctuations of the entropy S and of the quantity $\Delta S$ upon changing the length scale as well as the extent to which they are affected when shuffling randomly the consecutive events.

Since the publication of the aforementioned monograph [1], however, i.e., in the last decade, a lot of work on the natural time analysis of seismicity has been published which has shed light on the crucial importance of the order parameter of seismicity (${\kappa }_1$) and on the entropy change $\Delta S$ under time reversal, along with their fluctuations, in the identification of important changes before major EQs. It is the main aim of this work to review these precursory phenomena which may be of broad interest in the geosciences.

One of these phenomena is the finding that the minimum of the fluctuations of the order parameter of seismicity appears almost simultaneously with the initiation of Seismic Electric Signals (SES) activity, which is a series of transient low-frequency Earth electric field changes (whose main properties were recapitulated in the monograph [1] in 2011) preceding major EQs (see also Section 4). This simultaneous appearance, which was experimentally confirmed and reported in 2013, is of profound importance, because the former quantity (fluctuations of the order parameter of seismicity) comes from seismic data alone, while the second one (SES activity) solely from electric field measurements. The correlation of these two types of data has been suggested on the basis of a physical model from the 1980s in which it is underlined that SES is a critical phenomenon (cf. major EQs are also considered as critical phenomena, where a mainshock is the new phase).

In view of the above, the present feature article is structured as follows: In Section 2, we summarize the latter physical model on the SES generation, based on the thermodynamics of point defects, a strict foundation of which was completed in the beginning of the 1980s [12]. In the following two main parts, the precursory changes of the two quantities ${\kappa }_1$ and $\Delta S$ of seismicity—along with their fluctuations—are described in the Section 3, Section 4, Section 5, Section 6, Section 7 and Section 8 and Section 9, Section 10 and Section 11, respectively. Finally, a discussion is presented in Section 12 and our conclusions are given in the last section.

Note that the new method for estimating the seismic risk, termed nowcasting, by Turcotte, Rundle and coworkers [13,14,15,16,17,18,19], is also based on the concept of natural time.

2. The Physical Model for the Seismic Electric Signals Generation

The physical model for SES generation [20] (see Figure 1) also foresees (see [21]) that some additional transient multidisciplinary phenomena should be simultaneously generated and completed well before the EQ rupture, which was actually observed (see Figure 3 of Ref. [21]) before the M9 Tohoku EQ of 11 March 2011 which will be discussed later. This essentially differs from other precursory mechanisms proposed in [22] that usually exhibit anomalous behavior becoming more intense upon approaching the EQ failure as seen in Figure 1b of [21] which is reproduced here in Figure 2b. This physical model is termed ”pressure stimulated polarization currents (PSPC) model” ( [12]; summarized also in [23,24,25,26]), and suggests the following (see Figure 1): In the Earth, electric dipoles always exist [12] due to lattice imperfections (point and linear defects) in the ionic constituents of rocks. In the future focal region of an EQ, where the electric dipoles have initially random orientations (Figure 1c), the stress, $\sigma$, starts to gradually increase due to an excess stress disturbance (Figure 1a). Let us call hereafter this stage A. When this gradually increasing stress reaches a critical value (${\sigma }_{cr}$), the electric dipoles exhibit a cooperative orientation (Figure 1e) resulting in the emission of a transient electric signal which constitutes an SES (Figure 1b) with current density j. (Note that “cooperativity” is a hallmark of criticality [27]). We call this stage B. Uyeda et al. [28] pointed out that the PSPC model is unique among other models in that SES would be generated spontaneously during the gradual increase of stress without requiring any sudden change of stress such as microfracturing [29] (or faulting). It should be emphasized that in this model the concept that SES is a critical phenomenon emerges in view of the “cooperativity” mentioned above.

Several SES recorded within a short time are termed SES activity [30]. Intense SES activities precede major EQs along with clear magnetic field variations [31] mainly on the z-component [32,33].

3. Minimum of the Order Parameter Fluctuations of Seismicity before Major Earthquakes

For a time series comprising of N EQs, we define the natural time for the occurrence of the k-th EQ ${\chi }_k=k/N$ and $p_k=Q_k/{\sum }_{n=1}^N{Q}_n$ is the normalized energy of the k-th EQ of energy $Q_k$ calculated from seismic catalogs. The quantity ${\kappa }_1$ serves, as mentioned in the Introduction, as an order parameter [1,5,34,35]. We recall that

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

(1)

To obtain the fluctuation $\beta$ of ${\kappa }_1$, we need many values of ${\kappa }_1$ for each target EQ. For this purpose, we first take an excerpt comprising W successive EQs just before a target EQ from the seismic catalog. The number W is chosen to cover a period of a few months. Thus, W is the number of EQs that on average occur during the average lead time of an SES activity that ranges from a few weeks to $5\frac{1}{2}$ months (see Chapter 7 of Ref. [1]). For this excerpt, we form its sub-excerpts $S_j={\left\{Q_{j+k-1}\right\}}_{k=1,2\dots ,N}$ of consecutive N = 6 EQs (since at least six EQs are needed [5] for obtaining reliable ${\kappa }_1$) of energy $Q_{j+k-1}$ and natural time ${\chi }_k=k/N$ each. Further, $p_k=Q_{j+k-1}/{\sum }_{k=1}^N{Q}_{j+k-1}$, and by sliding $S_j$ over the excerpt of W EQs, $j=1,2,\dots ,W-N+1(=W-5)$, we calculate ${\kappa }_1$ using Equation (1) for each j. We repeat this calculation for $N=7,8,\dots ,W$, thus obtaining an ensemble of $\left[(W-4)(W-5)\right]/2(=1+2+\dots +W-5)$${\kappa }_1$ values. Then, we compute the average $\mu \left({\kappa }_1\right)$ and the standard deviation $\sigma \left({\kappa }_1\right)$ of thus obtained ensemble of $\left[(W-4)(W-5)\right]/2$${\kappa }_1$ values. The variability $\beta$ of ${\kappa }_1$ for this excerpt W is defined to be:

$$\beta \equiv \sigma \left({\kappa }_1\right)/\mu \left({\kappa }_1\right)$$

(2)

and is assigned to the $(W+1)$-th EQ in the catalog, i.e., the target EQ.

The time evolution of the $\beta$ value can be pursued by sliding the excerpt through the EQ catalog. Namely, through the same process as above, $\beta$ values assigned to $(W+2)\mathrm{th}$, $(W+3)\mathrm{th},\dots$ EQs in the catalog can be obtained. Note that $\beta$ may bear a subscript ${\beta }_W$ to clarify the window length W used for its calculation.

Since 2013 (see [36,37]) the calculations have been made by taking all natural time windows comprising 6 to W events following the procedure described above. On the other hand, before 2013 (as for example is the case in [38]), the calculations were made as follows: First, for an event considering natural time windows from 6 to 40 events (i.e., selecting the precise value 40 as an upper limit). Second, this process was performed for all the events sliding through the whole data set (e.g., an EQ catalog).

In [36], the Japan seismic catalog was analyzed in natural time from 1 January 1984 until the occurrence time of the M9 Tohoku EQ on 11 March 2011 considering a natural time window of fixed length comprising the number of events that would occur in a few months. It was found that the fluctuations of the order parameter of seismicity, computed by the procedure described above, exhibit minima a few months before all shallow EQs of magnitude 7.6 or greater during this 27-year period in the Japanese area (Figure 3). Among the minima, the one before the M9 Tohoku EQ was the deepest. In this analysis, the Japan Meteorological Agency (JMA) seismic catalog was used by considering all EQs in the period from 1984 to the occurrence time of the M9 Tohoku EQ, within the area ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{148}$ (Figure 3). The energy of EQs was deduced from $M_{JMA}$ after converting [39] to the moment magnitude $M_w$ [40]. Considering a threshold $M_{JMA}$ = 3.5 for the sake of data completeness, we are left with 47,204 EQs for the period under study, which is around 326 months. Thus, on average, $\sim {10}^2$ EQs occur per month. We considered the values W = 200, 300 and 400, which occur during a period of a few months, which is comparable with the average lead time of SES activities recorded in Japan [41] and Greece [26,30,42].

We first focus on the minimum of the variability $\beta$ observed before the M9 Tohoku EQ. Figure 4 depicts about 47,200 $\beta$-values calculated for W = 300 versus the target EQ number from 1984 to the day of the Tohoku EQ, 11 March 2011. EQs with $M_{JMA}\ge 6.9$ ($M_{JMA}$ in the right scale) are indicated by blue asterisks. The $\beta$-values show the deepest minimum (min) just before the occurrence of the Tohoku EQ (rightmost of Figure 4A). Figure 4B is an expanded version, in the conventional time, of the last 10-month period shown in yellow in Figure 4A. The red, blue and green curves depict what happened to $\beta$ for W = 200, 300 and 400. Hereafter, we use the symbols ${\beta }_W$ and ${\beta }_{W,min}$ as needed. After around 1 September 2010, a decrease of ${\beta }_W$ became evident and ${\beta }_W$ fell to a minimum in early January 2011, about 2 months before the Tohoku EQ. The results of this minimum of $\beta$ are summarized as follows (Figure 4A,B): First, such a minimum of $\beta$ was not observed again during the 27-year period studied. Second, the ratio ${\beta }_{300,min}/{\beta }_{200,min}$ is near unity. Third, the dates of ${\beta }_{W,min}$ for W = 200, 300 and 400 are 5 January, 5 January and 10 January 2011, respectively, i.e., almost the same. Fourth, this minimum is less clear for greater W corresponding to time intervals longer than a few months. It is almost invisible for W = 2000 and 3000 (Figure 4C). The same applies to all other ${\beta }_{W,min}$, as seen in Figure 4C. In what follows, we restrict ourselves to W = 200 and W = 300.

We now proceed to the $\beta$ minima before other major EQs in Japan during the 27-year period investigated. Six shallow EQs with $M_{JMA}$ 7.6 or larger (Figure 3) occurred:

EQa 1993-07-12: 1993 Southwest-Off Hokkaido EQ ($M_{JMA}$ = 7.8);

EQb 1994-10-04: 1994 East-Off Hokkaido EQ ($M_{JMA}$ = 8.2);

EQc 1994-12-28: 1994 Far-Off Sanriku EQ ($M_{JMA}$ = 7.6);

EQd 2003-09-26: 2003 Off Tokachi EQ ($M_{JMA}$ = 8.0);

EQe 2010-12-22: 2010 Near Chichi-jima EQ ($M_{JMA}$ = 7.8);

EQf 2011-03-11: 2011 Tohoku EQ ($M_w$ = 9.0).

It was examined whether $\beta$ minimum precede these EQs as well. Figure 5A–C are the expanded versions in the conventional time of Figure 4A in three 10-year periods. The EQs are marked by a–f. In view of the fact that these figures are still small, the time axis was expanded for each EQ, as shown in Figure 6A–E, just as in Figure 4B for the Tohoku EQ. We now see $\beta$ minima within 1–3 months before every one of the six mainshocks as shown in Table 1 of [36]. In this table, the ${\beta }_{300,min}/{\beta }_{200,min}$ ratio and $\Delta t_{200}$ of $\beta$ minima before each of these EQs are very similar to those before the Tohoku EQ, e.g., the ratios lie in the range 0.95–1.08. It was then concluded that these minima are precursory to the time-correlated EQs. Beyond these $\beta$ minima before the six $M_{JMA}\ge 7.6$ EQs, there were many more minima during the 27-year period, as seen in Figure 5. We examined whether they were also followed by EQs. Along these lines, the minima were chosen deeper than the shallowest one of the six ${\beta }_{200,min}$ in Table 1 of [36], which happened before EQb, the 1994 East-Off Hokkaido EQ ($M_{JMA}$ = 8.2), yielding 0.295 for ${\beta }_{200,min}$. We found that, out of the thus chosen 31 minima, 9 (labeled 1–9 in Figure 5) also displayed ${\beta }_{300,min}/{\beta }_{200,min}$ ratio (Figure 5, Table 2 of [36]) in the range 0.95–1.08, thus almost equal to those observed before the six $M_{JMA}\ge 7.6$ EQs (Table 1 of Ref. [36]). These nine minima have been followed by $M_{JMA}\ge 6.4$ EQs within 3 months (Figure 5 and $\Delta t_{200}$ in Table 2 of Ref. [36]). Such correspondences are naturally less certain, because of a larger number of EQs. In particular, there were 139 $M_{JMA}\ge 6.4$ EQs during the 27-year period studied.

In summary, we conclude that analysing in natural time the seismicity of Japan from 1 January 1984 to the occurrence time of M9 Tohoku EQ on 11 March 2011, using a natural time window comprising the number of events that would occur in a few months (this is called hereafter crucial scale for the reasons discussed in Ref. [43]), the following results were obtained.

Almost 2 months before the M9 Tohoku EQ, a minimum in the variability $\beta$ of the order parameter of seismicity ${\kappa }_1$ is observed which is the deepest in the whole study period. Distinct minima of $\beta$, but of shallower depth, were found also one month to a few months before the occurrence of all other Japanese major EQs ($M_{JMA}\ge 7.6$, depth <400 km) during 1984–2011. The $\beta$ minima, which seem to be precursory to sizeable EQs show the ${\beta }_{300,min}/{\beta }_{200,min}$ ratio in the range of 0.95–1.08.

The approximate coincidence of the minima of $\beta$ with the initiation of the SES activities identified by Varotsos et al. [38]—which will be discussed in Section 4—strengthens the physics of both phenomena.

3.1. Statistical Significance of the Minima of the Order Parameter Fluctuations

In summary, the results obtained by Sarlis et al. [36] are as follows: Upon analyzing the JMA seismic catalog [36] in natural time and considering all EQs of magnitude $M_{JMA}\ge 3.5$ from 1 January 1984 to 11 March 2011 (the time of the M9 Tohoku EQ) within the area ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{148}$, fifteen distinct minima—observed simultaneously (see, e.g., Appendix A in [44]) at ${\beta }_{200}$ and ${\beta }_{300}$ with a ratio ${\beta }_{300,min}/{\beta }_{200,min}$ in the range 0.95 to 1.08 and ${\beta }_{200,min}\le 0.295$—of the fluctuations of the order parameter of seismicity were found 1 to around 3 months before large EQs. All shallow EQs of magnitude 7.6 or larger during this 27-year period were preceded by six (out of 15) of these minima ${\beta }_{W,min}$. Remarkably, among the minima, the deepest minimum was observed around 5 January 2011, i.e., around two months before the M9 Tohoku EQ.

In 2016, Sarlis et al. [45] showed, using Monte Carlo calculations, that the probability to achieve the above results by chance is of the order of ${10}^{-5}$. This conclusion was also strengthened in 2020 by Christopoulos et al. [46], using the event coincidence analysis (ECA) [47]. In addition, in 2020 Sarlis et al. [48] found the same value (${10}^{-5}$) by means of the receiver operating characteristics (ROC) technique by focusing on the area under the ROC curve (AUC). They explained that, if this area, which is currently considered an effective way to summarize the overall diagnostic accuracy of a test, has the value 1, it corresponds to a perfectly accurate test. It was found [48] that the AUC in these results is around 0.95, which is evaluated as outstanding. This result indicates that the study of the minima ${\beta }_{W,min}$ of seismicity for an area covered by a well-developed seismological network such as Japan may be used for operational EQ forecasting.

3.2. Minimum of the Order Parameter of Seismicity Preceding the Ultra-Deep 2015 EQ ($M_{w}7.9$) in Japan

On 30 May 2015, a powerful EQ ($M_w$7.9, JMA reported magnitude $M8.1$) struck west of Japan’s remote Ogasawara (Bonin) island chain, which lies more than 800 km south of Tokyo. It occurred at 680 km depth in an area without any known historical seismicity and caused significant shaking over a broad area of Japan at epicentral distances in the range 1000–2000 km. Specifically, the maximum shaking intensity was reached at both Hahajima Island (above the hypocenter) and Kanagawa (near Tokyo) located over 800 km from the hypocenter. It was the first EQ felt in every Japanese prefecture since intensity observations began in 1884. By applying natural time analysis, Varotsos et al. [49] found that almost three and a half months before this powerful EQ, which is the strongest one after the M9 Tohoku EQ on 11 March 2011, the fluctuations of the order parameter of seismicity were minimized around 17 February 2015. We have to mention that no SES has been reported for this undersea ultra-deep EQ.

4. The Experimental Fact Showing the Physical Interconnection of SES Activity with Seismicity

In 2013, in [38] it was found that the fluctuations $\beta$ of the order parameter of seismicity exhibited a clearly detectable minimum approximately at the time of the initiation of a pronounced SES activity recorded by Uyeda and coworkers [41,50], around 2 months before the volcanic-seismic swarm activity in 2000 in the Izu Island region, Japan. Specifically, this SES activity started [41,50] on 26 April 2000, and the swarm was then characterized by JMA as being the largest EQ swarm ever recorded [51].

To meet this goal, we analyzed in natural time the EQs reported during this period in the JMA seismic catalog. In particular, we considered all EQs within the area ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{146}$ (see the map in Figure 3). Setting a threshold $M_{JMA}>3.4$ to assure data completeness, there exist 52,718 EQs in the period from 1967 to the time of Tohoku EQ. This reflects that we have on average $\sim {10}^2$ EQs per month.

This is the first time that, well before the occurrence of major EQs, anomalous changes were found to appear approximately simultaneously in two independent datasets of different geophysical observables (geoelectrical measurements, seismicity). In addition, it was shown [38] that these two phenomena are also linked closely in space as follows: By using a narrow $5^{\circ }\times 5^{\circ }$ spatial window sliding through the wide (${21}^{\circ }\times {21}^{\circ }$) area ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{146}$, we investigated the variability of the order parameter of seismicity corresponding to EQ events that would occur in an almost two-month period. We found that the characteristics of the fluctuations of the order parameter of seismicity (i.e., their lowest minimum) are linked to the seismicity occurring in the areas ${\mathrm{N}}_{33}^{38}{\mathrm{E}}_{135}^{140}$ and ${\mathrm{N}}_{31}^{36}{\mathrm{E}}_{135}^{140}$ that include the Izu Island region which became active a few months later.

5. Estimation of the Epicenter by Means of the Spatiotemporal Variations of the Minimum of the Seismicity Order Parameter Fluctuations

In an earlier study, it was found (as described in Section 3) that the fluctuations $\beta$ of the seismicity order parameter calculated for the entire Japanese region were minimized a few months before the shallow major EQs (magnitude larger than 7.6), which took place in the region from 1 January 1984 to 11 March 2011. Here, by dividing the Japanese region ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{148}$ into small areas, the $\beta$ calculation on them was carried out by following [37]. It was found that some small areas show $\beta$ minimum almost simultaneously with the large area and such small areas clustered within a few hundred kilometers from the actual epicenter of the related mainshocks. These results revealed [37] that this approach may estimate the epicentral location of forthcoming major EQs.

In short, the procedure is as follows (cf. for details [37]): The data source is the same JMA seismic catalog. We set circular areas with radius R = 250 km of which the center is sliding through the large area with steps of 0.1${}^{\circ }$ in longitude and latitude. We worked on the time changes of $\beta$ for these areas. From these small areas that showed “local” $\beta$ minimum, we selected the ones where the date of $\beta$ minimum coincided (i.e., ±2 days) with the one in the large area (whole Japanese region). We started our investigation at 5.5 months before each major EQ. The reason for this was that 5.5 months is the maximum lead time of SES activities observed to date. To assure that a “local” $\beta$ minimum is clearly recognizable, the criterion that it should differ more than 10% from the $\beta$ value of the events that occurred within 10 days before and after was imposed. When several “local” $\beta$ minima appeared simultaneously (±2 days) with the $\beta$ minima in the large area in many small areas, we investigated the spatial distribution of their centers. We found that, for all of the six shallow EQs of magnitude larger than 7.6 that occurred in Japan from 1 January 1984 to 11 March 2011, a large number of small areas exhibited $\beta$ minimum almost simultaneously with the large area. Such small areas are accumulated in a region that lies within a few hundred kilometers of the actual epicenter.

6. Identification of the Occurrence Time of a Mainshock by Means of the Minimum of the Order Parameter Fluctuations

In the time series analysis using natural time, the behavior of the normalized power spectrum [2,3,4]

$$\Pi \left(\omega \right)\equiv {|\Phi \left(\omega \right)|}^2$$

(3)

defined by

$$\Phi \left(\omega \right)=\sum _{k=1}^N{p}_{k}exp\left(i\omega {\chi }_k\right)$$

(4)

where $\omega$ stands for the angular frequency at $\omega$ 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 $\Phi \left(\omega \right)$ at $\omega \to 0$ (see p. 499 of [52]). For this purpose, a quantity ${\kappa }_1$ is defined from the Taylor expansion $\Pi \left(\omega \right)=1-{\kappa }_1{\omega }^2+{\kappa }_2{\omega }^4+\dots$. The relation for the critical state that has been shown for SES activities [1,2]:

$$\Pi \left(\omega \right)=\frac{18}{5{\omega }^2}-\frac{6cos\omega }{5{\omega }^2}-\frac{12sin\omega }{5{\omega }^3}.$$

(5)

for $\omega \to 0$, simplifies to

$$\Pi \left(\omega \right)\approx 1-0.07{\omega }^2$$

(6)

which shows that the second-order Taylor expansion coefficient of $\Pi \left(\omega \right)$, i.e., ${\kappa }_1$, is equal to 0.070. Equation (6) has been shown to hold also for EQ models like the Burridge–Knopoff EQ model [53,54] as well as in the Olami Feder Christensen (OFC) EQ model [1,55] when they approach the critical point.

In general, two procedures have been suggested in order to determine the occurrence time of an impending mainshock: A procedure, which solely for the sake of convenience has been called preliminary procedure, was followed for several cases (e.g., see References [2,5,7,33,56,57]) by starting the natural time analysis of seismicity in the candidate area A immediately after the SES initiation, because the latter signals that the system enters the critical stage and the time variation of parameters should be traced only on the single candidate area A. To identify, however, such a single area, one should know the properties of the preceding SES activity in general, or, in particular, the selectivity map [30,58] of the measuring station that recorded the SES activity. An alternative procedure, called “updated” procedure, was later developed [42], that considers the natural time analysis of the seismicity in all the possible subareas of A, instead of a single area identified by the SES selectivity map. This has been inspired as follows [42]: When area A reaches criticality, one expects in general that all its subareas have also reached criticality simultaneously. At that time, therefore, the evolution of seismicity in each of these subareas is expected to result in a ${\kappa }_1$ value close to 0.070. Assuming equipartition of probability among the subareas, the distribution Prob(${\kappa }_1$) of the ${\kappa }_1$ values of all subareas should peak at around 0.070, exhibiting also magnitude threshold invariance (for more details on how the subareas are selected, see page 300 of [1]). This usually enables the estimation of the occurrence time of major EQs with a time window of less than a week or so.

When geoelectrical measurements are not available, we cannot identify directly the initiation time of an SES activity. Since, however, this initiation occurs almost at the time of the minimum of the fluctuations of the order parameter of seismicity in a large area, we take advantage of the occurrence time of this minimum. In addition, a spatiotemporal study of this minimum described in the preceding section enables an estimation of the epicentral area of the impending major EQ [37]. In particular, we will present below the case of the identification of the occurrence time of the impending M9 Tohoku mega-earthquake that occurred on 11 March 2011 in Japan.

We start the natural time analysis of seismicity at the date of 5 January 2011, at which the fluctuations of the order parameter of seismicity exhibited [36] the deepest minimum ever observed during the period 1984–2011 in Japan. This date should almost coincide with the initiation of an SES activity [38], as explained in Section 4. Actually, anomalous magnetic field variations appeared during the period of 4 to 14 January 2011 on the z component at two measuring sites (Esashi and Mizusawa) lying at epicentral distances of around 130 km [59,60,61,62] (which reflect also the existence of an SES activity, see, e.g., [32,33]).

Since the SES selectivity maps of the aforementioned measuring stations are lacking, we proceed to the estimation of the epicentral location of the impending mainshock without making use of SES data, and follow the procedure described in Section 5 [37]: Dividing the entire Japanese region ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{148}$ (large area) into small areas, a calculation of the fluctuations of ${\kappa }_1$ of seismicity is carried out on them. Some small areas show a minimum of the fluctuations almost simultaneously with the minimum in the entire Japanese region (on 5 January 2011) and such small areas cluster within a few hundred km from the actual epicenter. Such a calculation results in an estimate of the candidate epicentral area shown by the blue-green area in Figure 7.

Starting from 5 January 2011, we now compute the ${\kappa }_1$ values of seismicity in the blue-green area and the results are depicted in Figure 8 for $M_{thres}$ = 4.2 to 5.0. Considering that at least six EQs are needed [5] for obtaining a reliable ${\kappa }_1$ value, which happens in the blue-green area on 16 February 2011 for $M_{thres}$ =4.2, we show the computed ${\kappa }_1$ values during the last four weeks before the M9 Tohoku EQ in Figure 8A. This figure reveals that the condition ${\kappa }_1=0.070$ is not satisfied for all magnitude thresholds at least until the M7.3 EQ on 9 March 2011. Since the results afterwards cannot be seen clearly in Figure 8A, we plot in expanded time scale in Figure 8B (see also Figure 9) the ${\kappa }_1$ values of seismicity from 00:00 LT on 9 March 2011 until the Tohoku EQ occurrence. This figure clearly shows that the condition ${\kappa }_1=0.070$ (which signals that the mainshock is going to occur within the next few days or so) is fulfilled for $M_{thres}$ = 4.2 to 5.0 in the morning of 10 March 2011 upon the occurrence of the EQs from 08:36 to 13:14 LT, i.e., almost one day before the Tohoku EQ (see the gray shaded area in Figure 8B), thus, the M7.3 EQ on 9 March 2011 was a foreshock.

We now return to the same example, but using the preliminary procedure according to which the criteria to assure a true coincidence of the EQ time series with that of critical state are [1]:

• The “average” distance $\langle D\rangle$ between the curves of $\Pi \left(\omega \right)$ of the evolving seismicity and Equation (5) for $\omega \le \pi$ should be $\langle D\rangle <{10}^{-2}$, for details, see [1,2];
• The final approach of the evolving $\Pi \left(\omega \right)$ to that of Equation (5) must be from below. This reflects that ${\kappa }_1$ gradually decreases with time before strong EQs finally approaching from above that of the critical state, i.e., ${\kappa }_1=0.070$, see Figure 7.1 on page 293 of [1];
• At the coincidence, both entropies S and $S_{-}$ must be smaller than $S_u$;
• Since this process (critical dynamics) is considered to be self-similar, the occurrence time of the true coincidence should not markedly vary upon changing the magnitude threshold $M_{thres}$.

In light of new advances (see Section 9), it was shown (see, e.g., [63]) that the second criterion for the true coincidence starts to be fulfilled on the date at which $\Delta S$ (cf. of the seismicity in the candidate epicentral area) minimizes. In other words, the quantity ${\kappa }_1$ after the $\Delta S$ minimization starts monotonically decreasing from values ${\kappa }_1>0.070$ and approaches finally from above the value ${\kappa }_1=0.070$ very close to the impending strong EQ occurrence date. Specifically, in Figure 10, we show the evolution of the $\Delta S$ values of seismicity (left scale) in the blue-green area of Figure 7 almost a month (a) and two-and-a-half days (b) before the M9 Tohoku EQ on 11 March 2011, i.e., the same periods as those depicted in Figure 8 for the evolution of the ${\kappa }_1$ value. An inspection of this figure shows that the fulfillment of the second criterion for the true coincidence of the EQ time series with that of critical state begins around 06:24 LT on 10 March 2011, thus agreeing with the time of the minimum of $\Delta S$ shown in Figure 10. In other words, it is of major importance that after the appearance of the $\Delta S$ minimum of seismicity, the order parameter begins to gradually diminish, thus approaching the critical value ${\kappa }_1=0.070$, which signals that a major EQ is going to shortly occur in the blue-green area of Figure 7.

Figure 7. Estimating the M9 Tohoku EQ epicenter: Color contours of the number of the small areas exhibiting the minimum of the fluctuations of ${\kappa }_1$ almost simultaneous with that in the entire Japanese region. The red diamond depicts the M9 Tohoku EQ epicenter and the red triangle the location of the Esashi measuring magnetic station (the other measuring station, Mizusawa, lies only a few tens of km away from Esashi). The epicenters of the EQs that occurred within the candidate epicentral area (blue-green) after 5 January 2011 and before the M9 Tohoku EQ for various magnitude thresholds $M_{thres}$ are also shown. Taken from Ref. [64].

Figure 7. Estimating the M9 Tohoku EQ epicenter: Color contours of the number of the small areas exhibiting the minimum of the fluctuations of ${\kappa }_1$ almost simultaneous with that in the entire Japanese region. The red diamond depicts the M9 Tohoku EQ epicenter and the red triangle the location of the Esashi measuring magnetic station (the other measuring station, Mizusawa, lies only a few tens of km away from Esashi). The epicenters of the EQs that occurred within the candidate epicentral area (blue-green) after 5 January 2011 and before the M9 Tohoku EQ for various magnitude thresholds $M_{thres}$ are also shown. Taken from Ref. [64].

Figure 8. Determination of the occurrence time of M9 Tohoku EQ. The evolution of the ${\kappa }_1$ values of seismicity (left scale) in the blue-green area of Figure 7 along with the EQs (right scale) shown with thick black vertical bars are plotted against the conventional time: (A) From 12 February 2011 until the M9 Tohoku EQ occurrence on 11 March 2011. (B) From 00:00 LT on 9 March 2011 until the Tohoku EQ occurrence. The shaded area shows the period in the morning of 10 March 2011 during which the condition ${\kappa }_1=0.070$ is fulfilled (see the text). Taken from Ref. [64].

Figure 8. Determination of the occurrence time of M9 Tohoku EQ. The evolution of the ${\kappa }_1$ values of seismicity (left scale) in the blue-green area of Figure 7 along with the EQs (right scale) shown with thick black vertical bars are plotted against the conventional time: (A) From 12 February 2011 until the M9 Tohoku EQ occurrence on 11 March 2011. (B) From 00:00 LT on 9 March 2011 until the Tohoku EQ occurrence. The shaded area shows the period in the morning of 10 March 2011 during which the condition ${\kappa }_1=0.070$ is fulfilled (see the text). Taken from Ref. [64].

Figure 9. The evolution of the ${\kappa }_1$ values of seismicity in a similar fashion as in Figure 8B, but plotted versus natural time $\chi$ ($\chi =1$ corresponds to the occurrence of the M9 Tohoku EQ, while $\chi \approx$ 0.2 to the M7.3 foreshock on 9 March 2011). Taken from Ref. [64].

Figure 9. The evolution of the ${\kappa }_1$ values of seismicity in a similar fashion as in Figure 8B, but plotted versus natural time $\chi$ ($\chi =1$ corresponds to the occurrence of the M9 Tohoku EQ, while $\chi \approx$ 0.2 to the M7.3 foreshock on 9 March 2011). Taken from Ref. [64].

Figure 10. The evolution of the $\Delta S$ values of seismicity (left scale) in the blue-green area of Figure 7 almost a month (a) and two-and-a-half days (b) before the M9 Tohoku EQ on 11 March 2011, i.e., the same periods as those depicted in Figure 8 for the evolution of the ${\kappa }_1$ value. Concerning the importance of the appearance of the minimum of $\Delta S$, which is evident in (b), see the fulfillment of the criteria for the approach to criticality discussed in Section 6.

Figure 10. The evolution of the $\Delta S$ values of seismicity (left scale) in the blue-green area of Figure 7 almost a month (a) and two-and-a-half days (b) before the M9 Tohoku EQ on 11 March 2011, i.e., the same periods as those depicted in Figure 8 for the evolution of the ${\kappa }_1$ value. Concerning the importance of the appearance of the minimum of $\Delta S$, which is evident in (b), see the fulfillment of the criteria for the approach to criticality discussed in Section 6.

In summary, the natural time analysis of seismicity results in the conclusion that the ${\kappa }_1$ values converge to 0.070 almost one day before the Tohoku EQ occurrence, i.e., the system approaches the critical point almost one day before the mainshock.

Alternatively, one may say that natural time analysis allowed for the recognition that the approach to the critical point happened almost a day after the occurrence of the M7.3 EQ of 9 March 2011, thus recognizing that this M7.3 EQ was a foreshock.

The identification of the occurrence time of a mainshock by using the minimum of the fluctuations of the order parameter of seismicity—instead of the initiation of an SES activity—was also achieved for the M8.2 EQ on 7 September 2017 [65], which was the strongest EQ in Mexico for more than a century, the deadly M7.1 EQ on 19 September 2017 in Mexico [63], and for the Ridgecrest M7.1 EQ on 6 July 2019 in California [66]. At this point, we have to mention that in the latter case [66], a natural time method introduced [67,68] for the estimation of the occurrence time of strong aftershocks could also provide an estimate of the occurrence time of the Ridgecrest M7.1 EQ.

7. Temporal Correlations in the Magnitude Time Series and the Minimum of the Order Parameter Fluctuations of Seismicity

Focusing on the minima of the order parameter fluctuations of seismicity preceding all EQs of magnitude 8 (and 9) class that occurred in a Japanese area from 1 January 1984 to the M9 Tohoku EQ occurrence on 11 March 2011 and applying Detrended Fluctuation Analysis (DFA) to the EQ magnitude time series, it was found [44] that each of these minima is preceded as well as followed by characteristic changes of temporal correlations between EQ magnitudes. The following three main features were identified: The minima are observed during periods when long-range correlations have been developed, but they are preceded by a stage in which an evident anti-correlated behavior appears. After the minima, the long-range correlations break down to an almost random behavior turning to anti-correlation. This is strikingly similar to the findings in other complex time series as, for example, the case of electrocardiograms. In this case, the long-range temporal correlations that characterize the healthy heart rate variability break down for individuals with a high risk of sudden cardiac death and often accompanied by the emergence of uncorrelated randomness [1,69].

The aforementioned minima ${\beta }_{W,min}$ of the order parameter fluctuations that precede $M\ge 7.8$ EQs can be distinguished from other minima that are either non-precursory or followed by smaller EQs on the basis of certain criteria developed in [44] (see its Appendix). They could be summarized as follows: In order to classify a ${\beta }_{W,min}$ value, it should be a minimum for at least W values before (bef) and W values after (aft). Further, to assure that ${\beta }_{200,min}$, ${\beta }_{250,min}$ and ${\beta }_{300,min}$ (where the subscript $W=200$, 250 and 300 events stands for the natural time window length used in the calculation) are precursory of the same mainshock and hence belong to the same critical process, almost all (in practice above 90% of) the events which led to ${\beta }_{200,min}$ should participate in the calculation of ${\beta }_{250,min}$ and ${\beta }_{300,min}$.

To distinguish the ${\beta }_{W,min}$ that are precursory to EQs of magnitude $M_{JMA}\ge 7.8$ from other minima (min) which are either non-precursory or may be followed by smaller magnitude EQs, we create separate studies for the larger and the smaller area, i.e., ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{148}$ and ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{146}$, respectively, see Figure 3, and the results obtained should be necessarily checked for their self-consistency. For example, a major EQ whose epicenter lies in both areas should be preceded by ${\beta }_{W,min}$ identified in the separate studies of these two areas approximately (in view of their difference in seismic rates) on the same date. In particular, we work in the following way (see the Appendix of [44]).

Let us assume that we start the investigation from the larger area where five EQs of magnitude 7.8 or larger occurred from 1 January 1984 until the M9 Tohoku EQ in 2011 (

Table 1. All $M_{JMA}\ge 7.8$ EQs since 1 January 1984 until the M9 Tohoku EQ on 11 March 2011 in the larger area ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{148}$. Reproduced under permission from Ref. [44]. ©2014. American Geophysical Union. All Rights Reserved.

Label	EQ Date	EQ Name	Lat. (${}^{\circ }$ N)	Lon. (${}^{\circ }$ E)	M
EQ1	12 July 1993	Southwest-Off Hokkaido EQ	42.78	139.18	7.8
EQ2 1	4 October 1994	East-Off Hokkaido EQ	43.38	147.67	8.2
EQ3	26 September 2003	Off Tokachi EQ	41.78	144.08	8.0
EQ4	22 December 2010	Near Chichi-jima EQ	27.05	143.93	7.8
EQ5	11 March 2011	Tohoku EQ	38.10	142.86	9.0

¹ This EQ lies outside of the smaller area ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{148}$.

). We first identify the $\beta$ minima that appear a few months before all these EQs and determine their ${\beta }_{300,min}/{\beta }_{200,min}$ values (see

Table 2. The values of ${\beta }_{W,min}$, ${\alpha }_{min,bef}$, ${\alpha }_{corr}$ and ${\alpha }_{min,aft}$ in the larger area investigated with sliding natural time windows of length 6 to W that preceded the $M_{JMA}\ge 7.8$ EQs listed in Table 1. Hereafter, the value of ${\alpha }_{corr}$ is given when ${\beta }_{250,min}$ appears and ${\alpha }_{min,bef}$ is the minimum of the DFA exponent up to three-and-a-half months (105 days) before ${\beta }_{250,min}$. Reproduced under permission from Ref. [44]. ©2014. American Geophysical Union. All Rights Reserved.

Label	${\mathit{\beta }}_{\mathit{200},\mathit{min}}$	${\mathit{\beta }}_{\mathit{250},\mathit{min}}$	${\mathit{\beta }}_{\mathit{300},\mathit{min}}$	$\frac{{\mathit{\beta }}_{\mathit{300},\mathit{min}}}{{\mathit{\beta }}_{\mathit{200},\mathit{min}}}$	${\mathit{\alpha }}_{\mathit{min},\mathit{bef}}$	${\mathit{\alpha }}_{\mathit{corr}}$	${\mathit{\alpha }}_{\mathit{min},\mathit{aft}}$
EQ1	0.292(19930523) 1	0.266(19930602)	0.270(19930607)	0.92	0.397(19930327)	0.515(19930602)	0.469(19930613)
EQ2	0.294(19940630)	0.274(19940722)	0.312(19940722)	1.06	0.451(19940627)	0.557(19940722)	0.476(19940724)
EQ3	0.289(20030703) 2	0.273(20030703)	0.307(20030714)	1.06	0.356(20030325)	0.553(20030703)	0.384(20030925) 3
EQ4	0.233(20101130)	0.221(20101130)	0.242(20101130)	1.04	0.467(20101116)	0.501(20101130)	0.385(20101201)
EQ5	0.156(20110105)	0.160(20110108)	0.154(20110105)	0.99	0.348(20101222)	0.524(20110108)	0.422(20110123)

¹ A comparable (within ± 0.0005) value ${\beta }_{200,min}$ = 0.293 is also found on 24 May 1993. ² The same (within ± 0.0005) value ${\beta }_{200,min}$ = 0.289 is also found on 21 June 2003. ³ The date of the M8 EQ is either 25 or 26 September 2003 depending on the use of LT or UT, respectively. Here, we use LT. This value of a is observed 15.5 h before the mainshock.

). These values are found to lie in a narrow range very close to unity [44], i.e., in the range 0.92–1.06. (This range slightly differs from the previously reported [36] range 0.95–1.08 mentioned in Section 3 since in [44], the numerical accuracy of the calculated ${\kappa }_1$ values for $W>100$ is somewhat improved.) This is understood in the context that these values correspond to similar critical processes, thus exhibiting the same dependence of ${\beta }_W$ on W. During the whole period studied, however, beyond the above-mentioned $\beta$ minima before all the $M_{JMA}\ge 7.8$ EQs, more minima exist. Among these minima, we choose those which are as deep or deeper than the shallowest one of the ${\beta }_{200,min}$ values previously identified (e.g., 0.294 in Table 2) and, in addition, they have ${\beta }_{300,min}/{\beta }_{200,min}$ values lying in the narrow range determined above. Thus, we now find a list of “additional” $\beta$ minima (see Table 3 of Ref. [44]) for which is must be checked whether they are either non-precursory or may be followed by EQs of smaller magnitude.

We now repeat the whole procedure—as described above—for the determination of $\beta$ minima in the smaller area. Thus, we obtain a new set of $\beta$ minima (with shallowest ${\beta }_{200,min}$ = 0.293 and ${\beta }_{300,min}/{\beta }_{200,min}$ range 0.97–1.09, see

Table 3. The values of ${\beta }_{W,min}$, ${\alpha }_{min,bef}$, ${\alpha }_{corr}$ and ${\alpha }_{min,aft}$ in the smaller area ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{146}$ investigated with sliding natural time windows of length 6 to W that preceded $M_{JMA}\ge 7.8$ EQs. Reproduced under permission from Ref. [44]. ©2014. American Geophysical Union. All Rights Reserved.

Label	${\mathit{\beta }}_{\mathit{200},\mathit{min}}$	${\mathit{\beta }}_{\mathit{250},\mathit{min}}$	${\mathit{\beta }}_{\mathit{300},\mathit{min}}$	$\frac{{\mathit{\beta }}_{\mathit{300},\mathit{min}}}{{\mathit{\beta }}_{\mathit{200},\mathit{min}}}$	${\mathit{\alpha }}_{\mathit{min},\mathit{bef}}$	${\mathit{\alpha }}_{\mathit{corr}}$	${\mathit{\alpha }}_{\mathit{min},\mathit{aft}}$
EQ1	0.288(19930526)	0.264(19930607)	0.278(19930607)	0.97	0.374(19930319)	0.557(19930607)	0.473(19930617)
EQ3 1	0.293(20030624)	0.311(20030709)	0.309(20030718)	1.06	0.374(20030408)	0.570(20030709)	0.434(20030925) 2
EQ4	0.228(20101130)	0.233(20101130)	0.250(20101130)	1.09	0.467(20101122)	0.516(20101130)	0.408(20101209)
EQ5	0.156(20110105)	0.160(20110108)	0.154(20110105)	0.99	0.354(20101222)	0.522(20110108)	0.439(20110117)

¹ The epicenter of EQ2—not included here—lies outside of the smaller area. ² The date of the M8 EQ is either 25 or 26 September 2003 depending on the use of LT or UT, respectively. Here, we use LT. This value of α is observed 24.5 h before the mainshock.

) that appear a few months before all the ${\mathrm{M}}_{JMA}\ge 7.8$ EQs (cf. there exist four such EQs in the smaller area, see Figure 1), as well as a new list of “additional” $\beta$ minima (see Table 5 of [44]) to be checked for whether they are non-precursory or may be followed by smaller magnitude EQs. Comparing these new $\beta$ minima with the previous ones, we investigate whether they:

(1) Appear in both areas practically on the same date (differing by no more than 10 days or so);

(2) Exhibit the three main features (i.e., the sequence (A) anti-correlated behavior/(B) correlated/(C) almost random behavior) identified from the results of DFA of EQ magnitude time series discussed in the main text of [44]. In Tables 3 and 5 of [44], corresponding to the “additional” minima, are marked in bold the values which do not satisfy at least one of the three inequalities given below which quantify these three main features (A), (B), and (C):

$${\alpha }_{min,bef}\le 0.47,{\alpha }_{corr}>0.50,{\alpha }_{min,aft}\le 0.50$$

In view of the fact that considerable errors are introduced in the estimation of the $\alpha$ exponent when using a relatively small number of points as the one used here, i.e., 300 (for application in other seismically active areas see [66,70]), we adopted the value $\alpha =0.47$ as the maximum one in order to assure anti-correlated behavior in the period before the appearance of $\beta$ minimum.

A summary of the main results obtained after carrying out this investigation is as follows:

First, the $\beta$ minima identified a few months before all the $M_{JMA}\ge 7.8$ EQs with epicenters inside of both areas obey the aforementioned requirements (1) and (2), see Table 2 and Table 3. Remarkably, in the remaining case, i.e., the East-Off Hokkaido M8.2 EQ labeled EQ2 in Table 1, with an epicenter outside of the smaller area but inside the larger, we find $\beta$ minimum not only in the study of the larger area (on 30 June 1994, see Table 2), but also in the relevant study of the smaller area (see the third $\beta$ minimum in Table 5 of [44] observed on 5 July 1994). Second, all the other “additional” $\beta$ minima resulting from the studies in both areas (see Tables 3 and 5 of [44]) violate at least one of the requirements (1) and (2). In other words, only the $\beta$ minima appearing a few months before the five $M_{JMA}\ge 7.8$ EQs obey all the aforementioned properties. Interestingly, this procedure could have been applied before the occurrence of the M9 Tohoku EQ, just after the identification of the deepest $\beta$ minimum on 5 January 2011 with ${\beta }_{300,min}/{\beta }_{200,min}$ almost in unity, thus lying inside the narrow ranges determined from previous $M_{JMA}\ge 7.8$ EQs (see Table 2 and Table 3), resulting in the conclusion that a $M_{JMA}\ge 7.8$ EQ was going to occur in a few months.

We clarify that the above procedure does not preclude of course that one of the “additional” minima may be of truly precursory nature, but corresponding to an EQ of magnitude smaller than the adopted threshold. As a first example, we mention the case marked FA7 in Table 3 of [44] referring to a ${\beta }_{200,min}$ observed on 12 April 2000. It preceded the M6.5 EQ that occurred on 1 July 2000 of the volcanic-seismic swarm activity in the Izu Island region in 2000 (which was preceded by a pronounced SES activity that started on 26 April 2000 recorded by Uyeda and coworkers, as already mentioned in Section 4), see the ${\beta }_{200}$ curve (red) in Figure 3C of [44]. A second example is the case marked EQc in Table 3 of [44], which refers to a ${\beta }_{200,min}$ on 15 October 1994 that preceded the M7.6 Far-Off Sanriku EQ on 28 December 1994 [36].

8. On a Unique Change of the Seismicity Order Parameter Fluctuations before the M9 Tohoku Earthquake in Japan in 2011

In Section 3, upon analyzing the Japan seismic catalog in natural time from 1 January 1984 to 11 March 2011, we computed the seismicity order parameter fluctuations $\beta$ by using a sliding natural time window comprising the number i of EQs that would occur on average within the crucial scale [43] of a few months, or so, which is the average lead time of SES activities. We explained that these fluctuations $\beta$ exhibited distinct minima ${\beta }_{min}$ a few months before all the shallow EQs of magnitude 7.6 or larger in the Japanese area ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{148}$. Among these minima, the minimum before the M9 Tohoku EQ observed at around 5 January 2011 was the deepest. In addition to this unprecedented seismicity order parameter fluctuation decrease, almost two weeks earlier, i.e., on 22 December 2010, a unique increase of the seismicity order parameter fluctuations was observed [71] in the following respect: Upon using a sliding natural time window comprising 300 events, which is the average number of EQs ($M\ge 3.5$) that occur within ≈2 months and plotting the fluctuations $\beta$ of the order parameter of seismicity in the entire Japanese region ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{148}$ versus the conventional time from 1 January 1984 until the Tohoku EQ occurrence on 11 March 2011 upon considering all $M\ge 3.5$ EQs irrespective of their depth, we find that there exist six prominent fluctuations of $\beta$ (greater than 1.2) before the M9 Tohoku EQ. These are accompanied by the following major EQs: Upon the occurrence of the M7.5 EQ on 15 January 1993 at 42.92${}^{\circ }$ N 144.35${}^{\circ }$ E as well as upon the occurrences of the shallow EQs of magnitude 7.6 or larger (EQa to EQe) already listed in Section 3 (see Figure 11 and Figure 12). Among these six prominent $\beta$ fluctuations, the following two are the largest and were observed in 2003 and in 2010: First, the large fluctuation of $\beta$ in 2003 appears upon the occurrence of the M8 Off Tokachi EQ on 26 September 2003 (attaining its maximum value on 28 September 2003, see Figure 12b). It is the largest one observed upon any other major EQ occurrence during the 27-year period of our study. Second, the $\beta$ fluctuation increase on 22 December 2010 is observed upon the occurrence on the same day of the M7.8 near Chichi-jima EQ. Concerning these two larger $\beta$ fluctuations, the following comments apply: The one in 2010 appears almost simultaneously with the minimum $\Delta S_{min}$ of the change $\Delta S$ of the entropy of seismicity (cf. that of the entire Japanese region) under time reversal, which will be discussed in detail in a later part of this work, showing that such a minimum is of precursory nature signaling that a large EQ is impending according to the conclusions deduced from the natural time analysis of the OFC model summarized in [71].

The above facts corroborate the suggestion that between the two larger $\beta$ fluctuations discussed here, only the second, i.e., the one on 22 December 2010, that precedes the Tohoku M9 EQ by almost three months, is likely to be of a precursory nature. It is very important to draw attention to the fact that the validity of this conclusion has been also checked, beyond the case of $i=$ 300 events presented in Figure 11, for several other lengths from $i=$ 150 to 500 events. Such a change of the length i also reveals the following behavior: Upon increasing i, it is observed (see Figures 2B and 4E of Ref. [36]) that the increase $\Delta {\beta }_i$ of the ${\beta }_i$ fluctuation on 22 December 2010 becomes distinctly larger, which does not happen (see Figure 4A–D of [36]) for the $\beta$ fluctuations increase upon the occurrences of all other shallow EQs in Japan of magnitude 7.6 or larger from 1 January 1984 to the time of the M9 Tohoku EQ. Such a behavior that obeys the interrelation $\Delta {\beta }_i=0.5ln(i/114.3)$, see Figure 12g,h, has a functional form strikingly reminiscent of the one discussed by Penrose et al. [73] in computer simulations of phase separation kinetics using the ideas of Lifshitz and Slyozov [74], see their Equation (33) which is also based on work by Lifshitz and Slyozov (cf. the seminal work by Lifshitz and Slyozov [74] and independently by Wagner [75] on phase transitions, which will be explained in more detail in Section 10). It is also very important to note that this functional form is lost in our calculation upon considering not all $M\ge 3.5$ EQs, but restricting ourselves to depths h either $h\le 70$ km or $h\le 300$ km, see the symbols in green in Figure 13. Hence, the $\beta$ fluctuation on 22 December 2010 accompanying the minimum $\Delta S_{min}$ is unique. By employing the most recent method of ECA [47], which considers a time lag $\tau$ and a window $\Delta T(>0)$ between the precursor and the event to be predicted, the profound statistical significance of this unique result can be further assured as follows: Assuming the M9 Tohoku EQ as the event to be predicted, we estimated the probability (p-value) that the fluctuation increase on 22 December 2010 is a chancy precursor with a lag $\tau$ and a window $\Delta T$, i.e., the EQ occurs within the time period from $\tau$ to $\tau +\Delta T$ days after. By employing the function CC.eca.es of the CoinCalc package [76] for R that implements ECA, we found that the probability of this happening by chance is very small, i.e., it varies from $p=$ 0.79% to 0.01% upon varying the window $\Delta T$ varies from 78 days to 1 day, respectively.

In short, on 22 December 2010, the order parameter fluctuations (that are minimized almost two weeks later) showed a unique increase exhibiting a functional form, which stems from the ideas of Lifshitz and Slyozov for the time growth of the characteristic size of the minority phase droplets in phase transitions. Remarkably, this unique increase occurred at the same date at which the entropy change of seismicity under time reversal was minimized as shown in [77], which will be discussed in the next section of this work.

9. Natural Time Analysis: Minimum of the Seismicity Entropy Change under Time Reversal before Major Earthquakes

In [77], it was shown that upon analyzing in natural time all EQs of magnitude M 3.5 or larger in Japan from 1 January 1984 until the occurrence of the super-giant M9 Tohoku EQ on 11 March 2011, a pronounced minimum of the entropy change of seismicity under time reversal is observed two-and-a-half months before this M9 EQ. Notably, a simultaneous minimum with an unusual low value ($\alpha \sim 0.35$) for the exponent $\alpha$ resulting from the detrended fluctuation analysis of the EQ magnitude time series indicates an evident anticorrelated behavior. These findings are consistent with an earlier result [78] that upon analyzing the OFC model for EQs in natural time, which is the most studied non-conservative self-organized criticality model for EQs, one finds a non-zero change of the entropy upon time reversal. This reveals a breaking of the time symmetry, and hence reflects the predictability in this model. Note that OFC model originated from a simplification of the Burridge–Knopoff spring-block model [79].

We recall from Section 3 that for a time series comprising N events, an index for the occurrence of the k-th event is defined by ${\chi }_k=k/N$, which is termed natural time. We, then, study the pairs $({\chi }_k,Q_k)$, or the pairs $({\chi }_k,p_k)$, where $p_k=Q_k/{\sum }_{n=1}^N{Q}_n$ is the normalized energy for the k-th event. The entropy S in natural time is defined [8,80] as

$$S=\langle \chi ln\chi \rangle -\langle \chi \rangle ln\langle \chi \rangle ,$$

(7)

where the bracket $\langle f\left(\chi \right)\rangle ={\sum }_{k=1}^N{p}_{k}f\left({\chi }_k\right)$ denotes the average value of $f\left(\chi \right)$ weighted by $p_k$, i.e., $\langle \chi ln\chi \rangle ={\sum }_{k=1}^N{p}_k(k/N)ln(k/N)$ and $\langle \chi \rangle ={\sum }_{k=1}^N{p}_k(k/N)$. It is dynamic entropy depending on the sequential order of events [81], thus changing upon the occurrence of each event. The entropy obtained by Equation (7) upon considering [9] the time-reversal $\widehat{T}$, i.e., $\widehat{T}{p}_k=p_{N-k+1}$, is labeled by $S_{-}$, i.e.,

$$\begin{array}{c}\hfill S_{-}=\sum _{k=1}^N{p}_{N-k+1}\frac{k}{N}ln\left(\frac{k}{N}\right)-\left(\sum _{k=1}^N{p}_{N-k+1}\frac{k}{N}\right)ln\left(\sum _{k=1}^N{p}_{N-k+1}\frac{k}{N}\right)\end{array}$$

(8)

(see also [56,57]).

The difference $S-S_{-}$ labeled $\Delta S$ may also have a subscript ($\Delta S_i$) meaning that the calculation is made (for each S and $S_{-}$) with a sliding natural time window of length i (=number of successive events), i.e., at scale i (see also below). As for the calculation of the time series of the entropy change under time reversal, a window of length i is sliding each time by one event, through the whole time series and the entropies S and $S_{-}$ (and therefrom their difference $\Delta S_i$) are calculated each time. Thus, we form a new time series comprising successive $\Delta S_i$ values. It was found [1] that $\Delta S_i$ is a key measure which may identify when the system approaches the critical point (dynamic phase transition). $\Delta S_i$ has been applied [10,82], for example, for the identification of the impending sudden cardiac death risk, which is the major cause of death in industrialized countries [83,84,85]. In addition, $\Delta S_i$ was employed as a useful tool [78] (see also Section 8.3.4 of [1]) to study the predictability of the OFC model, as mentioned above.

The analysis in natural time was conducted by using the JMA seismic catalog, e.g., see [36,37] and considering all EQs of magnitude $M\ge 3.5$ from 1984 until the Tohoku EQ occurrence on 11 March 2011 within the area ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{146}$ (see the map of Figure 3). The calculation was performed also for a second larger area ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{148}$ in order to avoid boundary effects and assure that the results do not depend on the selection of the area studied having on average ≈125 and ≈145 EQs per month for the smaller and larger area, respectively. The energy of EQs was obtained from the JMA magnitude $M_{JMA}$ after converting [39] to the moment magnitude $M_w$ [40].

The time evolution of $\Delta S_i$ was studied for a number of scales i of the seismicity with $M_{JMA}\ge 3.5$ occurring in both areas, larger and smaller, by selecting proper scales i as follows: An SES activity, whose lead time ranges from a few weeks up to around $5\frac{1}{2}$ months [1], exhibiting critical behavior [3,8,86], is observed during a period in which long-range correlations between EQ magnitudes prevail. On the other hand, before the initiation of the SES activity, and hence before ${\beta }_{min}$, another stage appears in which the temporal correlations between EQ magnitudes exhibit an anticorrelated behavior [44]. Hence, there exists a significant change in the temporal correlations between EQ magnitudes when comparing the two stages that correspond to the periods before and just after the initiation of an SES activity. This change is likely to be captured by the time evolution of $\Delta S_i$, thus, our study of $\Delta S_i$ started from the scale of $i\sim {10}^3$ events, which corresponds to the number of seismic events $M_{JMA}\ge 3.5$ that occur during a period around the maximum lead time of SES activities.

Plotting the $\Delta S_i$ values versus the conventional time for several scales, a careful inspection of the results depicted in Figures 2 and 4 of [77] for the larger and smaller area, respectively, reveals the following common feature: At shorter scales, i.e., from $i={10}^3$ to $3\times {10}^3$ events, a number of local minima appear, but leaving aside all these changes, we find that at longer scales, i.e., $4\times {10}^3$, $5\times {10}^3$ and $7\times {10}^3$ events, a pronounced minimum is observed on 22 December 2010 (on this date, the M7.8 Near Chichi-jima EQ occurred at $27.{05}^{\circ }$ N $143.{94}^{\circ }$ E [36,37]). It was shown that this minimum is statistically significant, for example in the larger area, by means of the following procedure: By randomly shuffling the EQ magnitude time series and assigning each magnitude to an existing EQ occurrence time, the calculations were repeated ${10}^2$ times and investigated the resulting $\Delta S_i$ time series of the longer scales, i.e., $\Delta S_{4000}$, $\Delta S_{5000}$ and $\Delta S_{7000}$, for minima occurring on the same date and deeper than or equal to those depicted in the aforementioned Figure 2 of Ref. [77]. Only three such cases were found out of the ${10}^2$ studied. Hence, the probability of obtaining minima comparable or deeper than those shown in Figure 2 of Ref. [77] by chance is approximately 3%, which demonstrates that our result is statistically significant.

This minimum on 22 December 2010 was found [77] to be robust upon changing either the EQ depth (from 0–300 km) or the magnitude threshold (considering also additional thresholds $M_{thres}=$ 3.7 and 4.0) as well as the size of the area, i.e., beyond the two areas 21${}^{\circ }\times {23}^{\circ }$ (larger, ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{148}$) and 21${}^{\circ }\times {21}^{\circ }$ (smaller, ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{146}$) studied, we repeated the calculations for thirteen additional areas.

In summary, the earlier finding in [78] that before a large avalanche in the OFC model for EQs a minimum of the entropy change of seismicity under time reversal is observed is confirmed. This is so because the natural time analysis of seismicity in Japan—during the almost 27-year period from 1 January 1984 until the occurrence of the M9 Tohoku on 11 March 2011—revealed that for longer scales, i.e, $i>3500$ events, the minimum of $\Delta S_i$ values was observed on 22 December 2010.

10. Fluctuations of the Entropy Change of Seismicity under Time Reversal before the M9 Tohoku Earthquake in Japan

In summary, in [87], it was shown that almost two-and-a-half months before the M9 Tohoku EQ that occurred in Japan on 11 March 2011, the complexity measure ${\Lambda }_i$ associated with the entropy change of seismicity in natural time under time reversal (that will also be defined below) exhibited an abrupt increase, which conforms to the Lifshitz–Slyozov–Wagner (LSW) theory for phase transitions (mentioned in short in Section 8). In particular, ${\Lambda }_i$ displayed a behavior similar to that of the characteristic size of the minority phase droplets through a scaling behavior in which time growth has the form $A{(t-t_0)}^{1/3}$ where the prefactors A are proportional to the scale i, while the exponent (1/3) is independent of i. As for the Tsallis entropic index q [88,89,90,91,92,93], it exhibits a simultaneous increase which interestingly obeys the same exponent (1/3), but the prefactors A are not proportional to the scale i. In addition, it was found that the complexity measures ${\Lambda }_{2000}$, ${\Lambda }_{3000}$ and ${\Lambda }_{4000}$ show a decrease just after the fulfillment of the criticality condition ${\kappa }_1=0.070$ almost one day before the Tohoku EQ occurrence.

Beyond the aforementioned M9 Tohoku EQ in 2011, in Japan, another recent example was studied in [94], which is the M8.2 EQ that occurred in the Chiapas region in Mexico on 7 September 2017. This was Mexico’s largest EQ in more than a century. Upon analyzing the seismicity in natural time, it was shown that almost three months before its occurrence, the following precursory behavior was identified: The entropy change under time reversal exhibited a minimum [94] together with increased fluctuations of the entropy change under time reversal, as well as by a simultaneous Tsallis entropic index q increase [95]. The variations of this index before strong EQs has been studied by several researchers [96,97,98,99,100,101,102,103]. In addition, by analyzing the seismicity during the 6-year period 2012–2017 in natural time in the Chiapas region where the M8.2 EQ occurred, it was found [95] that on the same date, i.e., on 14 June 2017, the complexity measure ${\Lambda }_i$ (see below) associated with the fluctuations of the entropy change under time reversal showed an abrupt increase together with a simultaneous increase of the Tsallis entropic index q [88,89,90,91,92,93].

As already explained in the preceding section, the entropy S in natural time, which is a dynamic entropy that exhibits [9] concavity, positivity and Lesche stability [104,105], is given by:

$$S\equiv \langle \chi ln\chi \rangle -\langle \chi \rangle ln\langle \chi \rangle$$

(9)

where $\langle f\left(\chi \right)\rangle ={\sum }_{k=1}^N{p}_{k}f\left({\chi }_k\right)$. Its value upon considering the time reversal $\widehat{T}$, i.e., $\widehat{T}{p}_k=p_{N-k+1}$, is labeled by $S_{-}$. The value of $S_{-}$ differs, in general, from S, see, e.g., [10], (see also [57] and references therein), and thus, S does satisfy the conditions to be “causal”. The physical meaning of the change of entropy $\Delta S\equiv S-S_{-}$ in natural time under time reversal is discussed in [1,10,57].

Using a window of length i (number of consecutive events) sliding through the time series of L consecutive EQs, the entropy in natural time has been determined for each position $j=1,2,\dots ,L-i$ of the moving window. Thus, a time series of $S_i$ is obtained. By considering the standard deviation $\sigma (\Delta S_i)$ of the time series of $\Delta S_i\equiv S_i-{\left(S_{-}\right)}_i$, we define [1,95,106], the complexity measure ${\Lambda }_i$

$${\Lambda }_i={\displaystyle \frac{\sigma (\Delta S_i)}{\sigma (\Delta S_{100})}}$$

(10)

when a moving natural time window comprising i consecutive events is sliding through the time series and the denominator has been selected [95] to correspond to the standard deviation $\sigma (\Delta S_{100})$ of the time series of $\Delta S_i$ of i = 100 events. This complexity measure quantifies how, upon changing the scale from ${10}^2$ events to a longer scale, e.g., $i={10}^3$ events, the statistics of the $\Delta S_i$ time series changes. The calculations are performed by means of a window of length i (=number of successive EQs) sliding, each time by one EQ. The difference $\Delta S_i$ of the entropies S and $S_{-}$ is calculated each time, hence, a new time series comprising successive $\Delta S_i$ values is formed and the complexity measure ${\Lambda }_i$ is determined in accordance with its definition in Equation (10).

Here, in Figure 14a–f, the ${\Lambda }_i$ values, for example for the scales i = 2000, 3000 and 4000 events, are plotted versus the conventional time from 1 January 1984 until the occurrence of the M9 Tohoku EQ on 11 March 2011. Furthermore, all EQs of magnitude 7.0 or larger are also inserted in the same figure with vertical lines ending at circles, which are read in the right scale. This figure reveals that, for all the three scales studied, almost $2\frac{1}{2}$ months before the M9 EQ, i.e., on 22 December 2010, we observe an abrupt increase of the ${\Lambda }_i$ values. This happens upon the M7.8 EQ occurrence with an epicenter at $27.{05}^{\circ }$ N $143.{94}^{\circ }$ E [37]. This can be better seen in Figure 15, which is a three-month excerpt of Figure 14f in expanded horizontal scale from 10 December 2010 until 11 March 2011 and depicts for readers’ better inspection the increase $\Delta {\Lambda }_i$ of the ${\Lambda }_i$ values after the M7.8 EQ on 22 December 2010 mentioned above.

Figure 15 reveals that all three complexity measures, ${\Lambda }_{2000}$, ${\Lambda }_{3000}$ and ${\Lambda }_{4000}$ (i.e., the ${\Lambda }_i$ values at the scales i = 2000, 3000 and 4000 events, respectively), show an obvious and abrupt increase on 22 December 2010 and after the aforementioned M7.8 EQ occurrence exhibit a scaling behavior of the form

$$\Delta {\Lambda }_i=A{(t-t_0)}^c$$

(11)

where the exponent c is approximately equal to $1/3$ independently of the scale i, while the pre-factors A are proportional to i (see Figure 16) and $t_0$ is almost 0.2 days after the M7.8 EQ occurrence. Equation (11) confirms the seminal work by Lifshitz and Slyozov [74] and independently by Wagner [75] on phase transitions (LSW theory), which predicts that the time growth of the characteristic size of the minority phase droplets grows with time as $t^{1/3}$.

To shed more light on whether this increase of ${\Lambda }_i$ is possibly associated with strong EQs in general, let us now further investigate the obvious increases (which are the most significant ones in Figure 14) of ${\Lambda }_i$ on 2 November 1989 and on 15 January 1993 that can be visualized in Figure 14b and study whether they exhibit a scaling behavior similar to that identified above for the abrupt ${\Lambda }_i$ increase on 22 December 2010. Concerning the ${\Lambda }_i$ increase on 2 November 1989, Figure 17 depicts the log–log plot of the changes $\Delta {\Lambda }_i$ of the complexity measures ${\Lambda }_{2000}$, ${\Lambda }_{3000}$ and ${\Lambda }_{4000}$ versus the elapsed time $(t-t_0)$ in days from the M7.1 EQ occurrence on the same date increased by 0.022 days. In addition, for the readers’ convenience, the straight black line with slope 0.5 is also plotted in Figure 17 in order to easily recognize that the value of c for the three complexity measures, ${\Lambda }_{2000}$, ${\Lambda }_{3000}$ and ${\Lambda }_{4000}$, is close to 0.5, thus differing markedly from the value $c=1/3$ predicted by the LSW theory. Furthermore, a careful inspection of Figure 17 reveals that the pre-factors are not proportional to i (e.g., see that for i = 3000 events, the green line is higher than the blue line that corresponds to i = 4000 events). Hence, the ${\Lambda }_i$ increase associated with the EQ on 2 November 1989 does not obey Equation (11), thus strongly deviating from the LSW theory of phase transitions. Along the same lines, we now give in Figure 18, as in Figure 17, the log–log plot that corresponds to the $\Delta {\Lambda }_i$ changes of the complexity measures ${\Lambda }_{2000}$, ${\Lambda }_{3000}$ and ${\Lambda }_{4000}$ versus the elapsed time $(t-t_0)$ in days from the M7.5 EQ occurrence on 15 January 1993 increased by 0.014 days. The exponent in Figure 18 is also close to 0.5—thus deviating from the value $c=1/3$ of LSW theory—and the relevant prefactors A are not proportional to the scale i, for example, the red crosses corresponding to the scale i = 2000 do not practically differ from the green symbols of the larger scale i = 3000.

We now proceed to the results deduced by means of nonextensive statistical mechanics [91], pioneered by Tsallis [88,93]. This has found application [89,90] in the physics of EQs through the Tsallis entropic index q. On the basis of the EQ magnitude distribution, one can obtain [92,107] the entropic index q and study how it varies with time as we approach a strong EQ (for a recent review, see [103]). In Figure 19, we present the q-value versus conventional time during 1984–2011 as it was estimated [107] for several sliding windows of i = 1000, 2000, 3000, 4000 and 5000 consecutive EQs for $M_{JMA}\ge 3.5$ in the Japanese area ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{148}$. We find that the q-value shows an abrupt increase upon the M7.8 EQ occurrence on 22 December 2010 and hence before the occurrence of the M9 Tohoku EQ. A three-month excerpt from Figure 19 is shown in expanded time scale in Figure 20 from 10 December 2010 until the Tohoku EQ on 11 March 2011. We finally present in Figure 21 as in Figure 16, Figure 17 and Figure 18, a log–log plot of the changes $\Delta q$ of the Tsallis entropic index versus the time $(t-t_0)$ in days elapsed from the M7.8 EQ occurrence on 22 December 2010 increased by approximately 0.2 days. Figure 21 reveals that the exponent c in Equation (11) is around $1/3$, in accordance with the LSW theory, but the prefactors A are not proportional to i.

Figure 15 also shows the following very interesting feature of the complexity measure ${\Lambda }_i$: On 9 March 2011, a M7.3 foreshock occurred, and the next day, on 10 March 2011, all three complexity measures, ${\Lambda }_{2000}$, ${\Lambda }_{3000}$ and ${\Lambda }_{4000}$, exhibited a simultaneous decrease. What is the origin of this ${\Lambda }_i$ decrease since the next day, i.e., on 11 March 2011, the M9.0 mainshock occurred? To answer this question, we recall Figure 8 (and in particular Figure 8B) which reveals that the condition ${\kappa }_1=0.070$ (which signals that the mainshock is going to occur shortly, i.e., within the next few days or so) is fulfilled for $M_{thres}$ = 4.2 to 5.0 in the morning of 10 March 2011 upon the occurrence of the EQs from 08:36 to 13:14 LT, i.e., almost one day before the Tohoku EQ. Strikingly, we now see clearly that the aforementioned decrease of $\Delta {\Lambda }_{2000}$, $\Delta {\Lambda }_{3000}$ and $\Delta {\Lambda }_{4000}$ after the M7.3 EQ occurs just after the gray shaded area in Figure 8B where the criticality condition ${\kappa }_1=0.070$ was obeyed, thus signaling that the mainshock was approaching.

11. Identifying the Occurrence Time of a Mainshock by Means of the Fluctuations of the Seismicity Entropy Change under Time Reversal

It is in the scope of this section to explain a procedure to identify the approach of the critical point (mainshock) without making use of the criticality condition ${\kappa }_1=0.070$ for the evolution of seismicity after the initiation of the SES activity in the candidate epicentral area that was used in Section 6. In particular, we shall show here, following [108], that in natural time analysis, the complexity measure ${\Lambda }_i$ that quantifies the fluctuations of the entropy change $\Delta S$ of the seismicity in the Japanese area under time reversal plays a key role in identifying the occurrence time of the Tohoku mega-EQ that occurred on 11 March 2011 in Japan with magnitude 9.0, which is the largest-magnitude event recorded in Japan.

The main point of the physical basis for the procedure followed here is that, as mentioned in the first part of Section 9, when analyzing in natural time (see below) the OFC model for EQs [55], $\Delta S$ exhibits an evident minimum [1] (or maximum upon defining, e.g., see [78] $\Delta S\equiv S_{-}-S$ instead of $\Delta S\equiv S-S_{-}$) before a large avalanche, i.e., a large EQ. We also consider the following two points: First, we recall from Section 9, that the study of the time evolution of $\Delta S_i$ should be conducted at scales of the order of $i\sim {10}^3$ events, since this scale corresponds to the number of $M_{JMA}\ge 3.5$ EQs that occur during a period of at least around the maximum lead time of SES activities (because ∼145 EQs occur per month in the entire Japanese region and consider a period of at least $5\frac{1}{2}$ months). Second, the following characteristic feature was identified in [109], where we analyzed the seismic data ($M_{JMA}\ge 3.5$) in the entire region of Japan in natural time and calculated the complexity measure ${\Lambda }_i$: When the calculation is made at scales markedly longer than 2000 events, e.g., $i=$ 3000, 4000 and 5000 events, the ending dates of the calculation show a tendency to be clearly clustered into two groups: First, the one group that comprises markedly larger ${\Lambda }_i$ values corresponds to ending dates later than the date 22 December 2010 at which the minimum of the entropy change of seismicity under time reversal $\Delta S_{min}$ has been observed (as discussed in Section 9), thus being closer to the occurrence date of the Tohoku EQ. Second, the other group that comprises appreciably lower ${\Lambda }_i$ values corresponds to earlier dates. Practically the same behavior is observed upon increasing the magnitude threshold to 4.0 and using scales that are smaller by a factor of 2.5 in view of the smaller number of EQs per month we have for this threshold (cf. Figure S6 in [77]).

In Figure 22, we plot the ${\Lambda }_i$ values versus the scale i (number of events) for all $M\ge 3.5$ EQs in the entire Japanese region ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{148}$. They have been computed by starting the calculation on 1 January 2010 (a) and 1 January 2006 (b) until the following dates: 30 November 2010 (pluses in red, just before a M7.1 EQ on this date), 1 December 2010 (crosses in green), 22 December 2010 (asterisks in blue, just before the M7.8 EQ that occurred on this date), 1 January 2011 (open squares in magenta, this date is very close to the initiation of an SES activity because—as explained in Section 6—anomalous magnetic field variations were recorded during the period of 4 to 14 January 2011 on the z-component at two measuring sites at epicentral distances ≈ 130 km), 1 February 2011 (solid circles in cyan), 1 March 2011 (open circles in olive green), 9 March 2011 (open triangles in orange, at 00:00 LT, thus almost 12 h before the M7.3 EQ on 9 March 2011), 10 March 2011 (gray filled triangles, at 00:00 LT thus almost 12 h after the M7.3 EQ on 9 March 2011), 11 March 2011 (solid circles in black, at 00:00 LT, thus almost 15 h before the mega-earthquake occurrence) and 11 March (inverted red triangles, almost 10 min before the M9 Tohoku EQ).

A careful inspection of this figure reveals that, for scales longer than ≈2000 events, the ${\Lambda }_i$ values after 1 January 2011—a date very close to the initiation of the SES activity (when long-range temporal correlations are developed, as discussed in Section 7—maximize and later start to decrease until the mainshock occurrence (see Figure 22a). This behavior, which is systematic (and robust since it is not affected if we change the scale in the denominator in Equation (10) from $i=100$ events to $i=50$ or 200 events, see Figure S1 in the Supplementary Material of [108]), disappears in the longer scales in Figure 22b when for $M_{JMA}\ge$ 3.5 the calculation starts on 1 January 2006 (the same feature is observed in Figure 23a for $M_{JMA}\ge$ 4.0 EQs, where the computation started on 1 January 2010, while it disappears in Figure 23b when starting in 2006). This reflects that for each scale, when the computation starts appreciably early in the past, a smaller percentage of events close to the mainshock occurrence are included in the ${\Lambda }_i$ calculation.

The above systematic behavior persists upon starting the calculation in 2010 even closer to Tohoku M9 EQ, i.e., on 1 March, or 1 May, or 1 July (instead of 1 January), see Figure 24 for M $\ge 3.5$ EQs, see also Figures 2 and 3 in [110] if we start the calculation on 1 September (a), 1 October (b) and 1 November 2010 (c). In these figures, we can see clearly that ${\Lambda }_i$ decreases gradually when we approach the M9 Tohoku EQ (see also Figure 25). However, when focusing on the ${\Lambda }_i$ values calculated by using only the EQs inside the candidate epicentral region of Figure 1 of [109] (i.e., the blue-green area in Figure 7 estimated, as mentioned, by the spatiotemporal study of $\beta$ developed in [37]), we find (see Figure 5 of [109], which is also reproduced here for readers’ better information, see Figure 26) that independently from the year at which the computation starts (i.e., 1 January 2006, 1 January 2008, 1 January 2009, and 1 January 2010, thus almost 5, 3, 2 and 1 year before the 2011 M9.0 EQ), for every scale larger than around $i=$ 400–500 EQs all ${\Lambda }_i$ values almost coincide, but after the 9 March 2011 M7.3 EQ, a monotonic increase is observed, which continues until 10 minutes before the M9 EQ on 11 March 2011.

Let us summarize concerning the evolution of the ${\Lambda }_i$ values obtained from the natural time analysis of the EQs inside the large area ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{148}$ (entire Japanese region), the following facts have been secured for (each) longer scale i (cf. These results are consistent with the PSPC physical model for SES generation as it will become clear in the next Section 12): The values of ${\Lambda }_i$ start to increase after 22 December 2010, when $\Delta S$ reached a minimum, and attain a maximum around the first days of January 2011, very close to the beginning of the strong SES activity that gave rise to the anomalous variations of the magnetic field of the Earth recorded predominantly on the vertical component (in a period when EQ magnitudes exhibit long-range correlations as evidenced by the DFA). Subsequently, ${\Lambda }_i$ decrease gradually up to the 11 March 2011 M9 Tohoku EQ.

Second, ${\Lambda }_i$ calculated by considering the seismicity inside the candidate epicentral area behaves differently, increasing abruptly after the M7.3 EQ on 9 March 2011 until the Tohoku M9 EQ.

The above observation that, when the M7.3 EQ occurred on 9 March 2011, ${\Lambda }_i$ of the seismicity inside the epicentral area exhibited a clear increase, while ${\Lambda }_i$ in the large area ${\mathrm{N}}_{25}^{46}{\mathrm{E}}_{125}^{148}$ continued to diminish, represents the important difference unveiled by natural time analysis leading, well in advance, to the characterization of the M7.3 EQ as a foreshock. This characterization agrees with the conclusion drawn in Section 6 by an independent procedure.

12. Discussion: Compatibility of the SES Generation Model with the Experimental Results of the Recent Decade

As mentioned in the Introduction, here, we present in short the main results of the recent decade by means of natural time analysis and investigate whether the experimental results are compatible with the SES generation model presented in Section 2. We start with the discussion of the results deduced in the preceding Section 11: In Figure 27a,b, we have plotted ${\Lambda }_i$ at various ending dates close to the end of December 2010 and the beginning of January 2011 in order to make it clear that they have reached a maximum value after starting to increase from their minimum on 22 December 2010 when the entropy of the Japanese seismicity $\Delta S_{min}$ exhibited its minimum. This behavior can be understood along the following lines.

Close to the start of 2011, in January, three phenomena appeared. First, the observation of an SES activity initiation that lasted during 4 to 14 January 2011 since Earth’s magnetic field exhibited [59,61,62] anomalous variations mainly on the vertical component at two sites—Esashi (ESA) and Mizusawa (MIZ), see Figure 1 of [109]—located at epicentral distances of around 130 km. In addition, two more phenomena appeared almost simultaneously: The establishment of long-range temporal correlations between EQ magnitudes (see Section 7 and [44]) and the occurrence of the deepest ${\beta }_{min}$ (since 1 January 1984) of the seismicity order parameter fluctuations $\beta$ on 5 January 2011 [36] (see also Section 3). The latter is in agreement with the experimental result that the seismicity order parameter ${\kappa }_1$ fluctuations are observed [21,38] to minimize when an SES activity starts. These three simultaneous phenomena around 5 January 2011 were preceded by a clear anticorrelated behavior in EQ magnitudes since, as mentioned in Section 7, the DFA exponent of the EQ magnitudes time series becomes $\alpha =0.35$ upon the occurrence of a M7.8 EQ at 27.05${}^{\circ }$ N 143.94${}^{\circ }$ E (Southern Japan) on 22 December 2010. Strikingly, on this date, the initially random horizontal GPS azimuths gradually started to orient toward the southern direction probably as a result of excess stress. These observations are fully compatible with the first stage (stage A) of the SES generation model of Section 2 (see also [21]) according to which as a result of an excess stress disturbance (cf. this excess disturbance seems to be associated with the unique increase of the $\beta$ fluctuations discussed in Section 8), $\sigma$ starts to increase until it reaches ${\sigma }_{cr}$ at which the prexisting electric dipoles of the future focal area cooperatively align (stage B of the physical model, see Section 2). Then, very interestingly, the most intense crust uplift was observed [21]. Later than the observation of ${\beta }_{min}$ on 5 January 2011, the GPS azimuth orientations returned [111] to become random, being consistent with an EQ magnitudes time series DFA exponent $\alpha \approx 0.5$ [44], around 13 January 2011.

Subsequently, anti-correlation was observed around 23 January 2011 with an exponent of DFA ${\alpha }_{min}=0.42$ and a shift of EQ-related stress disturbance was identified [111] where westward movements replaced the southward ones, in the sense that orientations of the residual displacements were re-aligned along the western direction and the crust depressed. Later, the stress disturbance gradually approached the threshold of the fault rupture, and the orientations of the residual surface displacements became random again [111], thus agreeing with the results of DFA for the EQ magnitudes which led to $\alpha \approx$ 0.5 up to approximately 10 February 2011 (see Figure 5 of [44]), indicating randomness.

We finally comment on the fact that a behavior close to random was observed before the Tohoku EQ. This is in accordance with findings in time series originating from other complex systems such as electrocardiograms (ECGs). In this case, long-range temporal correlations indicating heart rate variability of healthy subjects cease for those of high risk of sudden cardiac death. The latter is accompanied by the appearance of uncorrelated randomness [1,69] (cf. sudden cardiac death can be considered as a critical phenomenon, e.g., see [10,80,81]). Figure 28 shows that ${\Lambda }_i$ values for H exceed markedly those for SD. Thus, the fact that ${\Lambda }_i$, for the longer scales, first reached maximum in the beginning of January 2011 and later gradually decreased until 10 minutes before the M9 Tohoku EQ occurrence fully agrees with a similar behavior in other complex systems.

Other points explained in the frame of the SES generation physical model can be found in [21], in which a schematic diagram is presented (see Figure 29) that compiles the multidisciplinary changes (during the period of observations from 12 December 2010 until 23 January 2011) before the M9 Tohoku EQ that occurred in Japan on 11 March 2011.

Last but not least, we note that in a recent perspective [110] under the title of “Self-organized Criticality and Earthquake Predictability: A long standing question in the light of natural time analysis”, the following matter is discussed: After the Bak-Tang-Wisenfeld seminal work on self-organized criticality (SOC), the following claim appeared by other workers in the 1990s: Large earthquakes cannot be predicted, since the Earth is in a state of SOC, hence, any small EQ has some probability of cascading into a large event. In that perspective [110], we discussed that such claims do not stand in the light shed by natural time analysis, which was shown in the beginning of 2000s to extract the maximum information possible from complex systems time series, as well as unveiling hidden properties in complex time series.

Let us summarize the most important dates mentioned in this work upon applying natural time analysis to the Japanese seismic data from 1 January 1984 until the super-giant M9 Tohoku earthquake on 11 March 2011. We find that, leaving aside the details, precursory phenomena were mainly accumulated around two dates, i.e., 22 December 2010 and 5 January 2011, which strikingly concur with the two stages A and B of the SES generation PSPC physical model, respectively. These phenomena, accompanying other phenomena of multidisciplinary nature, include the following (see Figure 29):

(A) Around 22 December 2010:

• The entropy change of seismicity under time reversal is minimized along with increased fluctuations;
• There is fluctuation increase of the order parameter of seismicity;
• The DFA exponent diminished in the value ${\alpha }_{min,bef}=$ 0.35, which is the lowest value observed during the period 1984–2011 of our study, pointing to an evident anti-correlated behavior in the EQ magnitude time series;
• The horizontal GPS azimuths started to become gradually oriented toward the southern direction (while during the preceding period 12–22 December 2010, they had random orientations);
• Anomalous changes started in the groundwater (level drop, temperature decrease, and probably increase in radon concentration);
• There was an increase in the Tsallis entropic index q.

(B) Around 5 January 2011:

• Unprecedented minimum ${\beta }_{min}$ of the fluctuations of the order parameter of seismicity;
• Earth’s magnetic field started exhibiting anomalous variations primarily on the vertical component (which, according to Maxwell equations, should be accompanied by a strong SES activity). This agrees with the findings in Greece where EQs of magnitude $M\ge 6.5$ were predated by considerable SES activities with simultaneous variations of the magnetic field of the Earth [31] predominantly observed on the vertical component [32,33];
• Full alignment of the orientations of the GPS azimuths southwards that was accompanied by the most intense crust uplift;
• Long-range temporal correlations in the earthquake magnitude time series.

All the above phenomena were observed to begin and end well before the M9 Tohoku EQ occurrence, as schematically shown in Figure 2a in accordance with the SES generation PSPC model (Figure 1). This differs essentially from other proposed precursory mechanisms [22], which suggested that the anomalous behavior becomes more intense upon approaching the EQ failure, as seen in Figure 2b.

13. Conclusions

The main advances in geosciences during the last decade in natural time analysis are presented. Among these advances are the following:

• It has been experimentally demonstrated that there exists a physical interconnection of an SES activity with seismicity. In particular, the fluctuations of the order paramenter of seismicity exhibit a clearly detectable minimum at the time of the initiation of an SES activity. These two phenomena are also linked closely in space;
• Focusing on the minima of the order parameter fluctuations of seismicity preceding all EQs of magnitude 8 (and 9) class that occurred in the entire Japanese area from 1 January 1984 to 11 March 2011 and applying Detrended Fluctuation Analysis to the EQ magnitude time series, it was found [44] that each of these minima is preceded as well as followed by characteristic changes of temporal correlations between EQ magnitudes (cf. EQ magnitude correlations can be also identified by natural time analysis, see, e.g., Refs. [7,112,113,114,115,116]). The following three main features were identified: The order parameter fluctuation minima are observed during periods when long-range correlations have been developed, while they are preceded by a stage in which an evident anti-correlated behavior appears. After these minima, the long-range correlations break down to an almost random behavior, possibly turning to anti-correlation;
• The epicentral area as well as the occurrence time of an impending major EQ can be identified. A typical example is the case of the M9 Tohoku EQ that occurred in Japan on 11 March 2011;
• To achieve the above identification, the key quantities in natural time analysis are the order parameter of seismicity and the entropy change of seismicity under time reversal along with their fluctuations. The accurate achievement of the identification of these key quantities requires, however, that an area should have been covered by a well-developed seismological network such as Japan, for example;
• When an EQ of magnitude class 7.0 occurs, the crucial question emerges of whether it is a foreshock or a mainshock. This can be answered by means of the key quantities mentioned in the previous conclusion;
• The experimental findings identified before major EQs are consistent with a physical model proposed in the 1980s for the generation of Seismic Electric Signals (critical phenomenon) based on the thermodynamics of point defects in solids.