On Possible Electromagnetic Precursors to a Significant Earthquake (Mw = 6.3) Occurred in Lesvos (Greece) on 12 June 2017

This paper reports an attempt to use ultra-low-frequency (ULF) magnetic field data from a space weather monitoring magnetometer array in the study of earthquake (EQ) precursors in Greece. The data from four magnetometer stations of the HellENIc GeoMagnetic Array (ENIGMA) have been analyzed in the search for possible precursors to a strong EQ that occurred south of Lesvos Island on 12 June 2017, with magnitude Mw = 6.3 and focal depth = 12 km. The analysis includes conventional statistical methods, as well as criticality analysis, using two independent methods, the natural time (NT) method and the method of critical fluctuations (MCF). In terms of conventional statistical methods, it is found that the most convincing ULF precursor was observed in the data of ULF (20–30 mHz) depression (depression of the horizontal component of the magnetic field), which is indicative of lower ionospheric perturbation just 1 day before the EQ. Additionally, there are indications of a precursor in the direct ULF emission from the lithosphere 4 days to 1 day before the EQ. Further study in terms of NT analysis identifies criticality characteristics from 8 to 2 days before the EQ both for lithospheric ULF emission and ULF depression, while MCF reveals indications of criticality in all recorded magnetic field components, extending from 10 to 3 days before the EQ. Beyond the recordings of the geomagnetic stations of ENIGMA, the recordings of the fracto-electromagnetic emission stations of the HELlenic Seismo-ElectroMagnetics Network (ELSEM-Net) in Greece have been analyzed. The MHz recordings at the station that is located on Lesvos Island presented criticality characteristics (by means of both NT analysis and MCF) 11 days before the EQ, while a few days later (7–6 days before the EQ), the kHz recordings of the same station presented tricritical behavior. It is noted that the magnetosphere was quiet for a period of two weeks before the EQ and including its occurrence.


