Recent seismic activity at Cephalonia ( Greece ) : a study through candidate electromagnetic precursors in terms of non-linear dynamics

The preparation process of two recent earthquakes (EQs) that occurred in Cephalonia (Kefalonia), Greece, ((38.22 N, 20.53 E), 26 January 2014, Mw = 6.0, depth ∼ 20 km) and ((38.25 N, 20.39 E), 3 February 2014,Mw = 5.9, depth ∼ 10 km), respectively, is studied in terms of the critical dynamics revealed in observables of the involved non-linear processes. Specifically, we show, by means of the method of critical fluctuations (MCF), that signatures of critical, as well as tricritical, dynamics were embedded in the fracture-induced electromagnetic emissions (EMEs) recorded by two stations in locations near the epicentres of these two EQs. It is worth noting that both the MHz EMEs recorded by the telemetric stations on the island of Cephalonia and the neighbouring island of Zante (Zakynthos) reached a simultaneously critical condition a few days before the occurrence of each earthquake. The critical characteristics embedded in the EME signals were further verified using the natural time (NT) method. Moreover, we show, in terms of the NT method, that the foreshock seismic activity also presented critical characteristics before each event. Importantly, the revealed critical process seems to be focused on the area corresponding to the western Cephalonia zone, following the seismotectonic and hazard zoning of the Ionian Islands area near Cephalonia.


Introduction
The possible connection of the electromagnetic (EM) activity that is observed prior to significant earthquakes (EQs) with the corresponding EQ preparation processes, often referred to as seismo-electromagnetics, has been intensively investigated during recent years.Several possible EQ precursors have been suggested in the literature (Uyeda et al., 2009a;Cicerone et al., 2009;Hayakawa, 2013a, b;Varotsos, 2005;Varotsos et al., 2011b;Hayakawa et al., 2015a, b;Contoyiannis et al., 2016;Potirakis et al., 2016).The possible relation of the field-observed fracture-induced electromagnetic emissions (EME) in the frequency bands of MHz and kHz, associated with shallow EQs with magnitude 6 or larger that oc-Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.S. M. Potirakis et al.: Recent seismic activity at Cephalonia curred inland or near the coast, has been examined in a series of publications in order to contribute to a better understanding of the underlying processes (e.g.Eftaxias et al., 2001Eftaxias et al., , 2004Eftaxias et al., , 2008Eftaxias et al., , 2013;;Kapiris et al., 2004;Karamanos et al., 2006;Papadimitriou et al., 2008;Contoyiannis et al., 2005Contoyiannis et al., , 2013;;Eftaxias and Potirakis, 2013;Potirakis et al., 2011Potirakis et al., , 2012aPotirakis et al., , b, c, 2013Potirakis et al., , 2015;;Minadakis et al., 2012a, b;Balasis et al., 2011Balasis et al., , 2013;;Donner et al., 2015;Kalimeris et al., 2016).Additionally, a four-stage model for the preparation of an EQ by means of its observable EM activity has been recently put forward (Eftaxias and Potirakis, 2013;Contoyiannis et al., 2015, and references therein).In summary, the proposed four stages of the last part of the EQ preparation process and the corresponding EM observations, for which specific features have been identified using appropriate time-series analysis methods, appear in the following order (Donner et al., 2015, and references therein): first stage: valid MHz anomaly; second stage: kHz anomaly exhibiting tricritical characteristics; third stage: strong avalanche-like kHz anomaly; fourth stage: electromagnetic quiescence.It is noted that, according to the aforementioned stage model, the pre-EQ MHz EME is considered to be emitted during the fracture of a part of the Earth's crust that is characterized by high heterogeneity.During this phase the fracture is non-directional and spans a large area that surrounds the family of large high-strength entities distributed along the fault sustaining the system.Note that for an EQ of magnitude ∼ 6, the corresponding fracture process extends to a radius of ∼ 120 km (Bowman et al., 1998).
A pair of simultaneous MHz EM signals was recorded at 41 MHz, with a sampling rate of 1 sample s −1 , by two different stations prior to each one of the above-mentioned significant shallow EQs.On 24 January 2014, 2 days before the M w = 6.0 Cephalonia EQ (EQ1), two telemetric stations of our EM signal monitoring network (see Fig. 1), the station of Cephalonia, located on the same island (38.18 • N, 20.59 • E), and the station of Zante (Zakynthos), located on a neighbouring island belonging to the same (Ionian) island complex (37.77• N, 20.74 • E), simultaneously recorded the first pair of aforementioned signals.The same picture was repeated for the second significant Cephalonia EQ, M w = 5.9 (EQ2).Specifically, both the Cephalonia and Zante stations simultaneously recorded the second pair of aforementioned signals on 28 January 2014, 6 days prior to the specific EQ.Note that it has been repeatedly made clear that all the pre-EQ EME signals, which have been observed by our monitoring network, have been recorded only prior to strong shallow EQs that have taken place on land (or near the coastline); this fact, in combination with the recently proposed fractal geo-antenna model (Eftaxias et al., 2004;Eftaxias and Potirakis, 2013), can explain why these signals succeed in being transmitted in air.This model provides a good reason for the increased possibility of detection of such EM radiation, since a fractal geo-antenna emits significantly increased power, compared to the power that would be radiated by the same source, if a dipole antenna model was considered.It should also be noted that none of the recordings of the other monitoring stations of our network (except for the ones of Cephalonia and Zante) presented critical characteristics before these two specific EQs.
The analysis of the specific EM time series, using the method of critical fluctuations (MCF) (Contoyiannis and Diakonos, 2000;Contoyiannis et al., 2002Contoyiannis et al., , 2013)), revealed critical features, implying that the possibly related underlying geophysical process was at critical state before the occurrence of each one of the EQs of interest.The critical characteristics embedded in the specific time series were further verified by means of the natural time (NT) method (Varotsos et al., 2011a, b;Potirakis et al., 2013Potirakis et al., , 2015)).The presence of the "critical point" during which any two active parts of the system are highly correlated, theoretically even at arbitrarily long distances, in other words when "everything depends on everything else", is consistent with the view that the EQ preparation process during the period that the MHz EMEs are emitted is a spatially extensive process.Note that this process corresponds to the first stage of the aforementioned four-stage model.
Moreover, we analysed the foreshock seismic activity using the NT method; the obtained results indicate that seismicity also presented critical characteristics before each one of the two important events.This result implies that the observed EM anomaly and the associated foreshock seismic activity might be considered "two sides of the same coin".Last but not least, 1 day before the occurrence of EQ2, and 5 days after the corresponding critical EME signal, tricritical characteristics were revealed in the EME recorded by the Cephalonia station.
The remainder of this paper is organized as follows: a brief introduction to the MCF and NT analysis methods is provided in Sect. 2. The analysis of the EME recordings according to these two methods is presented in Sect.3. Section 4 presents the results obtained by the analysis of the foreshock seismic activity using the NT method, while Sect. 5 concludes the paper by summarizing and discussing the findings.

