Long-Lasting Patterns in 3 kHz Electromagnetic Time Series after the M L = 6.6 Earthquake of 2018-10-25 near Zakynthos, Greece

: This paper reports one-month 3 kHz EM disturbances recorded at Kardamas, Ilia, Greece after a strong M L = 6.6 earthquake occurred on 2018 / 10 / 25 near Zakynthos and Ilia. During this period 17 earthquakes occurred with magnitudes M L = 4.5 and M L = 5.5 and depths between 3 km and 17 km, all near Zakynthos and Ilia. A two-stage, fully computational methodology was applied to the outcomes of ﬁve di ﬀ erent time-evolving chaos analysis techniques (DFA, fractal dimension analysis through Higuchi, Katz and Sevcik methods and power-law analysis). Via literature-based thresholds, the out-of-threshold results of all chaos analysis methods were located and from these, the common time instances of 13 selected combinations per ﬁve, four, three and two methods. Numerous persistent segments were located with DFA exponents between 1.6 ≤ α ≤ 2.0, fractal dimensions between 1.4 ≤ D ≤ 2.0 and power-law exponents between 2.2 ≤ β ≤ 3.0. Out of the 17 earthquakes, six earthquakes were jointly matched by 13 selected combinations of ﬁve, four, three and two chaos analysis methods, four earthquakes by all combinations of four, three and two, while the remaining seven earthquakes were matched by at least one combination of three methods. All meta-analysis matches are within typical forecast periods. 17 Five different time-evolving