Introduction
There has been a significant progress in the field of electromagnetic (EM) phenomena associated with earthquakes (EQs), so-called seismo-electromagnetics, during the last few decades, and it is recently agreed that EM effects take place not only in the lithosphere, but also in the atmosphere and ionosphere before a major EQ [1][2][3][4]. The direct effect of EQs is EM emission in a wide frequency range from direct current (DC) to very high frequency (VHF) (or even higher), and the indirect effects of EQs are seismo-atmospheric and -ionospheric perturbations that are detected with the use of radio techniques, while one of the main objectives of seismo-electromagnetics is the elucidation of the lithosphere-atmosphere-ionosphere (LAI) coupling process, which is still controversial (e.g., [4]).
It is well known that the continuous monitoring of ultra-low frequency (ULF) (≤1 Hz) geoelectric signals (Seismic Electric Signals, SES) is in extensive operation for a long time in Greece (e.g., [5,6]), and that MHz-kHz fracto-electromagnetic emissions (fracto-EME or simply EME) are also being observed in Greece for long time [7][8][9][10]. Moreover, it has been reported that magnetic field variations accompany the appearance of SES signals [8]. However, to the best of our knowledge, an analysis of ULF magnetic field data from a space weather monitoring fluxgate magnetometer array within the framework of seismo-electromagnetics or EQ prediction studies has never been carried out in Greece, even though ULF magnetic field variations from such magnetometer arrays have been found to be of great importance in EQ precursor studies for short-term EQ prediction [2,3].
There are two types of anomalies in the ULF effects in possible association with EQs. The first is the well-known conventional ULF direct lithospheric emissions, which have been considered to be an important candidate of EQ precursors for a long time, since the 1988 Spitak EQ [11][12][13][14], reflecting the generation of electric currents due to the pre-EQ fracturing process in the lithosphere. On the other hand, the second is a rather recent phenomenon, named ULF depression, which appeared in the form of depression (decrease, reduction) of the horizontal component of the magnetospheric ULF Alfvén waves that were observed on the ground, which can be attributed to the changes in the lower ionosphere [15][16][17].
In the present article, we focus on the above-mentioned two types of ULF anomalies, as well as on MHz and kHz fracto-EME signals that were recorded prior to a rather strong EQ of magnitude Mw = 6.3, which occurred in Greece (south of Lesvos Island) on 12 June 2017 at 12:28:39 UT in a focal depth of 12 km. We find that multiple evidence corroborate the view that the specific EQ came with EM precursors.
Specifically, we first apply conventional statistical analysis to specifically defined daily-valued ULF quantities [15][16][17][18], as calculated from the magnetic field measurements of four magnetometer stations that are scattered around Greece. From this analysis, we find notable anomalies in the ULF depression, which are indicative of the lower ionospheric perturbation, at two stations, at the station closest to the epicenter and a more distant station, one day before the EQ. Moreover, there are indications of a precursor in the direct ULF emission from the lithosphere four days to one day before the EQ at the same two stations. In further investigating magnetometer data, we apply the natural time (NT) criticality analysis method [19] to the above-mentioned daily-valued ULF quantities for the two stations for which possible precursors to the EQ were identified. NT analysis identifies the criticality characteristics from eight to two days before the EQ both for lithospheric ULF emission and ULF depression. Additionally, we apply another criticality analysis method, referred to as the method of critical fluctuations (MCF) [20][21][22][23] to the raw (unprocessed) magnetometer recordings, revealing indications of criticality in all recorded magnetic field components, extending from ten to three days before the EQ.  The ENIGMA (HellENIc GeoMagnetic Array) [http://enigma.space.noa.gr/] ground-based fluxgate magnetometer network, as shown in Figure 1 by green triangles, acquired the geomagnetic data that were used in this paper. ENIGMA, which is operated by the National Observatory of Athens (NOA), is the first and it presented only permanent magnetometer network that ever operated in Greece and within a few years of operations achieved becoming a member of SuperMAG [http://supermag.jhuapl.edu/]. The array consists of four ground-based magnetometer stations in the areas of Trikala (Klokotos-THL), Attica (Dionysos-DIO), Lakonia (Velies-VLI), and Lasithi (Finokalia-FIN) in Greece that primarily targeted providing measurements for the study of geomagnetic pulsations, resulting from the solar wind-magnetosphere coupling. The respective coordinates, altitude, and instrumentation of the stations are summarized in Table 1. During the period of interest (1 May 2017-17 June 2017), ULF geomagnetic field data from all four magnetometer stations were available. The ULF data with a sampling frequency higher than 1 Hz were downsampled to 1 Hz. Ground-based magnetometers have proven to be the workhorse of magnetosphere-ionosphere coupling physics. The instruments that were used are survey vector (fluxgate) magnetometers, which means that three orthogonal sensors are required to measure the components of the magnetic field in all three dimensions, along with their fluctuations. ENIGMA measures the geomagnetic field in magnetic coordinates (i.e., H and D represent horizontal variations, while Z vertical variations of Earth's magnetic field). The ENIGMA magnetometer array enables the effective remote sensing of geospace dynamics and the study of space weather effects on the ground (i.e., Geomagnetically Induced Currents-GIC). ENIGMA contributes data to SuperMAG, which provides easy access to validated ground magnetic field perturbations in the same coordinate system, identical time resolution, and with a common baseline removal approach. The purpose of SuperMAG is to help scientists, teachers, students, and the general public have easy access to measurements of the Earth's magnetic field [24].
The here employed MHz and kHz signals were recorded by the ELSEM-Net (HELlenic Seismo-ElectroMagnetics Network) (http://elsem-net.uniwa.gr) ground-based fracto-EME station network, spanning across Greece, as shown in Figure 1 by blue circles (see also Table 2). Specifically, these were recorded at the Agia Paraskevi, Lesvos (M) station. The EM variations that are continuously monitored by ELSEM-Net at the Zante (Z) station correspond to narrowband (Q > 10) 3 kHz north-south, 3 kHz east-west, 3 kHz vertical, 10 kHz north-south, 10 kHz east-west, and 10 kHz vertical magnetic field, as well as 41, 54, and 135 MHz electric field recordings, while at the rest of the stations, including the Agia Paraskevi, Lesvos station, only the 3 kHz north-south, 3 kHz east-west, 10 kHz north-south, 10 kHz east-west magnetic, and 41, 46 MHz electric field are recorded. Both the kHz and MHz receivers are custom-made, providing, as outputs, the RMS values of kHz and the signal strength of MHz signals, after a <0.5 Hz low pass filtering. All of the receiver outputs are sampled at f s = 1 Hz. More details on the instrumentation of the fracto-EME stations of ELSEM-Net can be found in the supplementary downloadable material of [9].

Possible ULF EQ Precursors by Means of Conventional Statistical Analysis
In the following we try to detect any possible EQ related ULF anomalies by means of conventional statistical analysis. This can be defined as a method to identify any extreme values, i.e., local maxima or local minima of specific ULF quantities [15][16][17][18].
As the initial step of the ULF conventional statistical analysis, we have estimated the averaged night (00 h-04 h, UT) background spectra of H, D, and Z components for all magnetometers, after it was found that industrial interferences at this time interval are minimal. Figure 2 shows these spectra in the frequency band 10-100 mHz, averaged over the time period from 1 May to 17 June 2017, from which one observes that the spectra of the horizontal components are about the same. The background level is only slightly higher in THL. As expected, the vertical components present lower power spectral density (PSD) values when compared to the horizontal components.

Possible ULF EQ Precursors by Means of Conventional Statistical Analysis
In the following we try to detect any possible EQ related ULF anomalies by means of conventional statistical analysis. This can be defined as a method to identify any extreme values, i.e., local maxima or local minima of specific ULF quantities [15][16][17][18].
As the initial step of the ULF conventional statistical analysis, we have estimated the averaged night (00 h-04 h, UT) background spectra of H, D, and Z components for all magnetometers, after it was found that industrial interferences at this time interval are minimal. Figure 2 shows these spectra in the frequency band 10-100 mHz, averaged over the time period from 1 May to 17 June 2017, from which one observes that the spectra of the horizontal components are about the same. The background level is only slightly higher in THL. As expected, the vertical components present lower power spectral density (PSD) values when compared to the horizontal components. Based on these preliminary analysis results, we decided to try to detect any EQ precursors by using ULF data from all four stations: the closest to the EQ epicenter DIO as well as the other more distant stations to investigate the distance effect. Note that DIO is located at about 250 km, while the other three stations are located at about 400 km from the EQ epicenter (see Figure 1). Based on these preliminary analysis results, we decided to try to detect any EQ precursors by using ULF data from all four stations: the closest to the EQ epicenter DIO as well as the other more distant stations to investigate the distance effect. Note that DIO is located at about 250 km, while the other three stations are located at about 400 km from the EQ epicenter (see Figure 1).
When trying to detect possible ULF EQ precursors, we have to plot the temporal evolutions of geomagnetic and seismic activities in parallel to the studied ULF quantities, because we have to distinguish between the seismic and space weather effects. This information is plotted in the top panels of Figures 3-6, in which the Kp index indicates the geomagnetic activity and K LS (local seismicity) is the seismic index defined by K LS = 10 0.75M /(R + 100) (R: epicentral distance in km, M: EQ magnitude) [2,15,16]. Since a K LS value higher than 1 usually means a significant EQ that could have precursors, this 2017 Lesvos EQ is found to be a remarkable event.
In each one of Figure 3a,b, we plot a number of different ULF quantities for a 17-days-long, period starting from 1 June 2017, as calculated from the D, Z (Figure 3a) and the H, Z (Figure 3b) magnetic field components that were observed at VLI. These are the mean power of a horizontal magnetic field ∆T,∆ f (Figure 3a) or F h = H 2 ∆T,∆ f (Figure 3b), the mean power of the vertical magnetic field component in the third panel, F z = Z 2 ∆T,∆ f , and the ratio of vertical to horizontal mean power in the fourth panel, P z/d = F z /F d (Figure 3a) or P z/h = F z /F h (Figure 3b), [15][16][17][18]. In the fifth and sixth panel, respectively, we plot the absolute, Dep d (Figure 3a (Figure 3b), ULF depression (defined later) [15,16]. The above defined ULF quantities are estimated in a specific frequency band (∆ f ) averaged over the (local) midnight interval (∆T).
When trying to detect possible ULF EQ precursors, we have to plot the temporal evolutions of geomagnetic and seismic activities in parallel to the studied ULF quantities, because we have to distinguish between the seismic and space weather effects. This information is plotted in the top panels of Figures 3-6 M : EQ magnitude) [2,15,16]. Since a LS K value higher than 1 usually means a significant EQ that could have precursors, this 2017 Lesvos EQ is found to be a remarkable event.
In each one of Figure 3a   The quantities portrayed in the third and fourth panel in Figure 3a,b are considered to be a possible indicator of conventional lithospheric emission (e.g., [13][14][15][16][17]). While, the bottom two refer to the inverse of the horizontal component, Dep d or Dep h , and its relative variation, δDep d or δDep h . These two quantities are not so common among the scientific society and they need some explanation [15]. These indicate the depression (decrease, reduction) of the horizontal component j=i−N Dep h,j for the D or H component, respectively, where N is equal to the number of preceding days for averaging. The denominator indicates the average value and the numerator means the deviation from the mean. The above defined ULF quantities are known to indicate the absorption (or scattering) of downgoing Alfvén waves from the magnetosphere [2,17], and which is quite similar in nature to the lower ionospheric perturbations as detected by subionospheric VLF propagation anomalies [25,26]. Additionally, the statistical former result [15] indicated that the values of Dep and δDep are proportional to the magnitude of an imminent EQ.
Based on the analysis of the EM background (see, for example, Figure 2), we have adopted a time window (∆T) of 00-04 h UT (i.e., 02-06 h LT) with the lowest interferences, frequency band (∆ f ) of 20-30 mHz as the most remarkable band, and an averaging interval of N = 10 days for δDep d and δDep h . Figure 3 suggests very clear enhancements in Dep d and Dep h , as well as in δDep d and δDep h one day before the EQ (marked by red arrows). This is especially visible in D component ( Figure 3a). Moreover, the ratio of the vertical to horizontal mean power as the most important parameter for ULF lithospheric emission [18] shows an enhancement from four days (marked by blue arrows) to one day before the EQ for both horizontal components (P z/d , P z/h ).
Focusing on the possible lithospheric ULF emission, Figure 4 shows the temporal evolution of mean power ratios P z/d ( Figure 4a) and P z/h (Figure 4b) in all stations for the same time period as in Figure 3. This ULF quantity presents a remarkable increase in the data of both components for the nearest station DIO for four days before the EQ, while this increase, though noticeably less in level, is also observed at a more distant point in VLI. This indicates a probable observation of lithospheric signal. At the same time, this precursor was not detected in the rest distant stations, FIN and THL.
In more detail, it can be seen from Figure 4 at the station DIO closest to the EQ epicenter (distance about 250 km) that, during the whole period, there exists a significant enhancement, a double-peaked broad maximum from four days to one day before the EQ. On the other hand, for the farther station of VLI, we may find the corresponding maximum at this station as well, although at lower levels. As for the peak appearing one day before the EQ, it is possible that this enhancement in P z/d and P z/h may be related with the decrease in F d and F h (i.e., enhancement in Dep d and Dep h ) or local interferences on this day. However, the peak in both mean power ratios four days before the EQ is highly likely to be a precursor to this EQ. This ULF emission is a local effect and the previous works on ULF emission indicate that the nominal detection range is less than 100 km for the present EQ magnitude [13,14]. However, these results suggest that this effect is also likely to be observed over long distances.  identified just one day before the EQ mainly at the VLI and DIO stations. The station of VLI is about 400 km away from the EQ epicenter (see Figure 1), but this effect is most clearly noticed there. It is well known from our previous statistical and event studies that this ULF depression effect can be sensitively detected even at far distances from the EQ epicenter [15,16]. Of course, the EQ preparation zone may extend up to ~500 km (in radius) for this EQ magnitude [1,2]. It also has to be noted that the region of depression seldom is exactly above EQ epicenter. One cannot easily define the location of the region in the ionosphere that causes a depression. It is expected to be located over a place of geochemical changes, such as gas emanation. In order to suggest such a region, one should carefully combine information concerning the tectonics of Greece, as well as possible information regarding gas emanation observations. However, this is considered to be out of scope for the present work. In order to further investigate for possible precursors, we have plotted the temporal evolutions of the ionospheric ULF quantities, absolute ( Figure 5) and relative depression (Figure 6), for a much longer time period (about 1.5 months) from 1 May 2017 to 17 June 2017 for all sites and both of the horizontal components. One can observe that the most reliable detection of an anomaly, in both quantities of depression, is obtained for the D component that was observed at VLI, while much less reliable results are obtained for the other sites.
Specifically, a clear enhancement in Dep d ( Figure 5a) and the simultaneous very clear enhancement in δDep d (Figure 6a) (i.e., maximum depression of the horizontal component) can be identified just one day before the EQ mainly at the VLI and DIO stations. The station of VLI is about 400 km away from the EQ epicenter (see Figure 1), but this effect is most clearly noticed there. It is well known from our previous statistical and event studies that this ULF depression effect can be sensitively detected even at far distances from the EQ epicenter [15,16]. Of course, the EQ preparation zone may extend up to~500 km (in radius) for this EQ magnitude [1,2]. It also has to be noted that the region of depression seldom is exactly above EQ epicenter. One cannot easily define the location of the region in the ionosphere that causes a depression. It is expected to be located over a place of geochemical changes, such as gas emanation. In order to suggest such a region, one should carefully combine information concerning the tectonics of Greece, as well as possible information regarding gas emanation observations. However, this is considered to be out of scope for the present work. At this point, we have to comment on the important issue of the possible association of the above identified precursors with any space weather activity. The ENIGMA magnetometer array is, in general, orientated towards the study of space weather and extreme space weather events, such as magnetic storms. Nevertheless, during periods of low geomagnetic activity, its recordings can be utilized, as happens with other magnetometers around the world, for the tracking of signals of different frequencies, such as the ones that are discussed in the present study in order to detect the possible precursors of seismic activity. During the period of interest (1 May 2017-17 June 2017), the space weather conditions, and consequently the geomagnetic activity, were extremely low. The only exception was marked on the 28 May 2017 (15 days before the studied EQ), when a moderate At this point, we have to comment on the important issue of the possible association of the above identified precursors with any space weather activity. The ENIGMA magnetometer array is, in general, orientated towards the study of space weather and extreme space weather events, such as magnetic storms. Nevertheless, during periods of low geomagnetic activity, its recordings can be utilized, as happens with other magnetometers around the world, for the tracking of signals of different frequencies, such as the ones that are discussed in the present study in order to detect the possible precursors of seismic activity. During the period of interest (1 May 2017-17 June 2017), the space weather conditions, and consequently the geomagnetic activity, were extremely low. The only exception was marked on the 28 May 2017 (15 days before the studied EQ), when a moderate At this point, we have to comment on the important issue of the possible association of the above identified precursors with any space weather activity. The ENIGMA magnetometer array is, in general, orientated towards the study of space weather and extreme space weather events, such as magnetic storms. Nevertheless, during periods of low geomagnetic activity, its recordings can be utilized, as happens with other magnetometers around the world, for the tracking of signals of different frequencies, such as the ones that are discussed in the present study in order to detect the possible precursors of seismic activity. During the period of interest (1 May 2017-17 June 2017), the space weather conditions, and consequently the geomagnetic activity, were extremely low. The only exception was marked on the 28 May 2017 (15 days before the studied EQ), when a moderate geomagnetic storm took place. On that day, the Dst index (http://wdc.kugi.kyoto-u.ac.jp/dstdir/index.html), a proxy of the strength of the magnetospheric ring current, and thus an indicator of magnetic storm intensity lowered to −125 nT at 08:00 UT. This is also captured by the maximum value of the Kp index, as shown in the top panels of Figures 3-6. Clearly, this geomagnetic activity event precedes the identified precursors both in ULF depression one day before the EQ and ULF lithospheric emission 4-1 days before the EQ, and for this reason, these are unlikely to be of global origin (space-sourced disturbances). By means of Kp index, one may also find a simultaneous geomagnetic activity one day before the EQ, but this geomagnetic activity is too small to affect our conclusion (possibly just slightly reducing the ULF depression).

Possible ULF EQ Precursors by Means of Criticality Analysis
In the following, we try to identify any indications of critical dynamics in the ULF magnetic field data by means of two independent time series analysis methods that are known for their ability to uncover critical dynamics. These are the recently proposed methods referred to as the natural time (NT) analysis [19] and the method of critical fluctuations (MCF) [20][21][22][23], as briefly described in Sections 4.1 and 4.3, respectively. Both of them have been applied to ULF magnetic field data recorded prior to a number of EQs that occurred in Japan, revealing useful information. The first one is applied to the daily valued ULF quantities, as defined in Section 3 [27][28][29][30], while the second one is applied to the raw (unfiltered) magnetometer measurements [31,32].

Natural Time Analysis
The transformation of a time-series of "events" from the conventional time domain to the NT domain is performed by ignoring the time interval between events and only retaining its normalized order of occurrence [19]. Accordingly, the k-th event corresponds to a NT χ k = k/N, where N is the total number of successive events. The "energy" Q k of each event is retained. Subsequently, the transformed time series (χ k , Q k ) is studied. By defining p k = Q k /∑ N n=1 Q n , the energy of k-th event normalized by the total energy, a system is considered to approach criticality when the parameter converges to the value κ 1 = 0.070, and at the same time, both the entropy in NT, S nt = ∑ N k=1 p k χ k ln χ k − ∑ N k=1 p k χ k ln ∑ N k=1 p k χ k and the entropy under time reversal, S nt− , satisfy the condition S nt , S nt− < S u = (ln 2/2) − 1/4(≈ 0.0966), where S u stands for the entropy of a "uniform" distribution in NT [19,33]. This is the set of criteria for criticality, according to NT analysis, which is usually applied to time series [19], and it has been applied to different raw (unprocessed) EM recordings that are possibly related to EQ, such as the ultra-low frequency (≤1 Hz), Seismic Electric Signals (SES) [5,[34][35][36], and the fracto-EME MHz signals [8][9][10].
Moreover, NT analysis is also applied to daily-valued quantities, such as the ULF magnetic field quantities that are defined in Section 3 [27][28][29]37], and VLF subionospheric propagation quantities [30], although in such cases there is usually a limited number of available data, as happens with the NT analysis of foreshock seismicity. However, in such cases, criticality is checked differently following the paradigm of NT analysis of foreshock seismicity [33,34,36,38].

1.
Specifically, the temporal evolution of specific NT analysis parameters is studied by progressively including new events in the analysis and each time calculating κ 1 , S nt , S nt− , and D , based on the already included events. Note that D is the "average" distance D = |Π( ) − Π critical ( )| between the curves of normalized power spectra Π( ) = ∑ N k=1 p k exp(j χ k ) 2 ( is the natural angular frequency, = 2π ϕ, with ϕ standing for the frequency in NT, termed "natural frequency") of the evolving seismicity and the theoretical estimation of Π( ) for κ 1 = 0.070, Π critical ( ) ≈ 1 − κ 1 2 [19]. Of course, in each step the time series, (χ k , Q k ) is rescaled in the NT domain, since each time the k th event corresponds to a NT χ k = k/N, where N is the progressively increasing total number of the considered successive events; then, all of the parameters involved in the NT analysis are calculated for this new time series; this process continues until the time of occurrence of the main EQ event. In the resultant time evolution of κ 1 , S nt , S nt− , and D , criticality is considered to be truly achieved when, at the same time [19,34]: (i) κ 1 approaches κ 1 = 0.070 "by descending from above", (ii) S nt , S nt− < S u , (iii) D < 10 −2 , and (iv) since the underlying process is expected to be self-similar, the time of criticality does not significantly change by varying the "magnitude" threshold. Although the selection of thresholds is arbitrary in the case of the ULF magnetic field and VLF subionospheric propagation quantities (usually more than 20 threshold values equispaced between zero and a maximum threshold value that is larger than the 50% of the maximum value of the examined quantity are considered), if criticality conditions are met in close dates for more than one of the considered threshold values, and then this is considered to be an indication of the validity of the performed analysis. It may be worth mentioning that these criticality conditions, including κ 1 = 0.070, were not derived theoretically, but rather empirically from computing the κ 1 values for well-known phenomena and the actual data of seismicity in Greece.

NT Analysis of ULF Data
As already mentioned, NT analysis is applied to the daily-valued ULF quantities that are defined in Section 3 [27][28][29][30] for the time period 1 May-12 June 2017, while the criticality criteria are the ones mentioned in Section 4.1 for the case of foreshock seismicity and other cases of limited amount of available data. A typical example of the results that were obtained after the application of NT analysis is shown in Figure 7. It shows the temporal evolution of the NT parameters κ 1 , S nt , S nt− , and D for four threshold values of the ULF quantity F h at DIO. Each time that an F h value exceeds the corresponding threshold, a new event is included in the analysis and all of the studied NT parameters are re-calculated. Green patches indicate the time period during which NT analysis criticality conditions are satisfied for each threshold value. One can observe that, during these periods, κ 1 approaches the value κ 1 = 0.070 "by descending from above", S nt , S nt− < S u (≈ 0.0966), and D < 10 −2 , simultaneously. The four time periods are overlapping, intersecting on 4 June 2017, which indicates that critical state is truly achieved on that date.  After extensively studying the ULF quantities that are defined in Section 3 for DIO and VLI by means of the NT analysis, signatures of critical dynamics were identified in both stations from 4 June 2017 to 10 June 2017, i.e., from eight to two days before the 2017 Lesvos EQ (indicative results are shown in Figures 8-10). In summary, 1. concerning the mean power of horizontal, vertical components and their ratio at DIO, clear criticality indications were found for F h and F d on 4 June 2019 (see Figures 7 and 8a,b), as well as marginal evidence of criticality (just for one threshold value) on 4-7 June 2017 for P z/h ; 2. concerning the same ULF quantities at VLI, clear criticality indications were found for F d on 5-9 June 2019 (see Figure 9a,b). It is noted that clear criticality indications were also found on 25 May 2019 for P z/d (see Figure 9c,d) and P z/h ; however, this is more likely to be related with the moderate geomagnetic storm that took place on 28 May 2019 (15 days before the studied EQ); 3. concerning ULF depression at DIO, δDep d presents clear criticality indications on 7 June 2019 (see Figure 8c

Method of Critical Fluctuations
MCF is based on the theory of phase transitions and the principles of Physics of critical phenomena and starts with the hypothesis of a one-dimensional (1D) non-linear intermittent map of the form φ n+1 = φ n + uφ z n + n [20,21], where φ n is the order parameter value at time n, while the shift parameter n introduces a non-universal uncorrelated noise that is necessary for the establishment of ergodicity [21]. Note that the corresponding density ρ(φ) of such a map resembles, to a large extent, the order parameter distribution close to the critical point, which is characterized by a plateau region and a rapidly decaying tail. By analyzing a time series, supposedly conforming to the above-mentioned map, MCF's scope is to extract the scaling behavior for the order parameter fluctuations in time, if any, and estimate the critical exponent z (in NT analysis of seismicity, it has been found [39] that the order parameter fluctuations of seismicity are almost simultaneously minimized before major earthquakes [40], with the initiation of SES activities). Notice, that for thermal systems, the exponent z is connected with the isothermal critical exponent δ as z = δ + 1. The key idea behind the MCF is that criticality manifests itself by a power-law distribution of properly defined laminar lengths (defined in the following) l, P(l) ∼ l −p l [41], where the exponent p l is p l = 1 + 1 δ = z z−1 . More details on the theory behind the method can be found in [20][21][22][23].
In brief, the application of MCF can be described, as follows, while a detailed description of the step-by-step application algorithm can be found in [32]. The analysis starts with a time series excerpt of adequate length (>~5000 values) presenting, at least, local stationarity, for which the histogram of the order parameter φ (usually the original time series amplitude values play the role of the order parameter) is calculated, and from that histogram, an amplitude value (usually on the steepest slope) is determined as the fixed-point φ o , which will serve as the "start of laminar regions". Following that, for a number of different values within the φ amplitude range, which are called "ends of laminar regions" and are denoted as φ l , the distribution P(l) of the "laminar lengths" of each corresponding laminar region (φ o , φ l ) is calculated (one distribution per φ l value) and plotted on a log-log plot. Note that laminar lengths are the waiting times within each laminar region (φ o , φ l ), in other words, the number of successive φ-values obeying the condition φ o < φ < φ l . Finally, each one of these P(l) distributions is fitted using the function f (l) = p 1 · l −p 2 · e −lp 3 , leading to a set of exponents p 2 , p 3 for each laminar region.
If the exponent p 3 is zero, then the exponent p 2 is equal to the exponent p l . The relation p l = z z−1 suggests that the exponent p l (or p 2 ) should be greater than 1. On the other hand, from the theory of critical phenomena [42], it results that the isothermal exponent δ is greater than 1. Accordingly, from the relation p l = 1 + 1 δ , one obtains the condition 1 < p l (= p 2 ) < 2. In conclusion, the critical profile of the temporal fluctuations is restored for the conditions: p 2 > 1 and p 3 ≈ 0. As the system departs from the critical state, then the exponent p 2 decreases, while the exponent p 3 simultaneously increases, thus reinforcing, in this way, the exponential tail of the laminar lengths distribution. Note that a similar map with opposite signs for u and z, namely φ n+1 = φ n − uφ −z n + n , describes tricritical dynamics [23], which, in terms of MCF analysis, corresponds to the conditions: p 2 < 1 and p 3 ≈ 0. In the tricritical case, it has been shown that p l (= p 2 ) = z z+1 = δ+1 δ+2 .

MCF Analysis of ULF Data
As already mentioned, MCF analysis is directly applied to the 1 Hz sampled raw (unfiltered) magnetometer measurements [31,32]. An example of the results that were obtained during the successive steps of the application of MCF is shown in Figure 10. It refers to part of the Z component recordings of the magnetic field at DIO (see inset of Figure 10a) that was found to be stationary by checking its cumulative mean value and the corresponding standard deviation using nested time series excerpts of progressively wider length. By applying the MCF algorithm that is presented in detail in [32], it was found that it is necessary to add uniform noise to the original time series (after normalizing it in the range [0, 1]) for the establishment of ergodicity [21]. Specifically, uniform noise within the interval [−0.002, 0.002] was added and the resultant time series (Figure 10a) was considered to play the role of the order parameter in further applying MCF. From the histogram of this time series amplitude values (Figure 10b), the fixed point (start of laminar regions) φ o = 0.965 was determined according to the turning point method [43,44]. Subsequently, for a number of different ends of laminar regions, φ l , (see Figures 10b and 10a) the waiting times within each laminar region were calculated, and subsequently the distributions P(l) of these waiting times (laminar lengths) were plotted and fitted by the function f (l) = p 1 · l −p 2 · e −lp 3 . Figure 10c shows one such example for the laminar region φ o (= 0.965) < φ < φ l (= 0.986), while Figure 10d shows the exponents p 2 , p 3 that were obtained for the different laminar regions. From Figure 10c,d, one can see that the MCF criticality condition p 2 > 1, p 3 ≈ 0 is valid for a wide range of end point values, denoting that the distribution P(l) of laminar lengths follows a power law.

Possible Fracto-EME EQ Precursors
In parallel with the analysis of the magnetic field recordings of the ENIGMA array of magnetometers, the MHz and kHz fracto-EME recordings of the stations of the ELSEM-Net network were also analyzed in the search of possible EQ precursors. After the analysis of both MHz and kHz recordings at different stations, using multiple time series analysis tools, it was found that (preliminary results have been presented in the pre-publication [45]): 1. Critical signature was identified in the 41 MHz recordings of station M (located on the Island of Lesvos). Specifically, a~1.7 h long time window (06:23:20-08:03:20) UT on 1 June 2017 was found by MCF ( Figure 12) and NT analysis (Figure 13) to present critical characteristics.
2. A few days after the appearance of the criticality in the 41 MHz recordings, tricritical behavior was found by means of MCF in the kHz fracto-EME recordings of the same station (both at the 3 kHz and 10 kHz) from 5 June 2017 21:25:55 UT-6 June 2017 14:05:55 UT (~16.5 h long time window). An example for 3 kHz east-west is shown in Figure 14.
3. No very strong, avalanche like, precursory kHz fracto-EME signals of clearly high organization and persistency were possible to be validated (by means of different entropic/information metrics, Hurst exponent, detrended fluctuation analysis (DFA), spectral exponent analysis, fractal dimension etc.). This is probably due to the fact that the EQ of interest happened in the Sea as it is discussed in the following.
From Figure 12d, one can verify that the MCF criticality condition p 2 > 1, p 3 ≈ 0 is satisfied for a wide range of laminar regions, while from Figure 14d, it is obvious that the tricriticality conditions p 2 < 1 and p 3 ≈ 0 are satisfied for a wide range of end points φ l . It is noted that MCF analysis was performed in a similar way to the one that is presented in Section 4.4 for the ULF magnetic field recordings, e.g., Figure 10. On the other hand, NT analysis, as already mentioned in SubSection 4.1, was not performed in the same way as the NT analysis of the daily-valued ULF quantities that were analyzed in SubSection 4.2. In the MHz fracto-EME case, the criticality criteria according to NT analysis are that  On the other hand, NT analysis, as already mentioned in Section 4.1, was not performed in the same way as the NT analysis of the daily-valued ULF quantities that were analyzed in Section 4.2.
In the MHz fracto-EME case, the criticality criteria according to NT analysis are that parameter κ 1 converges to the value κ 1 = 0.070, and at the same time, both the entropy in NT, S nt , and the entropy under time reversal, S nt− , satisfy the condition S nt , S nt− < (ln 2/2) − 1/4 [8][9][10].
It should be noted that the MHz fracto-EME time series in most cases are not in the form of clearly distinguishable bursts. Therefore, there is not an easy way to define the involved fracto-EME events (the energy of which is attributed to the important NT analysis quantity Q k ) [8], because this involves the determination of a threshold value, corresponding to the background noise level, above which the amplitude is taken into account for the calculations. In this case, it has been proposed [8] that an exhaustive search should be performed for the determination of the threshold value. This way, one can exclude thresholds that may not lead to reliable κ 1 values, because of possible contamination by uneliminated noise [8]. The key idea behind this approach is that, if the time series excerpt under analysis presents criticality characteristics, then there should be at least one threshold value for which the above-mentioned NT criticality conditions are satisfied. If there is no such threshold, and then the specific excerpt does not present criticality [8].
From Figure 13, one can verify that for the threshold range marked by the vertical orthogonal magenta shaded area the NT analysis criticality conditions are satisfied, since κ 1 converges to κ 1 = 0.070, and at the same time, both S nt and S nt− satisfy the condition S nt , S nt− < (ln 2/2) − 1/4(≈ 0.0966). possible contamination by uneliminated noise [8].
The key idea behind this approach is that, if the time series excerpt under analysis presents criticality characteristics, then there should be at least one threshold value for which the above-mentioned NT criticality conditions are satisfied. If there is no such threshold, and then the specific excerpt does not present criticality [8].
From Figure 13, one can verify that for the threshold range marked by the vertical orthogonal magenta shaded area the NT analysis criticality conditions are satisfied, since 1  converges to  Our analysis reveals that first in the timeline appear critical features in the MHz fracto-EME, implying that the possibly related underlying fracture process that is involved in the preparation of the main shock is at critical state. The presence of the "critical point" during which any two active parts of the system are highly correlated, 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 fracto-EME precursory signals are emitted is a spatially extensive process. It is noted that, according to the four-stage model of EQ dynamics by means of fracto-EME, as suggested in [7], the pre-seismic critical MHz fracto-EME, which corresponds to the first stage of the specific model, is considered to originate during the fracture of the part of the Earth's crust that is characterized by high heterogeneity. During this phase, the fracture is non-directional and it spans over a large area that surrounds the family of large high-strength entities (asperities) that are distributed along the main fault sustaining the system. Note that, for an EQ of magnitude ~6.5, such as the 2017 Lesvos EQ of interest, the corresponding fracture process extends to a critical radius of ~120 km [46]. Thus, during this phase, the fracture process is extended up to the land of the neighboring islands. For the EQ of interest, the analysis also reveals that next in Our analysis reveals that first in the timeline appear critical features in the MHz fracto-EME, implying that the possibly related underlying fracture process that is involved in the preparation of the main shock is at critical state. The presence of the "critical point" during which any two active parts of the system are highly correlated, 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 fracto-EME precursory signals are emitted is a spatially extensive process. It is noted that, according to the four-stage model of EQ dynamics by means of fracto-EME, as suggested in [7], the pre-seismic critical MHz fracto-EME, which corresponds to the first stage of the specific model, is considered to originate during the fracture of the part of the Earth's crust that is characterized by high heterogeneity. During this phase, the fracture is non-directional and it spans over a large area that surrounds the family of large high-strength entities (asperities) that are distributed along the main fault sustaining the system. Note that, for an EQ of magnitude~6.5, such as the 2017 Lesvos EQ of interest, the corresponding fracture process extends to a critical radius of~120 km [46]. Thus, during this phase, the fracture process is extended up to the land of the neighboring islands. For the EQ of interest, the analysis also reveals that next in the timeline, after the critical MHz fracto-EME, appear tricritical fracto-EME (kHz in the specific case). These correspond to the second stage of the aforementioned model. Tricritical dynamics signify the departure from critical state, implying that the possibly related underlying fracture process that is involved in the preparation of the main shock evolves from the highly symmetrical and spatially expanded phase to a low symmetry, spatially focused on a preferred direction (along the fault) phase [10,23]. The observed EM silence that follows up to the EQ occurrence corresponds to the fourth stage of the aforementioned model.
As already mentioned, the crucial very strong avalanche, like precursory kHz fracto-EME signals, which are recorded in the tail of the preseismic kHz EM emission in the case of shallow EQs that happen inland (or near coast-line) and correspond to the third stage of the aforementioned four-stage model, were not recorded prior to the seismic event under study. It has been suggested that the lounge of the kHz EM activity shows the fracture of asperities sustaining the main fault that corresponds to the third stage of the four-stage model. In the case under study, the fracture process during this stage has been confined to a narrow area, i.e., along the main fault in a submerged area. This situation by itself justifies the observed absence. On the other hand, in terms of percolation theory, the "hydraulic threshold", x c , during which the transition impermeable-permeable occurs, as well as the "mechanical" or "damage threshold", x m (x m > x c ), during which the infinite cluster (IC) is formed and the solid disintegrates, precede the shear displacement along the fault plane (the EQ), which means that many transport properties are activated before the main event [47], and references therein. This fact, when combined with the governing role of local hydraulic conditions on the water injection-induced fracture behavior in rocks concerning the mechanical properties, fracture nucleation, and the geometry of the shear fracture zone [47], further enhances the absence of the final strong pulse, like kHz fracto-EME. We note that, due to the crucial character of this emission, an austere set of criteria have been established to characterize such a recorded kHz anomaly as a seismogenic one [7,48,49]. The multidisciplinary analysis of the kHz records before the under study seismic event in terms of these criteria verified the absence of the last emerged part of the fracture induced kHz fracto-EME corresponding to the fracture of asperities.

Summary and Conclusions
This paper probably reports the first use of ULF magnetic field data from a space weather monitoring magnetometer array in Greece in the study of seismo-electromagnetics or EQ precursor studies. Our target was a rather big EQ occurred on 12 June 2017 at Lesvos Island in Greece with magnitude of 6.3.
First, the data from four stations of the Greek ULF network ENIGMA have been examined, and the detailed analyses that are based on the conventional statistical method have been performed for a time period from 1 May to 17 June 2017. We can summarize the obtained observational results on the precursors to this EQ as follows: (1) In order to find any ULF precursors, we found that the reliable time interval was nighttime as usual and the most important frequency range was 20-30 mHz (0.02-0.03 Hz), as already confirmed before [13][14][15][16]18].
(2) The most reliable and conspicuous precursor was ULF depression (most enhanced depression of the horizontal components of magnetospheric ULF waves), which was detected one day before the EQ. The more reliable results were obtained at a farther station of VLI (distance about 400 km). This kind of high sensitivity of the ULF depression phenomenon provides a further confirmation of our previous findings [15][16][17] and this ULF depression is likely to be attributed to the perturbation of the lower ionosphere [17].
(3) The direct ULF emission from the lithosphere as the result EQ preparation process is likely to be detected at the closest station of DIO (distance~250 km from EQ epicenter) in the form of enhanced in P z/d and P z/h . This effect is also observed, albeit not so clearly, at the more distant station of VLI (distance~400 km from EQ epicenter) and it does not appear in the rest two distant stations THL and FIN. This ULF lithospheric emission is known to be a much more local effect than the above ULF depression. The ULF lithospheric emission appearing at DIO four days before the EQ is convincing, but that appearing one day before the EQ is not so reliable because it might be related with the above ULF depression one day before the EQ. The lead time of four days seems to be much smaller than the conventional lead time of 2-3 weeks [13,14].
(4) The geomagnetic activity during the whole period of analysis seems to be rather quiet; that is, the maximum geomagnetic activity appeared in the end of May and the geomagnetic activity was quiet when the above precursors were detected.
Based on these results, we further investigated the ULF magnetic field data by applying criticality analysis for the two stations (DIO and VLI) for which possible precursors to the EQ were identified. The results that were obtained by two independent criticality analysis methods, NT analysis and MCF, revealed several indications of critical dynamics eight to two days and ten to three days before the EQ, respectively. Specifically: (1) Clear criticality indications were identified in the mean power of the horizontal components of the magnetic field (ULF quantities F h and F d ) on 4 June 2019 at DIO and on 5-9 June 2019 (in F d ) at VLI by means of NT analysis. While marginal evidence of criticality were found on 4-7 June 2017 for P z/h at DIO.
(2) Clear criticality indications were also identified in terms of NT analysis in the ULF depression quantities In order to find out whether other EM effects that were possibly related to the specific EQ took place, the data from the ground-based fracto-EME station network ELSEM-Net spanning across Greece have been analyzed, leading to the following results: (1) The MHz recordings of 1 June 2017 (11 days before the EQ) at the station located in Lesvos Island, very close to the EQ epicenter, presented criticality characteristics (by means of both NT analysis and MCF).
(2) A few days later, on 5 June 2017 (7-6 days before the EQ), the kHz recordings of the same station presented tricritical behavior by means of MCF.
(3) No very strong, avalanche like, precursory kHz fracto-EME signals of clearly high organization and persistency was possible to be validated (by means of different entropic/information metrics, Hurst exponent, DFA, spectral exponent analysis, fractal dimension, etc.). This is probably due to the fact that the 2017 Lesvos EQ happened under the Sea, as discussed in Section 5.
The above summarized results suggest that three types of possible precursory signs appeared before the 2017 Lesvos EQ: (a) noticeable change in the mean power of appropriately defined ULF magnetic field quantities that are possibly related to either direct lithospheric emissions or ionospheric depression, (b) approach to critical state for some of the abovementioned quantities as well as for raw magnetic field measurements and MHz fracto-EME, and (c) tricritical dynamics for kHz fracto-EME.
We note that the criticality indications that have been identified for the different EM signals all appear a few days before EQ occurrence. The EQ preparation process has various facets that reflect corresponding different EM signals. Indifferent to the mechanisms by which they are produced, the EQ-related MHz fracto-EME, direct lithospheric ULF emissions, and ionospheric ULF depression seem to be compatible with each other in terms of criticality. These EM precursors emerge during a "critical epoch", when the "short-range" correlations evolve into "long-range" ones, suggesting a spatially extensive process for their generation. For example, the generation of a preseismic ionospheric anomaly requires physical and chemical transformations that occur in a spatially extended preparation (activation) zone of an impending EQ. Such a requirement is satisfied during the appearance of the "critical epoch".
The obtained results show that the ULF magnetic field data of the ENIGMA magnetometer array, beyond their primary use for the study of space weather and extreme space weather events, such as magnetic storms [50], can be efficiently employed for the detection of possible precursors of seismic activity by applying advanced analysis methods. Moreover, the fusion of the information obtained from these ULF magnetic field data with information from other kinds of EM data, such as subionospheric VLF/LF (low frequency) propagation data (ionospheric perturbations), ionosonde data, satellite data, fracto-EME (MHz, kHz) signals, SES signals, etc., but also gas emanation data, such as radon concentrations in boreholes, as well as tectonic information, seismic zonation, and foreshock seismic activity data, may: (a) help deciphering the physical mechanisms that are related to different seismo-electromagnetic phenomena (e.g., elucidating the LAI coupling mechanism) and (b) contribute to a future multi-parametric seismic risk assessment system.