27.5
° E Map with the distribution of stations of the telemetric network that monitors electromagnetic variations in the MHz and kHz bands in Greece, which were operating during the time period of interest.The locations of the Cephalonia and Zante stations are marked by the magenta square and triangle, respectively, while the rest of the remote stations are denoted by red circles and the central data recording server by a blue circle.The epicentres of the two significant EQs of interest are also marked, the first (EQ1, M w = 6.0) by a red cross and the second (EQ2, M w = 5.9) by a green X mark.The electromagnetic variations monitored at the Zante station are 3 kHz north-south, 3 kHz east-west, 3 kHz vertical, 10 kHz north-south, 10 kHz east-west, 10 kHz vertical, 41, 54, and 135 MHz, while at the rest of the stations, including the Cephalonia station, only the 3 kHz north-south, 3 kHz east-west, 10 kHz north-south, 10 kHz east-west, 41 and 46 MHz are recorded.More details on the instrumentation of the telemetric network can be found in the supplementary downloadable material of Potirakis et al. (2015).(For interpretation of the references to colours, the reader is referred to the online version of this paper.)

Critical dynamics analysis methods
Criticality has earlier been suggested as an EQ precursory sign (Chelidze, 1982;Chelidze and Kolesnikov, 1982;Chelidze et al., 2006;Rundle et al., 2012;Wanliss et al., 2015).Critical phenomena have been proposed as the likely framework to study the origins of EQ-related EM fluctuations, suggesting that the theory of phase transitions and critical phenomena may be useful in gaining insight into the mechanism of their complex dynamics (Bowman et al., 1998;Contoyiannis et al., 2004aContoyiannis et al., , 2005Contoyiannis et al., , 2015;;Varotsos et al., 2011a, b).One possible reason for the appropriateness of this model may be the way in which correlations spread thought a disordered medium/system comprised of subunits.From a qualitative/intuitive perspective, according to the specific approach, initially single isolated activated parts emerge in the system which, then, progressively grow and proliferate, leading to cooperative effects.Local interactions evolve to longrange correlations, eventually extending along the entire system.A key point in the study of dynamical systems that develop critical phenomena is the identification of the "critical epoch" during which the "short-range" correlations evolve into "long-range" ones.Therefore, the theory of phase transitions and critical phenomena seems to be appropriate for the study of dynamical complex systems in which local interactions evolve to long-range correlations, such as the disordered Earth's crust during the preparation of an EQ.It is worth noting that key characteristics of a critical point in a phase transition of the second order are the existence of highly correlated fluctuations and scale invariance in the statistical properties.By means of experiments on systems presenting this kind of criticality as well as by appropriately designed numerical experiments, it has been confirmed that right at the "critical point" the subunits are highly correlated even at an arbitrarily large "distance".At the critical state self-similar structures appear both in time and space.This fact is quantitatively manifested by power-law expressions describing the distributions of spatial or temporal quantities associated with the aforementioned self-similar structures (Stanley, 1987(Stanley, , 1999)).
The time-series analysis methods employed in this paper for the evaluation of the MHz EME recordings and the seismicity around Cephalonia in terms of critical dynamics are briefly presented in the following.Specifically, the method of critical fluctuations (MCF) is described in Sect.2.1, while the NT method is described in Sect.2.2.
It has been shown (Contoyiannis and Diakonos, 2000) that the dynamics of the order parameter fluctuations φ at the critical state for a second-order phase transition can be theoretically formulated by the non-linear intermittent map: where ϕ n is the scaled-order parameter value at the time interval n; u denotes an effective positive coupling parameter describing the non-linear self-interaction of the order parameter; and z stands for a characteristic exponent associated with the isothermal exponent δ for critical systems at thermal equilibrium (z = δ + 1).The marginal fixed point of the above map is the zero point, as expected from critical phenomena theory.However, it has been shown that in order to quantitatively study a real (or numerical) dynamical system, one has to add an unavoidable "noise" term, ε n , to Eq. ( 1), which is produced by all stochastic processes (Contoyiannis and Diakonos, 2007).Note that, from the intermittency mathematical framework point of view, the "noise" term denotes ergodicity in the available phase space.In this respect, the map of Eq. (1), for positive values of the order parameter, becomes Based on the map of Eq. (2), MCF has been introduced as a method capable of identifying whether a system is in critical state of intermittent type by analysing time series corresponding to an observable of the specific system.In a few words, MCF is based on the property of maps of intermittent type, like the ones of Eqs. ( 1) and ( 2), that the distribution of properly defined laminar lengths (waiting times) l follows a power-law P (l) ∼ l −p l (Schuster, 1998), where the exponent p l is p l = 1+ 1 δ (Contoyiannis et al., 2002).However, the distribution of waiting times for a real data time series which is not characterized by critical dynamics follows an exponential decay, rather than a power-law one (Contoyiannis et al., 2004a), due to stochastic noise and finite size effects.Therefore, the dynamics of a real time series can be estimated by fitting the distribution of waiting times (laminar lengths) to a function ρ(l) combining both power-law and exponential decay (Contoyiannis and Diakonos, 2007): The values of the two exponents p 2 and p 3 , which result after fitting laminar length distributions in a log-log scale diagram, reveal the underlying dynamics.Exact critical state calls for p 3 = 0; in such a case it is p 2 = p l > 1.As a result, in order for a real system to be considered to be at critical state, both criticality conditions p 2 > 1 and p 3 ≈ 0 have to be satisfied.Note that the choice of the function ρ(l) of Eq. ( 3), which combines both power-law and exponential decay, to model the distribution of waiting times was deliberately made in order to include both these fundamentally different behaviours, i.e. the critical dynamics (Contoyiannis et al., 2002) and the complete absence of specific dynamics (stochastic processes) (Contoyiannis et al., 2004b), respectively.In addition, the specific function also models intermediate behaviours (Contoyiannis and Diakonos, 2007).
In applying the MCF, the corresponding factors of ρ(l) appear to be competitive: any increase in the p 2 exponent value corresponds to a p 3 exponent value reduction and vice versa.However, this is expected, because, for example, any increase in the value of the p 3 exponent signifies the departure from critical dynamics and thus the reduction of a p 2 exponent value.Still, it is interesting to apply MCF analysis to observe this competition in the case of pre-earthquake EME time series and to explore whether the obtained exponent values are consistent with those of MCF analyses performed on other time series with large statistics which are considered references for the application of our method.This competition can be observed even within the critical windows as shown in Figs.2d and 3d.
Moreover, a special dynamics case is the one known as "tricritical crossover dynamics".In statistical physics, a tricritical point is a point in the phase diagram of a system at which the two basic kinds of phase transition, that is the first-order transition and the second-order transition, meet (Huang, 1987).The co-existence of three phases, specifi-cally, the high symmetry phase, the low symmetry phase, and an intermediate "mixing state", is a characteristic property of the area around this point.A passage through this area, around the tricritical point, from the second-order phase transition to the first-order phase transition through the intermediate mixing state, constitutes a tricritical crossover (Huang, 1987).
The specific dynamics is proved to be expressed by the map (Contoyiannis et al., 2015) where m stands for the order parameter.This map differs from the critical map of Eq. ( 2) in the sign of the parameter u and exponent z.Note that for reasons of unified formulation we use for these parameters the same notation as in the critical map of Eq. ( 2).At the level of MCF analysis this dynamics is expressed by the estimated values for the two characteristic exponent p 2 , p 3 values that satisfy the tricriticality condition p 2 < 1, p 3 ≈ 0. These values have been characterized in Contoyiannis and Diakonos (2007) as a signature of tricritical behaviour.Note that MCF analysis is possible only for time series which at least present cumulative stationarity.Therefore, a cumulative stationarity test is always performed before applying the MCF method; examples can be found in a number of publications (e.g.Contoyiannis et al., 2004aContoyiannis et al., , 2005Contoyiannis et al., , 2010;;Contoyiannis and Eftaxias, 2008;Potirakis et al., 2015).More details on the application of MCF can be found in several published articles (e.g.Contoyiannis et al., 2002Contoyiannis et al., , 2013Contoyiannis et al., , 2015)), as well as in Sect.3, where application to the MHz EM variations is presented.