Introduction
Earthquakes are natural phenomena that negatively impact society. Strong earthquakes are a major concern because they are destructive and occur inevitably when certain geophysical conditions chaos analysis methods are employed (detrended fluctuation analysis-DFA, analysis of fractal dimensions via Higuchi's, Katz's and Sevcik's methods and power-law spectral fractal analysis). All these methods have been used by the reporting team with success in several pre-earthquake EM and radon signals in Greece [1][2][3][4]9,[13][14][15][16][17][18]20,21]. The goal is to discover if long-lasting and fractal trends exist in time-evolving sections of the recorded EM time series that could be considered to be signs of predictability prior to each one of the 17 earthquakes that constitute the post-seismic activity of the M L = 6.6 Zante earthquake, assessing these as potential precursors. Via a modern two-stage computational approach (meta-analysis methodology) [49], the time instances are located and stored in which all chaos analysis results are out of literature-based thresholds and, simultaneously, all five methods or any possible combination of four, three or two methods coincide. Through this methodology, joint, highly persistent, fractional Brownian motion (fBm) EM segments with significant predictability and enhanced precursory value are located and isolated from the low-predictability ones. Associations are investigated among certain seismic events of the post-seismic period based on time-sequences of matching combinations of chaos analysis results. The potential geological source models are discussed and analyzed.

Geology and Seismic Significance of the Area
Ilia is situated on a large depression structure (graben) on the outer part of the Hellenic Trench delimited by the convergence of the Apulian, African and Aegean Sea plates ( Figure 1). Due its position, Ilia is associated with significant active seismic structures and important earthquakes in Greece [50,51], a country very prone to earthquakes due to its position on the convergence between the Eurasian, African, Aegean Sea and Anatolian plates ( Figure 1). The position and geological setting of Ilia and Kardamas make the study site very significant for investigating tectonic activity and especially for collecting disturbances related to earthquakes [18,52]. Ilia gave more than 600 earthquakes with M L ≥ 4.0 six of which were very destructive [53,54]: The Kato Achaia earthquake (greater of the above), evoked significant radon activity to Kardamas [13,18], which was found to be of similar patterns with the EM disturbances of the MHz [4,16,17] and kHz ranges [4]. Table 1 presents significant data for the 17 earthquakes of this paper. The earthquakes of Table 1 occurred between 2018/10/26 and 2018/11/26 with magnitudes M L ≥ 4.5, depths ≤ 20 km and epicenter's locations within a circle of radius of 101 km from the Kardamas station. All these earthquakes occurred after the great M L = 6.6 earthquake of 2018/11/25 near Zakynthos Island, Greece and comprise, as aforementioned, the post-seismic activity of this earthquake. Figure 2 presents all epicenters of the earthquakes of Table 1 in a map centered at the Kardamas station. As can be observed, the post-seismic activity is composed only by shallow earthquakes with most events occurring with depths ≤ 11 km. This serendipitous finding is deemed as important by others, e.g., [10,11] for kHz electromagnetic pre-earthquake recordings. It is also acknowledged by the authors of this paper in a recent publication [4] during significant seismic activity of Lesvos Island, Greece on 2017. The reader should note that the events are gathered both spatially and temporally into two main groups between 2018/10/26-2018/11/05 (JD:299-309) and 2018/11/11-2018/11/19 (JD:315-324), i.e., with a five-day pause. This important observation is rare to encounter and is extensively discussed in a recent publication of the authors [9] where a similar earthquake occurrence pattern was found. Moreover, the use of kHz radiation signifies further the earthquakes included in this paper. Indeed, as already mentioned in Section 1, numerous papers published during the last 20 years, indicate that the kHz electromagnetic emissions are important precursors of earthquakes, see, e.g., review of [10]. The kHz electromagnetic emissions have been used by the team with success in a very destructive earthquake in Lesvos Island, Greece [4,9]. Please note that the 3 kHz antennas used in this paper are selected after significant search so as to be clear from any artificial emissions in Greece according to the allocated Hellenic frequency band [16]. Table 1. Earthquakes included in the present study in chronological sequence. Actual dates and corresponding Julian Days (JD) of 2009 are also given. Abbreviations: M L : Local Magnitude, Lt: Latitude, Lg: Longitude, Depth: Epicentral depth, Dist: Distance of the earthquake epicenter from the Kardamas station. All data, except symbol and JD, are retrieved by [53]. As can be deduced from the above argumentation, the analyzed earthquakes of this paper are significant for studying the post-seismic activity of the M L = 6.6 earthquake which occurred near Zakynthos Island, Greece on 2018/11/25. In the following sections, several arguments will be presented that support the aspect of considering 3 kHz electromagnetic disturbances as seismic precursors of earthquakes, especially when combinations of different methods are employed that can detect hidden fractal and long-memory trends in time series.

Instrumentation
The electromagnetic disturbances of the 3 kHz range are continuously monitored in the Kardamas station by (i) circular magnetic field antennas synchronized properly at 3 kHz; (ii) Cambel CR-10 data-logger with 2-h buffer; (iii) telemetry equipment sending continuously the measurements to a personal computer at the rate of 1 Hz.
It should be noted that the Kardamas station is situated along the Hellenic Trench and the Outer non-volcanic Arc. According to several publications [13,14,[16][17][18] the measurement site is very sensitive for collecting disturbances related to earthquakes.

Fractal and Long Memory
There is a variety of physical systems in nature which can be described with fractals. The fractal behavior is characteristically found when the system or a part of it, is dilated, translated, or rotated in space. Depending on the mathematics of the transformations, the related system is self-affine or self-similar. Self-affine and self-similar natural systems are fractals, in the sense that any part of them is a small or large representation or imitation of the whole, however, at different scales. Due to this, a fractal system can be studied by focusing on its parts. In addition, the scaling properties of a fractal system is strongly related to its long memory [55][56][57] and its complexity [56], in the sense that a complex system is described by linear mechanisms and exhibits order [58,59]. Fractality, long memory and complexity is highly associated and, as a result, the analysis of the long memory of a system, yields the analysis of its fractal behavior and the delineation of its complexity and vice versa [1,3,4,17]. All these properties can reveal if strong links exist between the past, present and future of a system.
Among the different techniques that calculate the fractal properties of a system, the direct ones are more efficient. Concerning fractal dimension calculations the techniques of Katz, Higuchi and Sevcik provide very reliable estimations and for this reason they are employed in this paper. Fractal systems with long memory exhibit also power-law dependencies. These dependencies are outlined effectively with DFA and spectral power-law analysis. For this reason, these methods are employed here as well. From the above methods, DFA is considered to be the most robust of all. All methods can be compared via the related Hurst exponent. In the following sections, these methods will be described in detail. At the beginning Hurst exponent is presented and, thereafter, DFA, the methods for the calculation of fractal dimensions and the spectral power-law analysis.

Detrended Fluctuation Analysis (DFA)
During a system's long-memory behavior, the complex processes that are associated with the related phenomena, are described by long-range power-law associations, inconstant fluctuations, and behavior that is independent from scale [23,71]. Non-stationary features are embedded in the related time series, usually, associated with pseudosinusoidal patterns [22], repeated temporal trends [16], noise and other sources. If the related time series is non-stationary the usual methods cannot be used, e.g., the spectrum analysis and the techniques that are based on the autocorrelation of the time series [55,60,72].
DFA has been successfully applied in different subjects areas for which a behavior independent from scale is addressed. Characteristic examples are DNA [75], the dynamics of the heart [79,80], the day-to-day rhythms [81], meteorology [82], the variations of micro-climate temperature [83], economics [84], pre-seismic variations of radon in soil [14,16,17] and pre-earthquake activity [3,4,13,14,16,17,22,27,29,47]. Moreover, DFA as well as Hurst analysis have been shown in [43] to result in a better distinction of truly precursory signals from artificial noise if they are applied in the natural time domain [85]. The analysis in this new domain gives encouraging results in diverse fields including seismic electric signals [43] and cardiology [86].
At first the initial time signal is integrated. Then, the fluctuations, F (n), of the integrated signal are determined within a window of size n. Then the scaling exponent (self-similarity parameter), α, of the integrated time series is determined using a least-square fit to the linear transformation of log (F (n)) − log (n). Depending on the inherent dynamics of the system, the log (F (n)) − log (n) line may show one crossover at a scale n where the slope exhibits an abrupt change, two crossovers at two different scales n 1 , n 2 [16,18] or may not display a crossover at all.
DFA of a time signal in a single dimension, y i (i = 1, ...., N ), can achieved through the following process [3,4,16,17]: (i) First, the time series is integrated: In Equation (1), the symbol <...> indicates the total average value of the time series and k denotes the various time scales. (ii) Then, the integrated time series, y (k), is sub-divided into equal bins of length, n without overlapping. (iii) y (k) is then fitted to a function representing the trend in the box. Simple linear trends, or polynomials of 2 or higher order are used [3,4,16,17]. A linear function was employed here. The y coordinate of this linear function is symbolized as y n (k) in each box n. (iv) Then, the integrated time series y(k) is detrended. This is iterated in every box of length n, by subtracting the local linear trend, y n (k). In this way and for every bin, the detrended time series, y n d (k), is calculated as: (v) For every bin of size n, the root-mean-square (rms) of the fluctuations of the integrated and detrended time series is then calculated as where F(n) are the rms fluctuations of the detrended time series y n d (k). (vi) The procedure steps (i)-(v) are iterated for several sizes (n) of the scale boxes. This provides the type of link between F(n) and n. If there are long-term associations in the time series, the relationship between F(n) and n is exponential: In Equation (4), the scaling exponent α (DFA exponent) evaluates the power of the long-term associations of the time series. (vii) Via a logarithmic transformation of Equation (4), the linear relation between logF(n) and log(n) is determined the slope of which equals α. A good linear correlation indicates indicates the related fluctuations are long-lasting and, therefore, associated phenomenon has long memory.
In this paper, the goodness of the linear fit is quantified by the square of the Spearman's (r 2 ) correlation coefficient [3,4,16,17,49,87]. Good linear fits were considered those with r 2 ≥ 0.95.
The sliding window DFA was implemented according to the following steps [3,4,16,17,87]: (a) The time series were segmented in equal windows of 1024 samples each. This approximated one-month duration of the investigated segment of the time series; (b) A least-square fit of logF(n) versus log(n) was employed in every window in accordance to Equation (4). Following the approach of a recent paper of members of the team [88] the data were fitted to a straight line without seeking crossovers under the constraint that the slope of the fit exhibited square of Spearman's correlation coefficient above 0.95; (c) The window was forwarded one sample and the procedure (a)-(b) was iterated until the end of the signal; (d) DFA slopes α were finally plotted versus time and the corresponding plot data were extracted to ASCII output files for further use.

Katz's Method
The method of Katz can calculate the fractal dimension, D, of a time series. At first, the transpose array [s 1 , s 2 , ..., s N ] of the time series s i , i = 1, 2, ..., N, is determined, where s i = (t i , y i ) and y i represents the measured series values at the time instances t i [89,90].
Two subsequent points of the time series s i and s i+1 correspond to the value pairs (t i , y i ) and (t i+1 , y i+1 ), for which the Euclidean distance equals to: The total length of the curve that is generated from the distances of Equation (5) equals to: If this curve does not cross itself, it will extend in the planar to d, where: The fractal dimension, D, according to the Katz's method is calculated by Equations (5)-(7) as where n = L/a and a equals to the average value of the distances of the points.

Higuchi's Method
As with the Katz's method, the method of Higuchi also determines the fractal dimension D of a time series y(1), y(2), y(3), ..., y(N) that is recorded at intervals i = 1, 2...N [91,92]. For the application, the time series of Equation (9) are converted to a new sequence, y k m , which is constructed as follows [91][92][93]: According to [93], the length of the curve associated with the time series is given by: In Equation (11), m and k are integers that determine the time lag between the investigated samples and which are related as m = 1, 2...k. The symbol [...] in (11) denotes Gauss notation, i.e., the bigger integer part of the included value. The term is a normalization factor. For fractal curves of dimension D, the average value of lengths L(k) of Equation (7), exhibits a power law of the form: Therefore, from the linear regression of log-log transformation of L(k) versus k, k = 1, 2, ..., k max , the Higuchi' s fractal dimension, D, can be determined as corresponding slope. The reader should note here that the time intervals are k = 1, .., k max for k max ≤ 4, i.e., k = 1, 2, 3, 4, for k max = 4 and k = 2 ( j−1)/4 , j = 11, 12, 13..., for k > 4 (k max > 4) where [...] is notation of Gauss [90].

Sevcik's Method
The method of Sevcik estimates the fractal dimension D of time series as well. Following the method of Sevcik [94], the fractal dimension of a time series is approximated from the Hausdorff dimension, D h , of the related curve as [90]: In Equation (14), N( ) is the number of segments of length which make up the curve. If the curve has length L, then N( ) = L/2 [90] and, therefore D h can be written as: By employing a double linear transformation, the N points of the curve L can be corresponded to a unit square of N × N cells of the normalized metric space. With this transformation Equation (11) becomes [90,94] : The approximation of the fractal dimension according to the Sevcik method from Equation (16), improves as N → ∞.

Computational Methodology of Fractal Dimension
The fractal dimensions of the electromagnetic time series of this paper were calculated computationally according to the following methodology: (i) The time series was segmented in windows of 1024 samples each, i.e., of approximately 20 min span). (ii) In reference to each method, the fractal dimensions were calculated:

•
Katz's method: Equal to D of Equation (8) for n = 1024 and a = 1, a value that corresponds to the distance between the points of the series that constitute the parameter L and to the sampling rate of the electromagnetic time series (1 Hz).

•
Higuchi's method: Equal to the slope D of the first order least-square fit of the log-log transformation of Equation (8), namely the relation of log( L(k) ) versus log(k), for kmax = 16. (iii) Each window was forwarded one sample (sliding window technique) and the procedure (i)-(ii) was iterated until the end of the time series. (iv) Time-evolution plots of the fractal dimensions in accordance to the Katz's, Higuchi's and Sevcik's methods were generated, and the partial data were extracted to ASCII files for further use.
If a time series is a temporal fractal, its power spectral density, S( f ), will follow a power law of the form where f is a frequency of a transform. In this paper, and in accordance to the previous publications of the reporting team [9,13,14,20,21] this frequency was selected to be equal to the central frequency of the Morlet wavelet. In Equation (17), the exponent β evaluates the strength of the power-law connection whereas a (spectral amplification) quantifies the power of the contribution of each spectral component. By applying a log-log transform Equation (17) becomes: log S( f ) = log a + β · log f Equation (18) is a straight line and, hence, β and a can be calculated via the least-square fit to the corresponding data. As in previous publications [1][2][3][4]13,14,[16][17][18], the goodness of fit of the least-square fit of Equation (18) was quantified by the square of the Spearman's (r 2 ) coefficient under the constraint r 2 ≥ 0.95. In the above publications the method was also characterized as spectral fractal analysis. Hereafter, the simpler term fractal analysis will be employed.

Computational Methodology of Fractal Analysis
To implement fractal analysis of the electromagnetic time series of this paper the next methodology was followed: Similar approach in EM and radon time series [1,2,4,13] and in radon time series [3,13,14,[16][17][18].
(b) Class II: this class contains the windows of the time series segments with DFA's r 2 < 0.95 (i.e., they do not follow the prominent fBm class) or with 0 < α < 1 (i.e., they follow the fractional Gaussian noise (fGn) class).

Chaos Analysis Outcomes Comparisons
According to previous publications [4,[15][16][17], the results of the chaos analysis methods can be compared to each other, but the best approach is to compare all results through the Hurst exponent.
For the Class-I segments that are characterized by predictability and precursory value, the Hurst exponent (H) is calculated from the fractal analysis parameters as follows, e.g., [4]: (1) From (DFA exponent) (α) as: (2) From fractal dimension (D) as: (Berry's equation) (3) From power-law β as: It should be emphasized that according to extensive argumentation given in recent references [3,4,16], deviations are observed from the simple linear association of Equations (19)-(21) in the analysis from in situ measurements. As explained in the above publications, the relation between the chaos analysis parameters remains linear, possibly, of a slightly different type.

Meta-Analysis
The results from the application of all five methods of Sections 3.3-3.5 (DFA, Higuchi's, Katz's and Sevcik's fractal dimensions and spectral fractal analysis), are extracted in ASCII output for the purpose of meta-analysis in a two-stage procedure: Through this iterative procedure, 13 different combinations of techniques per five, four, three and two are generated. Similar procedure has been followed with success in a recent publication [49]. According to extended argumentation and discussion of recent publications [3,4,16], the important issue when analyzing pre-seismic time series to identify hidden pre-earthquake trends, is not just to locate some critical out-of-thresholds values, but rather to locate common areas with different techniques. When such common areas are found, the scientific evidence regarding the possibility of a pre-earthquake warning hidden in the time series, is increased and, hence a claim of pre-seismicity is stronger. Figure 3 presents the EM signal in parallel to time evolution of the DFA scaling exponent α and the evolution of the corresponding square of the Spearman's correlation coefficient. The profile of the DFA scaling exponent is completely different from the one of the time series. This has been acknowledged in previous publications of the reporting team [3,4,14,16,17,87] and is due to the fact that DFA identifies stationary and non-stationary patterns hidden in the time series with robustness [22,23,98]. Several DFA α exponent values lie in the Class-I value range (Section 3.6.1), namely the corresponding 1024-length segments are successful (Spearman's r 2 ≥ 0.95) fBm ones. As explained already, these segments are of notable pre-seismic precursory value [1][2][3][4]9,[13][14][15][16][17][18]20,21]. Figure 4c shows numerous segments with DFA exponents changing between 1.35 < α < 1.5 (anti-persistency) and 1.5 ≤ α < 2.0 (persistency). As mentioned in Section 3.6.1, these segments correspond to EM areas that are potentially associated with earthquakes of the period. On the contrary, non-successful (r 2 < 0.95) fBm segments, as well as fGn segments, are of low precursory value. Such low predictability and low-precursory areas are observed around 2018/10/29 (day 5 from start, day 0 at 26 October 2018) and 2018/11/04 (day 11 from start) and many other after 2018/11/18 (day 29 from start) and are the first that are neglected from the meta-analysis. The most important segments are the persistent Class-I ones. Several such Class-I segments are observed. Investigators (see e.g., the reviews of [10,11] and the references therein) have declared these segments as noteworthy signs of pre-earthquake activity. Several EM segments are spotted with distinct changes between anti-persistency (1.35 < α < 1.5) and persistency (1.5 ≤ α < 2.0). Several publications of the reporting team justify that these EM segments are of noteworthy pre-seismic precursory value, e.g., [1][2][3][4]9,[13][14][15][16][17][18]20,21]. Numerous persistent EM segments (1.5 ≤ α < 2.0) are observed. These EM segments are declared by others, e.g., [10,16] as undoubtedly footprints of pre-seismic activity. Similar observations as those of Figure 4, have been derived from pre-earthquake EM time series [1][2][3][4]9,16,17,20,21] and pre-earthquake time series of radon in soil [1,3,4,[13][14][15][16][17][18]. The meta-analysis in the ASCII outcomes of DFA (Section 3.8) is important, because it identifies computationally the segments which exhibit over-threshold DFA exponents. Accounting the argumentation given in Sections 3.3 and 3.8 and in the discussion of the related papers of the reporting team [3,4,10,11,14,16,17,87], the strict threshold of α = 1.6 is set for the DFA exponent, namely a threshold value bigger than the critical value of α = 1.5 which discriminates persistency from anti-persistency. Considering that acceptable DFA exponents are below or equal to 2.0, the corresponding value range becomes 1.6 ≤ α ≤ 2.0. The meta-analysis of the DFA slopes, yields a total of 22,943 DFA segments with acceptable persistent values between 1.6 ≤ α ≤ 2.0 at various intervals between 2018/10/26 (day 0) and 2018/11/25 (day 31). These segments correspond to critical long-lasting fractal epochs of the geo-system that generated the EM variations of Figure 3a. According to the argumentation presented so far, these EM time series segments recorded at the Ilia station are, most possibly, pre-seismic i.e., they are linked to earthquakes of the near area.  Figure 4 shows the time evolution of the fractal dimensions estimated by the methods of Katz, Higuchi and Sevcik. Noteworthy variations are observed. Deviations are also observed in the fractal dimension values calculated by the three algorithms. All discrepancies can be attributed to the different calculation approach of Katz's, Higuchi's and Sevcik's methods. Two recent publications [3,49] acknowledged that as well. Katz's and Higuchi's methods estimate higher fractal dimensions than the estimations of Sevcik's method. Fractal dimensions of Figure 4 and DFA exponents of Figure 3 can be associated from Equations (19) and (20) for precursory Class-I fBm segments as D = 3 − α ⇔ α = 3 − D. All fractal dimensions are within the value range calculated from relation D = 3 − α for Class-I segments (1.0 < α < 2.0), since this α-value range yields to fractal dimensions 1.0 < D < 2.0 as those of Figure 4. The opposite procedure (α = 3 − D), yields also to predictable Class-I DFA exponents in the range 1.0 < α < 2.0 from the results of the Katz's and Higuchi's methods and 1.5 < α < 2.0 from the results of Sevcik's method. The lower estimations of the Sevcik's method yield, hence, an estimation of pure persistent Class-I DFA exponents. EM Segments with distinct changes between anti-persistency (1.5 < D < 1.65) and persistency (1.0 ≤ D < 1.5) can be spotted for fractal dimensions calculated via the Katz's and Higuchi's methods. According to several publications these segments are of noteworthy pre-seismic precursory value, e.g., [1][2][3][4]9,[13][14][15][16][17][18]20,21]. Several EM segments corresponding to persistent behavior (1.0 ≤ D < 1.5) are observed with all methods of fractal dimension calculation. As aforementioned, these Class-I segments are considered by others, e.g., [10,16] as unambiguous pre-earthquake footprints. As with the meta-analysis of the DFA exponents, the corresponding meta-analysis threshold for the data of Figure 4, is of importance. For consistency with the meta-analysis of Figure Figure 5 presents the results from the fractal analysis method. As with Figures 3 and 4, the time evolution of power-law exponent, β, differs from the one of the time series. This is due to the fact that the fractal analysis identifies the fractal and long-memory trends hidden in the time series [1-4,9,13-17, 20,21,24,25,48,95-97]. Considering the points given in Section 3.2 the following categorization is valid for the comparison of the results of Figures 3 and 5, namely the comparison of DFA α and β exponents according to Equations (19) and (21): 1.

Results and Discussion
If 1.0 < β ≤ 3.0, the time series constitute a temporal fractal and follow the precursory Class-I category; • If 1.0 < β < 2.0, the time series are anti-persistent; • If 2.0 < β < 3.0, the time series are persistent;

2.
If −1.0 ≤ β < 1.0, the time series follow the Class-II category, i.e., they are of low predictability and precursory value; Especially: • If β = 1.0, the fluctuations of the processes do not grow and the related system is stationary; • If β = 2.0,the system follows random dynamics of no memory (random-walk); Most of the power-law values are successful (r 2 > 0.95) with β > 1.0 and, therefore, they correspond to predictable Class-I EM segments. Several successful segments are spotted with changes between anti-persistency and persistency. As emphasized already, the matched EM segments correspond to pre-seismic epochs of significant predictability and precursory value. Several segments have β > 2.
As analyzed above, the corresponding EM segments are considered by others [10,11], as precursory signs of the inevitable phase of the earthquake occurrence. Regarding the first phase of meta-analysis of the fractal analysis results, the threshold of β = 2.2 is set in accordance to the α = 1.6 of the DFA method (from Equations (19) and (21): β = 2 · α − 3). This β threshold is well above the value β = 2.0 which discriminates persistency from anti-persistency. With this threshold, a total of 62,294 EM segments are over-threshold Class-I segments of high predictability. These EM segments are associated with critical epochs of strong fractal behavior recorded by the Ilia station.  From the above argumentation, it can be supported that a significant number of EM disturbance segments recorded by the Ilia station between 2018/10/26 and 2018/11/26 are out of thresholds and bear significant signs of impeding seismic activity in the surrounding area due to the following reasons: (i) A total of 22,943 EM segments are persistent with α ≥ 1.6 according to the DFA. The robustness of DFA, its fundamental property to locate hidden long-memory trends in time series, together with its extensive use in studies pre-seismic activity from geosystems, e.g., [10,11,16,23,99], provide strong clues on the pre-seismic nature of the related EM segments. (ii) A significant portion of EM segments are below-threshold and recognized as signs of pre-seismic activity via three different fractal dimension calculation algorithms. A total of 564,082 are identified by the Katz's method, 142,725 with the Higuchi's method and 652,603 with the Sevcik's method. These segments are directly linked through relation D = 3 − α (Equations (19) and (20)) to several out-of-threshold EM segments identified from DFA. The out-of-threshold EM segments (common with DFA or not) have low fractal dimensions and high Hurst exponents both indicating high predictability of the related time series and significant precursory value of these segments as regards their pre-seismic nature. In addition, all fractal dimension algorithms have been used with success in radon in soil pre-earthquake disturbances [3]. (iii) A total of 62,294 EM segments are recognized as of high predictability and of significant pre-earthquake fractal nature according to the findings of the fractal analysis technique. The fractal methods are very important in the study of pre-earthquake geosystems, because these exhibit intense fractal activity, both in space and time, according to extensive literature reports, e.g., [8,10,16].
From the argumentation given so far and the logic of Sections 3.7 and 3.8, Figure 6 presents the stage 2 of the meta-analysis (subsection b of Section 3.8) in parallel to the earthquakes of Table 1 and Figure 2. As mentioned, Figure 6 presents all 13 possible combinations of fractal and long-memory methods (Section 3) per two, three, four and five techniques (stage 2 of meta-analysis) versus the 17 earthquakes of Table 1. Figure 6 is generated through GNU Octave ® based on the stage 2 meta-analysis results of all methods. It should be emphasized that the meta-analysis (Section 3.8) is a fully computational method and thus, Figure 6 is a computer-generated plot in accordance to the results of the meta-analysis. It should be noted though that Figure 6 is an effort to visually present altogether a great amount of data that correspond to all combinations of the 13 different ASCII files of 1 s −1 rate each. As observed from Figure 6, earthquakes 1, 2 and 12 are concurrent with the black '+' marks, which correspond to the combination of DFA versus all methods (combination of five methods). Earthquakes 7 and 13, 14 are almost concurrent with the black '+' mark. A computational search within the corresponding meta-analysis ASCII files of DFA versus all methods, shows that all six earthquakes (earthquakes 1, 2, 7, 12, 13, 14) emit pre-seismic signs close in time, from some hours to less than an hour before their occurrence and close in space (all earthquakes have epicenters close to the EM station with the maximum distance 100.5 km for earthquake 1). Regarding the kHz radiation other investigators [8,10,11], have claimed that the kHz radiation is emitted from days up to some hours prior to earthquake occurrence and that when emitted, the material's final catastrophe has started and the occurrence of the earthquake is inevitable. The findings of this paper for earthquakes 1, 2, 7, 12, 13 and 14 seem to support such an interpretation. In addition, earthquakes 3, 4, 5 and 6 are concurrent with the marks of all meta-analysis combinations except the black '+' one, namely they are concurrent with all chaos analysis methods combinations per four, three and two. Via these combinations, earthquakes 3, 4, 5 and 6 emit signs of pre-seismicity just before their occurrence. The computational search within all combinations of meta-analysis ASCII files shows that earthquakes 3, 4, 5 and 6 emit pre-earthquake disturbances some hours prior to their occurrence. Also, earthquakes 8, 9, 10, 15, 16 and 17 are concurrent with the black ' ' mark. This means that the combination of the methods of fractal analysis versus Higuchi's and Sevcik's (3 methods) support the view that earthquakes 8, 9, 10, 15, 16 and 17 emit pre-seismic signs shortly before their occurrence. According to the ASCII file of the stage 2 meta-analysis for the combination of the above methods, earthquakes 8, 9, 10, 15, 16 and 17 emit warnings some hours before their occurrence. Finally, earthquake 11 is concurrent with the magenta ' ' which corresponds to the combination of DFA versus Katz's and Sevcik's methods with similar interpretation.  Table 1 and Figure 2. The stem lines (black with green circle and red outline) correspond to the numbering of the earthquakes of Table 1 (1-17).
The above findings provide noteworthy evidence that all earthquakes of Table 1 and Figure 2, emit pre-seismic warnings from some hours to several minutes before their occurrence time according to the recordings of the 3 kHz antennas. It is very important that the validity of the outcomes presented above, is supported by the meta-analysis of combination of at least three techniques. Significant is also that the methodological approach with the meta-analysis, provides much stronger arguments than the results of any of the methods separated. Several other papers have published results from separated methods (one or more), e.g., [8,10,11]. Advanced approaches based on methods comparison and significant theoretical background have been published as well, e.g., [22,23]. The approach of the present paper is, however, quite different and novel. This is because it is based on well-studied methods with extensive use in analysis of pre-seismic activity, importantly, in combination. As analyzed in several parts of this paper, by combining more than two chaos analysis and long-memory methods, the scientific evidence is much stronger when investigating the aspect that some kind of recorded pre-earthquake activity is possibly associated with earthquakes of the same time period and geographical area. Re-organizing the results under this view, very strong clues support the aspect that pre-seismic 3 kHz EM disturbances are emitted prior to earthquakes 1, 2, 7, 12, 13 and 14, since all combinations of methods are matched. The arguments for the 3 kHz EM activity prior to earthquakes 3, 4, 5 and 6, are quite strong since the combination of four, three and two methods support the pre-seismicity of these. Noteworthy evidence is given for earthquakes 8,9,10,11,15,16 and 17 by three combinations of methods.
EM activity of the 3 kHz frequency emitted long before the occurrence of earthquakes (up to several days) has been reported also elsewhere [2,4,9,21]. Similar findings have been published for the MHz EM radiation [1,2,10,13,16,17,20,100] and for radon in soil variations [1,3,11,[13][14][15][16][17][18]. After-effects have been reported by other investigators [5,6,8,10,18, please see reviews of]. Under this view, the meta-analysis results of Figure 6 that do not coincide with the earthquakes of Table 1 and Figure 2, might be probably pre-seismic effects emitted well before these earthquakes or can be post-seismic effects. As discussed in several papers [1][2][3][4]9,[13][14][15][16][17][18]20,21,100], no one-to-one correspondence can be established between certain recorded activity and an ensuing earthquake. Moreover, as mentioned by [5], it is a serendipitous fact to record a very strong earthquake near a monitoring station and when found the evidence of an association are strong. Under these views, it is very hard to identify if the methods match in between the earthquakes of Table 1 and Figure 2, are post-seismic or pre-seismic. This is a limitation of the present methodology. On the other hand, the evidence for the precise or near matches are strong to noteworthy and this is a significant advantage of this methodology.
According to the argumentation given throughout the text, all earthquakes of this paper exhibited characteristic critical epochs of fractality and long-memory. As expressed in several publications [1][2][3][4]9,[13][14][15][16][17][18]20,21], these epochs can be linked to the propagation of micro-cracks and cracks during the preparation phase of these earthquakes. At this phase, the micro-cracks are generated continuously and form larger cracks in a self-organizing and fractal manner. In this way, the small cracks constitute small scale fractal imitations of larger cracks and hence generate effective pathways that allow the propagation of pre-seismic EM anomalies [3,4,100]. During this process, critical Class-I fBm-profile EM disturbances are addressed. The results of Figures 3-6, indicate this Class-I process since all critical segments exhibit persistent DFA exponents, fractal dimensions and power-law β-exponents. These are related to fBm modelling [3,4,16] which are produced by the 3 kHz EM generating geo-system of the Zakynthos area. During the preparation of the studied earthquakes, the focal area consists of a backbone of strong and large asperities that sustain the system and are modelled as fBm profiles. At a first stage, the fracture of the heterogeneous system in the focal area obstructs the backbone of asperities, but when the critical persistent meta-analysis matches of Figure 6 occur, the 'siege' of the asperities begin. Thereafter, the fracture starts and the unavoidable evolution of the process starts towards global failure. Finally, all critical warnings are of ensuing earthquakes and are revealed with the employed methods of this paper from the presented EM disturbances.

1.
This paper focuses on the post-seismic activity of a strong M L = 6.6 earthquake occurred on 2018/10/25 in Zakynthos Island, Greece. The post-seismic period extends over one month and is based on 3 kHz EM disturbance measurements derived by a ground-station located at Kardamas, Ilia, Greece. Seventeen earthquakes are included in the study with magnitudes between M L = 4.5 and M L = 5.5 and depths between 3 km and 17 km with all epicenters near Zakynthos Island and Ilia.

2.
Five different time-evolving chaos analysis methods are employed in the analysis. These methods are the detrended fluctuation analysis, the fractal dimension analysis with the methods of Higuchi, Katz and Sevcik and the power-law spectral fractal analysis. All these methods have been used with success in several pre-earthquake EM and radon signals in Greece.

3.
A novel fully computational methodology (meta-analysis) is applied to the time-evolution ASCII outcomes of all five chaos analysis techniques. Via a two-stage process, all out-of-threshold ASCII data values are computationally searched and the common time instances of 13 possible combinations of five, four, three and two techniques are noted. Through this process combination results of significant value are produced.

4.
Several persistent segments are found through DFA with exponents between 1.6 ≤ α ≤ 2.0. Higuchi's, Katz's and Sevcik's methods identify numerous segments with fractal dimensions 1.4 ≤ D ≤ 2.0. Many segments with 2.2 ≤ β ≤ 3.0 are recognized by the fractal analysis method. All these thresholds refer to persistent fBm Class-I segments of high predictability and pre-seismic value.

5.
Numerous combined meta-analysis segments are located with fractal behavior, dynamical complexity and long-memory. All these correspond to persistent fBm Class-I segments and are considered to be pre-earthquake footprints of high reliability. 6.
Six of the 17 post-earthquakes are matched by all 13 selected combinations of five, four, three and two chaos analysis methods. Four earthquakes are matched by all combinations of four, three and two methods from the 13 combinations. The remaining seven earthquakes are matched by at least one combination of three methods. Activity within typical time windows among or after these earthquakes is reported as well. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.