Natural time (NT) method
The natural time method was originally proposed for the analysis for a point process like DC (direct current) or ultralow frequency (≤1 Hz) SES (seismic electric signals) (Varotsos et al., 2002;Varotsos, 2005), and has been shown to be optimal for enhancing the signals in the time-frequency space (Abe et al., 2005).The transformation of a time series of "events" from the conventional time domain to the natural time domain is performed by ignoring the time stamp of each event and retaining only their normalized order (index) of occurrence.Explicitly, in a time series of N successive events, the natural time, χ k , of the kth event is the index of occurrence of this event normalized by dividing by the total number of considered events, χ k = k/N .On the other hand, the "energy", Q k , of each, kth, event is preserved.We note that the quantity Q k represents different physical quantities for various time series: for EQ time series it has been assigned to a seismic energy released (e.g.seismic moment) (Varotsos et al., 2005;Uyeda et al., 2009b), and for SES signals that are of a dichotomous nature it corresponds to SES pulse duration (Varotsos, 2005), while for MHz electromagnetic emission signals that are of a non-dichotomous nature, it has been attributed to the energy of fracto-electromagnetic emission events as defined in Potirakis et al. (2013).The transformed time series (χ k ,Q k ) is then studied through the normalized power spectrum ( ) where is the natural angular frequency, = 2π ϕ, with ϕ standing for the frequency in natural time, termed "natural frequency", and p k = Q k / N n=1 Q n corresponds to the kth event's normalized energy.Note that the term "natural frequency" should not be confused with the rate at which a system oscillates when it is not driven by an external force; it defines an analysis domain dual to the natural time domain in the framework of a Fourier-Stieltjes transform (Varotsos et al., 2011b).
The study of ( ) at close to zero reveals the dynamic evolution of the time series under analysis.This is because all the moments of the distribution of p k can be estimated from ( ) at → 0 (Varotsos et al., 2011a).Aiming at that, by the Taylor expansion , i.e. the variance of χ k weighted for p k characterizing the dispersion of the most significant events within the "rescaled" interval (0, 1].Moreover, the entropy in natural time, S nt , is defined (Varotsos et al., 2006) as and corresponds (Varotsos et al., 2006(Varotsos et al., , 2011b) ) to the value at q = 1 of the derivative of the fluctuation function F (q) = χ q − χ q with respect to q (while κ 1 corresponds to F (q) for q = 2).It is a dynamic entropy depending on the sequential order of events (Varotsos et al., 2006).The entropy, S nt-, obtained upon considering (Varotsos et al., 2006) the time reversal T , i.e.T p m = p N−m+1 , is also considered.
A system is considered to approach criticality when the parameter κ 1 converges to the value κ 1 = 0.070 and at the same time both the entropy in natural time and the entropy under time reversal satisfy the condition S nt , S nt-< S u = (ln 2/2) − 1/4 (Sarlis et al., 2011), where S u stands for the entropy of a "uniform" distribution in natural time (Varotsos et al., 2006).
It should be finally clarified that in the case of seismicity analysis, the temporal evolution of the parameters κ 1 , S nt , S nt-, and D is studied as new events that exceed the magnitude threshold M thres are progressively included in the analysis.Specifically, as soon as one more event is included, first the time series (χ k , Q k ) is rescaled in the natural time domain, since each time the kth event corresponds to a natural time χ k = k/N , where N is the progressively increasing (by each new event inclusion) total number of the considered successive events; then all the parameters involved in the natural time analysis are calculated for this new time series; this process continues until the time of occurrence of the main event.
More details on the application of NT to MHz EME as well as to foreshock seismicity can be found in Potirakis et al. (2013Potirakis et al. ( , 2015)), as well as in Sects.3 and 4, where its application to the MHz EM variations and foreshock seismicity is presented, respectively.
Note that in the case of NT analysis of foreshock seismicity, the introduction of a magnitude threshold, M thres , excludes some of the weaker EQ events (with magnitude below this threshold) from the NT analysis.On the one hand, this is necessary in order to exclude events for which the recorded magnitude is not considered reliable; depending on the installed seismographic network characteristics, a specific magnitude threshold is usually defined to ensure data completeness.On the other hand, the use of various magnitude thresholds, M thres , offers a means of more accurate determination of the time when criticality is reached.In some cases, it happens that more than one time point may satisfy the rest of NT critical state conditions; however, the time of the true coincidence is finally selected by the last condition that "true coincidence should not vary upon changing (within reasonable limits) either the magnitude threshold, M thres , or the area, used in the calculation".

Electromagnetic emissions analysis results
Some of the MHz recordings of the Cephalonia station associated with the M w = 6.0 EQ (EQ1) are shown in Fig. 2a.This was recorded on day of year 24, that is ∼ 2 days before the occurrence of EQ1.This stationary time-series excerpt, having a total length of 2.8 h (10 000 samples) starting at 24 January 2014 (12:46:40 UT), was analysed by the MCF method and was identified as a "critical window" (CW).CWs are time intervals of the MHz EME signals presenting features analogous to the critical point of a second-order phase transition (Contoyiannis et al., 2005).
The main steps of the MCF analysis (e.g.Contoyiannis et al., 2013Contoyiannis et al., , 2015) ) of the specific time series are shown in Fig. 2b-d.First, a distribution of the amplitude values of the analysed signal was obtained from which, using the method of turning points (Pingel et al., 1999), a fixed point, that is the start of laminar regions, φ o of about 700 mV was determined.Figure 2c portrays the obtained distribution of laminar lengths for the end point φ l = 655 mV, that is the distribution of waiting times, referred to as laminar lengths l, between the fixed point φ o and the end point φ l , as well as the fitted function f (l) ∝ l −p 2 e −p 3 l with the corresponding exponents p 2 = 1.35, p 3 = 0.000 with R 2 = 0.999.Note that the distribution of laminar lengths is directly fitted to the specific model using the Levenberg-Marquardt algorithm, while the fitting criterion is the chi-square minimization.The fitting is not done in log-log space; the axes of Fig. 2c are logarithmic for the easier depiction of the distribution of laminar lengths.Finally, Fig. 2d shows the obtained plot of the p 2 ,p 3 exponents vs. φ l .From Fig. 2d it is apparent that the criticality conditions, p 2 > 1 and p 3 ≈ 0, are satisfied for a wide range of end points φ l , revealing the power-law decay feature of the time series that proves that the system is characterized by intermittent dynamics; in other words, the MHz time-series excerpt of Fig. 2a is indeed a CW.Some of the MHz recordings of the Zante station associated with EQ1 are shown in Fig. 3a.This was also recorded on day of year 24, that is ∼ 2 days before the occurrence of Cephalonia EQ1.This stationary time-series excerpt, having a total length of 2.8 h (10 000 samples) starting at 24 January 2014 (12:46:40 UT), was also analysed by the MCF method and was identified as a CW.
The application of the MCF analysis to the specific time series (cf.Fig. 3) revealed that the criticality conditions, p 2 > 1 and p 3 ≈ 0, are satisfied for a wide range of end points φ l , for this signal too.In other words, this signal has also embedded the power-law decay feature that indicates intermittent dynamics, rendering it a CW.
After the M w = 6.0 (EQ1), ∼ 1 week later, the second, M w = 5.9 (EQ2), occurred on the same island with a focal area a few kilometres away from the first one.Six days earlier, both the Cephalonia and Zante stations simultaneously recorded MHz EME.Specifically, a stationary time-series excerpt, having a total length of 3.3 h (12 000 samples) starting at 28 January 2014 (05:33:20 UT), from the Cephalonia station, and a stationary time-series excerpt, having a total length of 5 h (18 000 samples) starting at 28 January 2014 (03:53:20 UT), from the Zante station, were analysed by the MCF method, and both of them were identified as CWs.Note that the Cephalonia CW was emitted within the time frame in which the Zante CW was emitted.Figures 4 and 5 show the results of the corresponding analyses.
In summary, we conclude that, according to the MCF analysis method, both stations recorded MHz signals that simultaneously presented critical state features 2 days before the first main event and 6 days before the second main event.
In order to verify this finding, we proceeded to the analysis of all the corresponding MHz signals by means of the NT analysis method, according to the way of application pro- .Temporal evolutions of the four natural time (NT) analysis parameters (κ 1 , S nt , S nt-, and D ) for the foreshock seismic activity recorded prior to EQ1: (a) for the activity of the whole investigated area of the Ionian Sea for M L threshold 2.5, during the period from 13 December 2013 00:00:00 to 26 January 2014 13:55:44 UT (just after the occurrence of EQ1); (b) for the activity of the whole investigated area of the Ionian Sea for M L threshold 2.3, during the period from 11 January 2014 04:13:00 (just after the M L = 4.7 that occurred in Zante) to 26 January 2014 13:55:44 UT; (c) for the activity of both Cephalonia (eastern and western) zones combined for M L threshold 2.1, during the period from 13 December 2013 00:00:00 to 26 January 2014 13:55:44 UT; and (d) for the activity of western Cephalonia for M L threshold 2.1, during the period from 13 December 2013 00:00:00 to 26 January 2014 13:55:44 UT.Note that the events employed depend on the considered threshold.Moreover, the time (x) axis is not linear in terms of the conventional date of occurrence of the events, since the employed events appear equally spaced relative to the x axis, as the natural time representation demands, although they are not equally spaced in conventional time.The horizontal solid light green, solid grey, the grey dashed lines denote the same quantities as in Fig. 6, while the horizontal solid brown line denotes the 10 −2 limit for D .(For interpretation of the references to colours, the reader is referred to the online version of this paper.)posed in Potirakis et al. (2013).According to the specific procedure for the application of the NT method to the MHz signals, we performed an exhaustive search, seeking at least one amplitude threshold value (applied over the total length of the analysed signal) for which the corresponding fracto-EME events satisfy the natural time method criticality conditions.The idea is that if the MCF gives valid information, and as a consequence the analysed time-series excerpt is indeed in critical condition, then there should be at least one threshold value for which the NT criticality conditions (cf.Sect.2.2) are satisfied.Indeed, as apparent from Fig. 6, all four signals satisfy the criticality conditions according to the NT method for at least one of the considered threshold values; therefore, the results obtained by the MCF method are successfully verified.
On 2 February 2014, i.e. 1 day before the occurrence of EQ2, MHz EME presenting tricritical characteristics was recorded by the Cephalonia station.This signal emerged 5 days after the CWs that were identified in the MHz EME si-multaneously recorded by the Cephalonia and Zante stations.The specific MHz time-series excerpt from the Cephalonia station, having a total length of 7.5 h (27 000 samples) starting at 2 February 2014 (07:46:40 UT), was analysed by means of the MCF method, yielding the results shown in Fig. 7.As is apparent from the results, this signal satisfies the tricriticality conditions p 2 < 1, p 3 ≈ 0 (cf.Sect.2.1) for a wide range of end points φ l , revealing the intermediate "mixing state" between the second-order phase transition and the first-order phase transition.Unfortunately, during the time that the Cephalonia station recorded a tricritical MHz signal, the Zante station was not in operation; actually, it was out of operation for several hours during that specific day.
It has been recently found (Contoyiannis et al., 2015) that such a behaviour is identified in the kHz EMEs which usually emerge near the end of the MHz EMEs when the environment in which the EQ preparation process evolves changes from heterogeneous to less heterogeneous, and before the emergence of the strong avalanche-like kHz EMEs which  have been attributed to the fracture of the asperities sustaining the fault.Actually, this has been proposed as the second stage of the four-stage model for the preparation of an EQ by means of its observable EM activity (Eftaxias and Potirakis, 2013;Contoyiannis et al., 2015, and references therein;Donner et al., 2015).The identification of tricritical behaviour in MHz EMEs is a quite important finding, indicating that the tricritical behaviour, attributed to the second stage of the aforementioned four-stage model, can be identified either in kHz or in MHz EMEs, leading thus to a revision of the specific four-stage model in order to include this case too.As a conclusion, after the first stage of the EQ preparation process where MHz EME with critical features are emitted, a second stage follows where MHz or kHz or both MHz and kHz EME with tricritical features are emitted.As already mentioned (cf.Sect.2.1), in terms of statistical physics the tricritical behaviour is an intermediate dynamical state which is developed in the region of the phase diagram of a system around the tricritical point, which can be approached either from the edge of the first-order phase transition (characterizing the strong avalanche-like kHz EME attributed to the third stage of the four-stage model) or from the edge of the secondorder phase transition (characterizing the critical MHz EME attributed to the first stage of the four-stage model).Therefore, although it is expected that the tricritical behaviour will be rarely observed, as has already been discussed in Contoyiannis et al. (2015), it can be found either in MHz time series, following the emission of a critical MHz EME, or in kHz time series preceding the emission of avalanche-like kHz EME.

Foreshock seismic activity analysis results
As already mentioned in Potirakis et al. (2013Potirakis et al. ( , 2015)): "seismicity and pre-fracture EMEs should be two sides of the same coin concerning the EQ generation process.If the MHz EMEs and the corresponding foreshock seismic sequence are observable manifestations of the same complex system at critical state, both should be possible to be described as a critical phenomenon by means of the natural time method".Therefore, we also proceeded to the examination of the corresponding foreshock seismic activity around Cephalonia before each one of the significant EQs of interest in order to verify this suggestion.However, we did not apply the NT method to concentric circles around the epicentre of each EQ, as in Potirakis et al. (2013Potirakis et al. ( , 2015)), but instead we decided to study seismicity within areas determined according to seismotectonic and earthquake hazard criteria.
Following the detailed study presented in Vamvakaris et al. (2016), we incorporated the seismic zones proposed there for our area of study.Thus, as presented in Fig. 8, we defined five separate seismic zones, based on the criteria explored in Vamvakaris et al. (2016) and the seismic zonation proposed by them.Since the study area comprises the most seismically active zone in Greece, assigned also the highest value in the Earthquake Building Code for the country, a large number of source, stress, and strain studies have been used in their study to establish such a definition of zoning.Hence, it was found to be well justified to follow their zone definition.In Fig. 8, from east to west and north to south, one can identify the zones of Akarnania (area no.1), Lefkada (area no. 2), eastern Cephalonia (area no.3), western Cephalonia (area no.4), NT analysis Between 1st and 2nd EQ (after 29/01) Figure 12.Natural time (NT) analysis results for the seismicity in the partition of western Cephalonia during the time period from 29 January 2014 00:00:00 to 03 February 2014 03:08:47 UT (between EQ1, M w = 6.0, and EQ2, M w = 5.9): (a-d) temporal evolutions of the four natural time analysis parameters (κ 1 , S nt , S nt-, and D ) for the different M L thresholds 2.2, 2.6, 2.8, and 3.0, respectively.Note that the events employed depend on the considered threshold.Moreover, the time (x) axis is not linear in terms of the conventional date of occurrence of the events, since the employed events appear equally spaced relative to the x axis, as the natural time representation demands, although they are not equally spaced in conventional time.The horizontal solid light green, solid grey, and grey dashed lines denote the same quantities as in Fig. 6, while the horizontal solid brown line denotes the 10 −2 limit for D .(For interpretation of the references to colours, the reader is referred to the online version of this paper.) and Zante (area no.5), respectively, covering the area of the Ionian Sea near Cephalonia.
Before we proceed to the NT analysis of seismicity, the seismic activity prior to EQ1 as well as between EQ1 and EQ2 is briefly discussed in relation to the above-mentioned seismic zones.Earthquake parametric data have been retrieved from the National Observatory of Athens online catalogue (http://www.gein.noa.gr/en/seismicity/earthquake-catalogs), while for all the presented maps and calculations the local magnitude (M L ), as provided by the specific earthquake catalogue, is used.The foreshock seismic activity before EQ1 for the whole investigated area of the Ionian Sea region from 13 December 2013 up to the time of occurrence of the main event is shown in the map of Fig. 9a.As can be easily observed from this map, there was a high seismic activity mainly focused on two specific zones: western Cephalonia and Zante.Notably, an EQ of M L = 4.7 occurred in Zante on 11 January 2014 04:12:58, indicated by the black arrow in Fig. 9a.No EQs were recorded in Akarnania, while very few events were recorded in Lefkada and eastern Cephalonia.The events which occurred in western Cephalonia are also shown in a separate map in Fig. 9b for later reference.
Applying the natural time analysis to seismic data (cf.Sect.2.2), the evolution of the time series (χ k , Q k ) was studied for the foreshock seismicity prior to EQ1, where Q k is in this case the seismic energy released during the kth event.The seismic moment, M 0 , as proportional to the seismic energy, is usually considered (Varotsos et al., 2005;Uyeda et al., 2009b;Potirakis et al., 2013Potirakis et al., , 2015)).Our calculations were based on the seismic moment M 0 (in dyn cm) resulting from the corresponding M L as (Varotsos et al., 2005;Potirakis et al., 2013Potirakis et al., , 2015) ) M 0 = 10 0.99M L +11.8 .First, we performed an NT analysis on the seismicity activity of the whole investigated Ionian Sea region during the period from 13 December 2013 00:00:00 to 26 January 2014 13:55:44 UT, i.e. just after the occurrence of EQ1, for different magnitude thresholds, M thres , for which all earthquakes with M L > M thres were included in the analysis.Note that only M thres ≥ 2 was considered in order to ensure data completeness (Chouliaras et al., 2013a, b).
For all the considered threshold values, the result was the same: no indication of criticality was identified (see for example Fig. 10a).Since, as we have already mentioned, the whole investigated area was mainly dominated by the seismic activity in western Cephalonia and the seismic activ-ity in Zante, while an EQ of M L = 4.7 occurred in Zante, we decided to start the NT analysis after the occurrence of the specific Zante EQ, in order to exclude from our analysis possible foreshock activity related to the specific event.As a result, we performed NT analysis for the time period 11 January 2014 04:13:00 (just after the M L = 4.7 Zante EQ) to 26 January 2014 13:55:44 UT, for different magnitude thresholds in three successively enclosed areas, namely, the whole investigated area of the Ionian Islands region, both Cephalonia (eastern and western) zones combined, and the zone of western Cephalonia.Representative examples of these analyses are depicted in Fig. 10b-d.The analysis over the whole investigated area of the Ionian Islands region indicates that seismicity reaches criticality on 19 and 20 January, while the other two progressively narrower areas indicate that the criticality conditions according to the NT method are satisfied on 19 and 22 January.These results imply that seismicity was also in critical condition a few days prior to the occurrence of the first studied significant Cephalonia EQ (EQ1).Actually, in the specific case, the critical condition of seismicity was reached before, but quite close to the emission of the corresponding MHz signals for which critical behaviour was identified (cf.Sect.3).Note that a very recent analysis of the foreshock seismic activity before EQ1, in terms of a combination of multi-resolution wavelets and NT analysis, which was performed in concentric areas of 50 and 30 km radii around the epicentre of EQ1, also found that NT analysis criticality requirements are met a few days before EQ1 (at approximately 20 January) (Vallianatos et al., 2015).
Before the application of the NT method to the seismic activity prior to EQ2, one should first study the time evolution of the activity between the two significant events of interest, in order to minimize if possible the influence of the first EQ aftershock sequence on the NT analysis.Our first observation about the EQs which occurred during the specific time period was that all but one had epicentres in western Cephalonia.Only one M L = 2.3 EQ occurred in Zante, at (37.79 • N, 21.00 • E) on 28 January 2014 02:08:27 UT.
Figure 11a shows all the events that were recorded in the whole investigated area of the Ionian Islands region vs. time from just after EQ1 (M w = 6.0) up to the time of EQ2 (M w = 5.9), including EQ2.As can be seen, if one considers that both significant EQs of interest were main events, it is quite difficult to separate the seismic activity of the specific time period into aftershocks of the first EQ and foreshocks of the second one.However, we observe that up to a specific time point, there is a rapid decrease in the running mean magnitude of the recorded EQs, while after that the long-range (75 events) running mean value seems to be almost constant over time, with the short-range (25 events) one varying around it.We arbitrarily set 29 January 00:00:00 UT as the time point after which the recorded seismicity is no longer dominated by the aftershocks of EQ1; this by no means implies that the aftershock sequence of the EQ1 stops after that date.It should also be underlined that changing this, arbi-trarily selected, date within reasonable limits, does not significantly change the results of our corresponding NT analysis which are presented next.On the other hand, a significant shift of this limit towards EQ1, i.e. to earlier dates, results in severe changes, indicating the domination of the recorded seismicity by the aftershock sequence of EQ1.Accordingly, the considered foreshock seismic activity before EQ2, i.e. from 29 January 2014 00:00 UT up to the time of occurrence of EQ2, is presented in the map of Fig. 11b for western Cephalonia and is analysed in the following.
Next, we applied the NT method to the seismicity of western Cephalonia for the time period from 29 January 2014 00:00:00 to 3 February 2014 03:08:47 UT.Note that we also applied the NT method to the whole investigated area of the Ionian Islands region, obtaining practically the same results.As we have already mentioned, only one M L = 2.3 EQ occurred outside the western Cephalonia zone, so on the one hand for magnitude threshold values M thres ≥ 2.3 this event was excluded, while, on the other hand, even for lower threshold values (2 ≤ M thres < 2.3), its inclusion does not change the results significantly.Figure 12 shows the NT analysis results for some threshold values, proving that seismicity reaches criticality on 1 or 2 February 2014, that is 1 or 2 days before the occurrence of the second significant EQ of interest (M w = 5.9).Actually, in the specific case, the critical condition of seismicity was reached after, but quite close to, the emission of the corresponding MHz signals for which critical behaviour was identified (cf.Sect.3).

Discussion and conclusions
Based on the methods of critical fluctuations and natural time, we have shown that the fracture-induced MHz EMEs recorded by two stations in our network prior to two recent significant EQs that occurred in Cephalonia present criticality characteristics, implying that they emerge from a system in critical state.
There are two key points that render these observations unique in the research up to now on the pre-EQ EME.
(i) Although not previously reported, we have to note here that, as has been observed, the Cephalonia station is insensitive to EQ preparation processes happening outside of the wider area of Cephalonia, as well as to EQ preparation processes leading to low magnitude EQs within the area of Cephalonia.Note that the only signal that has been previously recorded refers to the M = 6 EQ that occurred on the specific island in 2007 (Contoyiannis et al., 2010).
(ii) Prior to each one of the studied significant EQs, two MHz EME time series presenting critical characteristics were recorded simultaneously in two different stations very close to the focal areas, while no other station of our network (cf.Fig. 1) recorded such signals prior to the specific EQs.This indicates that the revealed criticality was not associated with a global phenomenon, such as critical variations in the ionosphere, but was rather local to the area of the Ionian Islands region, enhancing the hypothesis that these EME were associated with the EQ preparation process taking place prior to the two significant EQs.This feature, combined with the above-mentioned sensitivity of the Cephalonia station only to significant EQs occurring on the specific island, could have been considered to be an indication of the location of the impending EQs.
Electromagnetic emissions, as a phenomenon rooted in the damage process, should be an indicator of memory effects.Laboratory studies verify that during cyclic loading, the level of EME increases significantly when the stress exceeds the maximum previously reached stress level (Kaiser effect).The existence of the Kaiser effect predicts the EM silence during the aftershock period (Eftaxias et al., 2013;Eftaxias and Potirakis, 2013, and references therein).Thus, the appearance of the second EM anomaly may reveal that the corresponding preparation of a fracture process has been organized in a new barrier.
We note that, according to the view that seismicity and pre-EQ EM emissions should be "two sides of the same coin" concerning the earthquake generation process, the corresponding foreshock seismic activity, as another manifestation of the same complex system, should be at critical state as well, before the occurrence of a main event.We have shown that this really happens for both the significant EQs we studied.Importantly, the revealed critical process seems to be focused on an area corresponding to the western Cephalonia zone, one of the parts according to the seismotectonic and hazard zone partitioning of the wider area of the Ionian Islands.
To be more detailed, the foreshock seismicity associated with the first (M w = 6.0)EQ reached critical condition a few days before the occurrence of the main event.Specifically, it came to critical condition before, but quite close to, the emission of the corresponding MHz signals for which critical behaviour was identified.The seismicity that was considered to be a foreshock of the second (M w = 5.9) EQ also reached criticality a few days before the occurrence of the main event.
In contrast to the first EQ case, it reached criticality after, but quite close to, the emission of the corresponding MHz signals for which critical behaviour was identified.
One more outcome of our study was the identification of tricritical crossover dynamics in the MHz emissions recorded just before the occurrence of the second significant EQ of interest (M w = 5.9) at the Cephalonia station.Note that, unfortunately, the Zante station was out of order for several hours during the specific day, including the time window during which the tricritical features were identified in the Cephalonia recordings.As a result, we could not cross check whether tricritical signals simultaneously also reached Zante.This is considered a quite important finding, since it verifies a theoretically expected situation, namely the approach of the intermediate dynamical state of tricritical crossover, either from the first-or second-order phase transition state.In terms of pre-EQ EME, this leads to a revision of the four-stage model for the preparation of an EQ by means of its observable EM activity.That is, after the first stage of the EQ preparation process where MHz EMEs with critical features are emitted, a second stage follows where MHz or kHz or both MHz and kHz EMEs with tricritical features are emitted.Specifically, the tricritical crossover dynamics can be identified either in MHz time series, following the emission of a critical MHz EME, or in kHz time series preceding the emission of avalanche-like kHz EME.In summary, the proposed four stages of the last part of EQ preparation processes and the associated, appropriately identified, EM observables appear in the following order: first stage: valid MHz anomaly; second stage: MHz or kHz or MHz and kHz anomaly exhibiting tricritical characteristics; third stage: strong avalanche-like kHz anomaly; fourth stage: electromagnetic quiescence.Note that the specific four-stage model is a suggestion that seems to be verified by the MHz-kHz observation data available up to now and corresponding time-series analyses, while a rebuttal has not appeared in the literature.However, the understanding of the physical processes involved in the preparation of an EQ and their relation to various available observables is an open scientific issue.A lot of efforts still remain to be undertaken before one can claim clear understanding of EQ preparation processes and their associated possible precursors.
As has been repeatedly pointed out in previous works (e.g.Eftaxias et al., 2013;Eftaxias and Potirakis, 2013, and references therein), our view is that such observations and the associated analyses offer valuable information towards the understanding of processes in the Earth system taking place prior to the occurrence of a significant EQ.As is known, a large number of other precursory phenomena are also observed, both by ground and satellite stations, prior to significant EQs.Only a combined evaluation of our observations with other well-documented precursory phenomena could possibly render our observations useful for a reliable short-term forecast solution.Unfortunately, in the cases of the Cephalonia EQs under study, this requirement was not fulfilled.To the best of our knowledge, only one paper reporting the emergence of VLF seismic-ionospheric disturbances 4 days before the first Cephalonia EQ (Skeberis et al., 2015) has been published up to now.It is very important that the specific disturbances, which also correspond to a spatially extensive process as happens with the MHz EME, were recorded during the same time window with the herein presented MHz critical signals.However, more precursory phenomena could have been investigated if appropriate observation data had been available.For example, if groundbased magnetic observatories in the area of Greece had available magnetometer data for the time period of interest, EQwww.nonlin-processes-geophys.net/23/223/2016/ Nonlin.Processes Geophys., 23, 223-240, 2016 S. M. Potirakis et al.: Recent seismic activity at Cephalonia related ULF magnetic field variations could also be investigated.Such ULF variations, either of lithospheric or ionospheric origin, are also a result of spatially extensive processes and have been shown to present critical characteristics prior to significant EQ occurrence (Hayakawa et al., 2015a, b;Contoyiannis et al., 2016;Potirakis et al., 2016).

Data availability
The used earthquake parametric data can be retrieved from the National Observatory of Athens online catalogue (http:// www.gein.noa.gr/en/seismicity/earthquake-catalogs), while the electromagnetic recordings data can be obtained by email request to S. M. Potirakis (spoti@teipir.gr).
Figure1.Map with the distribution of stations of the telemetric network that monitors electromagnetic variations in the MHz and kHz bands in Greece, which were operating during the time period of interest.The locations of the Cephalonia and Zante stations are marked by the magenta square and triangle, respectively, while the rest of the remote stations are denoted by red circles and the central data recording server by a blue circle.The epicentres of the two significant EQs of interest are also marked, the first (EQ1, M w = 6.0) by a red cross and the second (EQ2, M w = 5.9) by a green X mark.The electromagnetic variations monitored at the Zante station are 3 kHz north-south, 3 kHz east-west, 3 kHz vertical, 10 kHz north-south, 10 kHz east-west, 10 kHz vertical, 41, 54, and 135 MHz, while at the rest of the stations, including the Cephalonia station, only the 3 kHz north-south, 3 kHz east-west, 10 kHz north-south, 10 kHz east-west, 41 and 46 MHz are recorded.More details on the instrumentation of the telemetric network can be found in the supplementary downloadable material ofPotirakis et al. (2015).(For interpretation of the references to colours, the reader is referred to the online version of this paper.)

Figure 2 .
Figure 2. (a) The 10 000-sample long critical window of the MHz EME that was recorded before the Cephalonia M w = 6.0 EQ at the Cephalonia station.(b) Amplitude distribution of the signal of Fig. 2a.(c) Distribution of laminar lengths for the end point φ l = 655 mV, as a representative example of the involved fitting.The solid line corresponds to the fitted function (cf.text in Sect.2.1) with the values of the corresponding exponents p 2 , p 3 also noted.(d) The obtained exponents p 2 , p 3 vs.different values of the end of laminar region φ l .The horizontal dashed line indicates the critical limit (p 2 = 1).

Figure 3 .
Figure 3. (a)The 10 000-sample long critical window of the MHz EME that was recorded prior to the Cephalonia M w = 6.0 EQ at the Zante station, while (b), (c), and (d) are similar to the corresponding parts of Fig.2.From (b), a fixed point (start of laminar regions) φ o of about 600 mV results, while in (c), the distribution of laminar lengths is given for the end point φ l = 665 mV for which the exponents p 2 = 1.49, p 3 = 0.000 with R 2 = 0.999 were obtained.

Figure 4 .
Figure 4. (a)The 12 000-sample long critical window of the MHz EME that was recorded before the Cephalonia M w = 5.9 EQ at the Cephalonia station, while (b), (c), and (d) are similar to the corresponding parts of Fig.2.In (c), the distribution of laminar lengths is given for the end point φ l = 660 mV.

Figure 5 .
Figure 5. (a)The 18 000-sample long critical window of the MHz EME that was recorded before the Cephalonia M w = 5.9 EQ at the Zante station; (b), (c), and (d) are similar to the corresponding parts of Fig.2.In (c), the distribution of laminar lengths corresponds to the end point φ l = 400 mV.

Figure 6 .
Figure 6.Natural time analysis results obtained for the MHz EME signals shown in (a) Fig. 2a, recorded at the Cephalonia station prior to EQ1, (b) Fig. 3a, recorded at the Zante station prior to EQ1, (c) Fig. 4a, recorded at the Cephalonia station prior to EQ2, and (d) Fig.5a, recorded at the Zante station prior to EQ2.The quantities κ 1 (solid curve), S nt (dash-dot curve), and S nt-(dot curve) vs. amplitude threshold for each MHz signal are shown.The entropy limit of S u (≈ 0.0966), the value 0.070, and a region of ±0.005 around it are denoted by the horizontal solid light green, solid grey, and grey dashed lines, respectively.(For interpretation of the references to colours, the reader is referred to the online version of this paper.)

Figure 7 .
Figure 7. (a) The 27 000-sample long tricritical excerpt of the MHz EME that was recorded before the Cephalonia M w = 5.9 EQ at the Cephalonia station; (b), (c), and (d) are similar to the corresponding parts of Fig. 2. In (c), the distribution of laminar lengths corresponds to the end point φ l = 675 mV.

Figure 8 .Figure 9 .
Figure 8. Seismic zonation in the Ionian Islands area.The locations of the Cephalonia and Zante stations, as well as the epicentres of the two significant EQs of interest, are marked using the same signs presented in Fig. 1.

Figure 11
Figure11.(a) Seismic activity from the time immediately after EQ1 (M w = 6.0) up to the time of EQ2 (M w = 5.9) for the whole investigated area of the Ionian Sea.The moving averages of the recorded earthquake local magnitudes vs. time for calculation windows of 25 and 75 successive events are shown by the dashed magenta and solid grey curve, respectively.The vertical solid red line denotes the time point 29 January 00:00:00 UT.(b) The seismic activity considered to be foreshocks before EQ2 (from 29 January 2014 00:00 UT up to the time of occurrence of EQ2) for western Cephalonia.All presented magnitudes are local magnitudes (M L ). (For interpretation of the references to colours, the reader is referred to the online version of this paper.)