Criticality features in ultra-low frequency magnetic fields prior to the 2013 M 6 . 3 Kobe earthquake

The nonlinear criticality of ultra-low frequency (ULF) magnetic variations is investigated before a particular earthquake (EQ) occurred in Kobe on April 12, 2013, by applying the “natural time” analysis on a few ULF parameters: Fh, Fz and Dh. The first two refer to radiation from the lithosphere, and the last parameter corresponds to depression of horizontal component as a signature of ionospheric perturbation. A recent paper of our team has indicated, using the same data as in this paper but by means of conventional statistical analysis, a clear effect of depression in the horizontal component as an ionospheric signature. But there seems to be no convincing signature of lithospheric ULF radiation according to the specific analysis, so this paper aims at extending our study on the electromagnetic data recorded prior to the specific EQ by trying to find any significant phenomenon in ULF effects (both lithospheric radiation and the depression of horizontal component) using the critical, natural time analysis. The natural time analysis has yielded that criticality at Shigaraki (SGA), as the station closest to the EQ epicenter, is reached on March 27-29 for Fh and March 27 to April 1 for Fz (about two weeks before the EQ). But, the criticality for Dh was not observed at SGA probably due to high noise, on the other hand such criticality was observed at Kanoya (KNY) because of its known property of a wider range of detection of ULF depression.


Introduction
Based on extensive studies during the last few decades it is recently plausible that electromagnetic (EM) phenomena do appear prior to an earthquake (EQ) [Hayakawa and Molchanov 2002, Pulinets and Boyarchuk 2004, Molchanov and Hayakawa 2008, Hayakawa 2009, Hayakawa 2012, Hayakawa 2013, Eftaxias and Potirakis 2013, Eftaxias et al. 2013].Such possible EQ precursors include lithospheric phenomena such as DC geoelectric field, ultra-low-frequency (ULF) radiation, fracto-EM MHz -kHz emissions, and seismoatmospheric and -ionospheric perturbations.There are already a few EQ precursor signatures, which seem to be statistically correlated with EQs.One is geoelectric signals, for which Varotsos [2005] found a close correlation of seismic electric signals (SESs) with EQs.The other is perturbations of the ionosphere, both in the lower D/E layer and in the upper region (F region).Using sub-ionospheric VLF/LF propagation data during a 7-year observation, Hayakawa et al. [2010] have established a significant statistical correlation between the lower ionospheric perturbations and EQs with magnitude greater than 6.0 and with depth smaller than 40-50 km.Perturbations in the ionospheric F region as a change in f o F 2 are also confirmed to be correlated with EQs based on long-term data [Liu 2009].
On the contrary, the number of events is not so abundant for conventional lithospheric ULF electromagnetic emissions, even though they are found to be very promising for short-term EQ prediction [Hayakawa et al. 2007, Fraser-Smith 2009, Kopytenko et al. 2009].
The first possible ULF events were observed for the Spitak EQ in 1988 [Kopytenko et al. 1990, Molchanov et al. 1992, Kopytenko et al. 1993], while evidence of ULF signature for the 1989 Loma Prieta EQ may have been found [Fraser-Smith et al. 1990, Bernardi et al. 1991].Both of these EQs were huge with magnitude around 7. Then, ULF emissions in possible association with the 1993 Guam EQ (M=8.0) were found [Hayakawa et al. 1996, Kawate et al. 1998].A number of other studies reporting lithospheric ULF electromagnetic emissions have also been published [e.g., Kopytenko et al. 1994, Molchanov and Hayakawa 1998, Kopytenko et al. 2001, Ismaguilov et al. 2001, Ismaguilov et al. 2002, Kopytenko et al. 2003, Kopytenko et al. 2007, Molchanov 2011, Kopytenko et al. 2012].Though there have recently been published a few papers casting a doubt that those ULF emissions were not seismogenic, but just an effect of geomagnetic storms [e.g., Campbell 2009], their arguments were not adequate to deny the presence of seismogenic ULF emissions.Later ULF studies have been summarized in Hattori [2004], Hayakawa et al. [2007], and Molchanov and Hayakawa [2008], while recently Hattori [2013] has summarized ULF emissions observed with the Kanto (Tokyo) ULF network during the long term of 2000-2011, with a statistical finding of significantly higher probabilities of ULF anomalies before an EQ than after the EQ.
It is a fact that geomagnetic storms, magnetosphere and ionosphere radiations are the main sources of the magnetic field variations observed by groundbased magnetic observatories, and may have significant contribution to the ULF band too.However, there is a list of published articles (some of them already cited above) reporting that, during the preparation period of a strong EQ, ULF (<10 Hz) lithospheric magnetic emissions with noise-like character originate in a forthcoming EQ hearth.The intensities of the lithospheric ULF disturbances are very weak, but when we can observe these emissions, the epicentral distances can extend up to some hundreds of km [Hayakawa et al. 2007, Fraser-Smith 2009] before strong earthquakes (M >6).Possible generation and propagation mechanisms have also been proposed/studied [e.g., Molchanov and Hayakava 1998, Kopytenko et al. 2001, Hattori 2004, Molchanov and Hayakawa 2008, Surkov and Hayakawa 2014].
There have been proposed an extensive number of signal processing techniques to identify precursor signatures of seismic events, thereby increasing the number of convincing ULF events.We can list only a few here.The first one is the use of polarization (the ratio of vertical to horizontal magnetic field component) [Hayakawa et al. 1996] to distinguish between local seismogenic ULF emissions and those of space (magne-tospheric) origin.Another example is "direction finding" of different principles (either goniometer [Hayakawa et al. 2007], gradient method [Kopytenko et al. 2009], etc.) to locate the source of ULF emissions and to compare the source location with the epicenter of a future EQ.More studies on the development or application of new sophisticated signal processing techniques, like the one presented in this paper, should be pursued in the future in order to acquire more evidence before one can accept or reject a ULF anomaly as a possible EQ precursor.
An additional new phenomenon in ULF magnetic field changes is the non-conventional depression (reduction of fluctuations) in the horizontal ULF magnetic field (its value is estimated as the inverse of the average power of the horizontal magnetic field component, please also see Section 3) before an EQ [Schekotov et al. 2006, Schekotov et al. 2013a, Schekotov et al. 2013b].Seismogenic ULF electromagnetic waves mentioned above are "emission or radiation" from the lithosphere, but this non-conventional effect is vice versa: the depression (or a decrease) in amplitude of ULF down-going waves of magnetospheric origin.The mechanism of this ULF depression is not fully understood at the moment, but a few hypotheses have been proposed [Schekotov et al. 2006, Molchanov andHayakawa 2008].For any hypothesis, such ULF depression is likely to be related with the perturbation in the lower ionosphere [Schekotov et al. 2013b], such as those observed by sub-ionospheric VLF/LF propagation anomalies [Hayakawa et al. 2013a].
In addition to the above signal processing techniques used to obtain definite EQ precursor signature, there is another direction of signal processing based on the nonlinear critical evolution or self-organized criticality.We will be able to monitor microfractures which occur in final break in the local zone of an EQ by looking at the conventional lithospheric ULF emissions [Hayakawa et al. 2007, Eftaxias andPotirakis 2013].Furthermore, using the same ULF magnetic changes it may be possible for us to study as well the criticality in ULF magnetic field depression as an indicator of seismo-ionospheric signatures.
This kind of nonlinear pre-EQ fracture evolution (either lithospheric radiation or ULF depression) can be monitored by a few possible methods: One is fractal analysis [Hayakawa et al. 1999, Smirnova et al. 2001, Gotoh et al. 2003, Kapris et al. 2004, Ida et al. 2005, Potirakis et al. 2012], and the other is "natural time" analysis recently developed by Varotsos [2005], which has been applied to time series data of a point process like SES events [Varotsos et al. 2005] as well as pre-fracture lithospheric MHz electromagnetic emissions [Potirakis et al. 2013, Potirakis et al. 2015a, Potirakis et al. 2015b].After the recent successful application of the natural time method on the ULF parameters related to the 2011 Tohoku EQ [Hayakawa et al. 2015], this paper is an attempt to apply this natural time method to our ULF geomagnetic data observed in Japan during a period of a few months prior to a particular EQ occurred in Kobe on April 12 (UT=20h 33m), 2013.This EQ was selected as a case study on the reason that this is a significant EQ (magnitude M=6.3) in the western part of Japan where the seismic activity is not so high, even after the 2011 Tohoku EQ, while it took place at an area very close to the 1995 Kobe EQ.Our aim is to investigate the pre-EQ fracture process of the EQ, and, in turn, to quantitatively examine whether there exists a clear precursor in various ULF parameters, especially since the conventional statistical study yielded no clear precursor.The results obtained here are compared with the ones already published by Schekotov et al. [2015] in terms of conventional statistical analysis applied on the same ULF data.

ULF geomagnetic data and EQ information
ULF geomagnetic data are available from three reference stations of Kakioka (KAK), Memambetsu (MMB) and Kanoya (KNY) belonging to Japan Meteorological Agency.An additional station near the EQ epicenter was established by Kyoto University, which was kind enough to provide us with ULF data at Shigaraki (SGA) in Kyoto prefecture.The magnetometer at Kakioka station is of the fluxgate type, with a constant sampling rate of 1 sample per sec ( f s =1 Hz).In the frequency range below 0.1 Hz, which is the band of our interest, the frequency responses of amplitude (gain) and phase are very flat.The ULF sensors at the three other stations of Memambetsu, Kanoya and Shigaraki are of the same type, with the same sampling frequency and all have nearly the same characteristics (http://www.kakiokajma.go.jp/).The data are provided in the conventional format of IAGA (International Association of Geomagnetism and Aeronomy) 2000, where the magnetic field is represented by four time series: horizontal (H) component, declination (D), vertical component (Z), and total field (F).The geographic coordinates of the considered ground-based stations are (36°13'56'' N, 140°11'11'' E) for KAK, (43°54'36'' N, 144°11'19'' E) for MMB, (31°25'27'' N, 130°52'48'' E) for KNY and (31°51'9'' N, 136°6'11.4''E) for SGA, respectively.
The EQ treated here is the one occurred on April 12, 2013 (UT=20h 33m) (April 13, 5h 33m a.m.JST), 2013 KOBE EQ: CRITICALITY IN ULF MAGNETIC FIELDS Figure 1.The relative locations of ULF magnetic observatories (Shigaraki -SGA, Kakioka -KAK, Memambetsu -MMB, and Kanoya -KNY, all in black boxes) and EQs with M ≥ 5.5 (as circles with their size indicating the EQ magnitude and their color referring to depth).Date format: dd/mm/yy.with M=6.3 and depth of 15 km.This EQ was an inland, fault-type one taken place in the land of Awaji-island, and the epicenter of this EQ was located at the geographical coordinates (34°25.1'N, 13°449.7'E).We know that there are many faults in this area, including Rokko-Awaji fault region, Senzan fault region, etc. [Matsuda 1995] and the 1995 Kobe EQ ( January 17, 1995, M=7.2) [Nagao et al. 2002] was the consequence of movement of the Nojima fault belonging to the Rokko-Awaji fault region.The EQ studied in this paper is estimated to be generated at the south-west edge of aftershock region of the 1995 Kobe EQ.
Figure 1 illustrates the locations of three ULF observatories (KAK, MMB and KNY) indicated by black rectangle and the additional station of SGA (again by a black rectangle).The epicenter of the recent EQ of our interest is indicated by a red circle just below the sign "SGA"; being located west of SGA we call it tentatively Kobe EQ, with epicentral distance of ~50 km.Unfortunately, during the time period of our interest, there happened a large number of EQs with considerable magnitudes (with M≥5.5) in the eastern part of Japan (which are also plotted by circles in Figure 1), in which the circle size is proportional to magnitude and its color refers to depth.
By analyzing the raw magnetic field data at differ-ent times of the day in terms of short Fourier Transform spectrograms, we observed that lowest industrial noise is emitted during specific nighttime hours, see for example Figure 2. Based on this analysis of the background electromagnetic noise at the locations of the ground-based geomagnetic stations we decided to use the following interval T=03h (±0.5h)JST (=UT+9h) for our analysis.The date/time is indicated in UT hereafter.For the chosen time interval, we further studied the relation between the level of the ULF parameters of interest (especially of the level of the depression of horizontal magnetic field component, cf.Section 3) and the level of the interferences (electromagnetic noise) in terms of Fourier Transform spectral analysis concluding that the best possible signal to noise ratio, and accordingly maximum useful dynamic range, is achieved for the narrow frequency band 0.005~0.01Hz (5-10 mHz).Therefore, before any analysis, for each day we extract the abovementioned one hour long excerpt of the raw magnetic field time series and we narrowband (5-10 mHz) filter it in the time-domain.

Parameters of ULF magnetic fields to be analyzed
As will be seen in the next section, we know that the concept of "natural time" is suitable for time series data of point processes, including seismicity, DC geo- electric potential etc.Because the DC geoelectric potential observed by the Varotsos group is known to exhibit a succession of rectangular waveforms, this natural time analysis is well suited to those geoelectric signal data.On the other hand, ULF data are not a point process, but they are continuous in nature, so that we have to think of what kind of ULF parameters should be used for the natural time analysis.
Here we are not interested in the raw ULF data, but we apply the following pre-processing procedure on the ULF data to be used for the natural time analysis as our first step, so we obtain one datum per day [Schekotov et al. 2013a], starting from the raw H and Z time series as provided in IAGA 2000 format: (i) First, we extract the data corresponding to the one hour long time series excerpt around 18:00 (UT) and remove any possible spikes and gaps from the time series excerpt.The resulting one hour long time series excerpts for the horizontal and vertical component of the geomagnetic field are denoted H DT and Z DT , respectively.
(ii) Second, we perform a zero-phase digital filtering (in the time-domain) by processing the extracted time series excerpt data in both the forward and reverse directions by means of a 3rd order bandpass Butterworth filter (passing the band 5-10 mHz) obtaining the narrowband filtered one hour long time series excerpts H DT,f and Z DT,f .
(iii) Third, we calculate the square of each amplitude value of the one hour long filtered data, which are proportional to the instantaneous power values: H 2 〉, we define the following ULF variables on which we will apply the natural time analysis: (1) Average power of the horizontal magnetic field component: The parameters F h and F z are indicators of conventional ULF radiation from the lithosphere.It should be mentioned at this point that an increase of the ULF parameters F h , F z cannot be simply related to lithospheric electromagnetic emissions, since such an increase may be the result of a magnetic storm.This is why in their use for the study of ULF variations of possible lithospheric origin, one should first try to focus on time periods when ULF variations from other sources (e.g., from magnetic storms) are not expected.This is why we always study geomagnetic activity (at least Dst index) in search for such phenomena, in parallel to the study of the above mentioned ULF parameters.On the other hand, as already mentioned in Section 1, there have been proposed an extensive number of signal processing techniques to distinguish between local seismogenic ULF emissions and those of space (magnetospheric) origin.In the case of methods, like the natural time method presented in this article, that detect characteristics in the analyzed ULF parameters F h , F z which may be features of either magnetoshperic or lithospheric phenomena, the possible origin of the revealed features should be carefully examined.
The third ULF parameter is indicator of another phenomenon: the non-conventional parameter D h is an inverse of the average power of the horizontal magnetic field component, and is used in order to investigate the depression of ULF waves (of magnetospheric origin) observed on the ground as the ionospheric signature.Schekotov et al. [2006] studied this effect based on the data during rather long periods (4 years in Russia and two years in Japan), which was proved to be an important parameter in identifying a precursor to an EQ.This phenomenon is not completely understood, but it can be interpreted in terms of enhanced absorption of down-going Alfvén waves through the perturbed lower ionosphere [Hayakawa et al. 2013b], indicating the presence of perturbations in the lower ionosphere.Note that, both Z and H components of magnetic field may decrease due to the phenomenon of the seismo-ionospheric depression, so one could think that could employ either 1 〉, for the study of this phenomenon.However, we should pay attention to the minimum levels of the fields which are limited not only by their nature but also industrial interferences, noise of sensors and other causes.Due to variety of reasons the vertical component Z is less suitable for this purpose.Sensors of both components have equal noise but the vertical field component has lower values, leading to narrower dynamic range for a possible ULF parameter that would be based on the 〉, based on the horizontal component of the magnetic field H is used for the study of seismo-ionospheric depression.
Figure 3 illustrates temporal evolutions of some physical parameters at the station of SGA: (from top to bottom) geomagnetic (Dst index) and seismic (Kls index) activities, Kp and AE indices, F h , F z , D h , and ULF rela-tive depression, dDep, which for the i-th date is defined as ), with N being the filter parameter equal to the number of preceding days for the involved averaging; here N=10 (not analyzed in terms of the natural time ), all in the same frequency range.
The Dst index is an index of magnetic activity derived from a network of near-equatorial geomagnetic observatories that measures the intensity of the globally symmetrical equatorial electrojet (the "ring current").It is a function of magnitude, M, and epicentral distance (in km), R, namely, Kls =10 0.75M /(R+100) [Molchanov and Hayakawa 2008], while an EQ with Kls>1 is significant.
Based on the comparison of this figure at SGA with the corresponding figures at KAK and KNY, Schekotov et al. [2015] have concluded that there is no clear evidence of the generation of lithospheric ULF radiation, but a very clear depression (in terms of the relative depression parameter dDep) is observed a few days before this EQ especially at SGA.We can see from Figure 3 that there was a moderate geomagnetic storm on March 17 with the variation of about −150 nT, while we can also identify a number of more weak substorms.Nevertheless, there weren't any noticeable magnetic field disturbances observed in the vicinity of the studied earthquake.Note also that the March 17 moderate geomagnetic storm is distant enough in time so as not to influence the natural time analysis results associated with the EQ in focus of our study which took place on April 12, 2013.However, the analysis results presented in Section 5 are evaluated taking also into account the specific magnetic storm.The corresponding variations on F h and F z are detected clearly at KAK and KNY (though not shown as figures), but to a much less extent than at SGA as in Figure 3.
Only by looking at the temporal evolutions of F h and F z [Schekotov et al. 2015], we cannot come to any definite conclusion on the absence (or presence) of lithospheric ULF radiation.Then, in order to answer this question, we will apply the natural time analysis to those ULF parameters.Further we want to confirm our previous finding on the precursor of ULF depression [Schekotov et al. 2015] with the help of natural time analysis.

Natural time analysis method
The natural time method was originally proposed for the analysis for a point process like DC or ultra-low frequency (≤1Hz) SES [Varotsos et al. 2002, Varotsos 2005], and has been shown to be optimal for enhancing the signals in the time-frequency space [Abe et al. 2005].The transformation of a time-series of "events" from the conventional time domain to natural time domain is performed by ignoring the time-stamp of each event and retaining only their normalized order (index) of occurrence.Explicitly, in a time series of N successive events, the natural time, | k , of the k th event is the index of occurrence of this event normalized, by dividing by the total number of the considered events, | k =k/N.On the other hand, the "energy", Q k of each k th event is preserved.We note that the quantity Q k represents different physical quantities for various time series: for EQ time series it has been assigned to a seismic energy re-leased (e.g., seismic moment) [Varotsos et al. 2005], and for SES signals that are of dichotomous nature it corresponds to SES pulse duration [Varotsos 2005], while for MHz electromagnetic emission signals that are of non-dichotomous nature, it has been attributed to the energy of fracto-electromagnetic emission events as defined in Potirakis et al. [2013].The transformed time series (| k ,Q k ) is then studied through the normalized power spectrum P (s)=|∑ N k=1 p k exp( js| k )| 2 , where s is the natural angular frequency, s =2r{, with { the natural frequency, and p k =Q k /∑ N n=1 Q n corresponds to the k th event's normalized energy.
The study of P (s) at s close to zero reveals the dynamic evolution of the time series under analysis.This is because all the moments of the distribution of p k can be estimated from P (s) at s → 0 [Varotsos et al. 2011a].Aiming to that, by the Taylor expansion P (s)=1−l 1 s 2 +l 2 s 4 +..., the quantity l 1 is defined, where , the variance of | k weighted for p k characterizing the dispersion of the most significant events within the "rescaled" interval (0,1].Moreover, the entropy in natural time, S nt , is defined [Varotsos et al. 2006] as and corresponds [Varotsos et al. 2006, Varotsos et al. 2011b] to the value at q=1 of the derivative of the fluctuation function fl(q)=〈| q 〉−〈|〉 q with respect to q (while l 1 corresponds to fl(q) for q=2).It is a dynamic entropy depending on the sequential order of events [Varotsos et al. 2006].The entropy, S nt− , obtained upon considering [Varotsos et al. 2006] the time reversal T, i.e., Tpm = p N−m+1 , is also considered.
A system is considered to approach criticality when the parameter l 1 converges to the value l 1 =0.070 and at the same time both the entropy in natural time and the entropy under time reversal satisfy the condition S nt , S nt− <S u =(ln2/2)−1/4 [Sarlis et al. 2011], where S u stands for the entropy of a "uniform" distribution in natural time [Varotsos et al. 2006].
In the special case of natural time analysis of foreshock seismicity [Varotsos et al. 2001, Varotsos et al. 2005, Varotsos et al. 2006, Sarlis et al. 2008], the seismicity is considered to be in a true critical state, a "true coincidence" is achieved, when three additional conditions are satisfied: (i) The "average" distance 〈D〉 between the curves of normalized power spectra P (s) of the evolving seismicity and the theoretical estimation of P (s), P critical (s) = (18/5s 2 ) − (6coss/5s 2 ) − (12sins/5s 3 ), P critical (s) ≈ 1−l 1 s 2 for l 1 =0.070 should be smaller than 10 −2 , i.e., 〈D〉=〈|P(s) − P critical (s)|〉<10 −2 (this is a practical criterion for signaling the achievement of spectral coincidence) [Varotsos et al. 2011b]; (ii) the parameter l 1 should approach the value l 1 =0.070 "by descending from above", i.e. before the main event the l 1 should gradually decrease until it reaches the critical value 0.070 (this rule was found empirically) [Varotsos et al. 2001, Varotsos et al. 2011b]; (iii) Since the underlying process is expected to be self-similar, the time of the true coincidence should not vary upon changing (within reasonable limits) either the magnitude threshold, M thres , or the area, used in the calculation.
It should be finally clarified that in the case of seismicity analysis, the temporal evolution of the parameters l 1 , S nt , S nt− , and 〈D〉 is studied as new events that exceed the magnitude threshold M thres are progressively included in the analysis.Specifically, as soon as one more event is included, first the time series (| k ,Q k ) is rescaled in the natural time domain, since each time the k th event corresponds to a natural time | k =k/N, where N is the progressively increasing (by each new event inclusion) total number of the considered successive events; then all the parameters involved in the natural time analysis are calculated for this new time series; this process continues until the time of occurrence of the main event.
Note that in the case of natural time analysis of foreshock seismicity, the introduction of magnitude threshold, M thres , excludes some of the weaker EQ events (with magnitude below this threshold) from the natural time analysis.On one hand, this is necessary in order to exclude events for which the recorded magni-tude is not considered reliable; depending on the installed seismographic network characteristics, a specific magnitude threshold is usually defined to assure data completeness.On the other hand, the use of various magnitude thresholds, M thres , offers a means of more accurate determination of the time when criticality is reached.In some cases, it happens that more than one time-points may satisfy the rest of natural time critical state conditions, however the time of the true coincidence is finally selected by the last condition that "true coincidence should not vary upon changing (within reasonable limits) either the magnitude threshold, M thres , or the area, used in the calculation."

Natural time analysis results of ULF parameters
We consider each daily value which is above a certain threshold as an event.In our ULF cases (F h , F z , and D h ), the "energy" of k th event that is the value of the quantity Q k , is considered to be equal to the corresponding value of each one of above quantities, provided that this is above a certain threshold such as F h,thres , F z,thres , and D h,thres , respectively.We follow here the way of application of the natural time analysis on ULF parameters which has recently been suggested in Hayakawa et al. [2015], checking for the corresponding criticality criteria as presented in Section 4 for the case  (a-d) variations of the natural time analysis parameters (l 1 , S nt , S nt− , and 〈D〉) for the different thresholds F h,thres (in arbitrary unit) 0.0000, 0.0025, 0.0050, and 0.0075, respectively.Note that the events employed depend on the considered threshold.Moreover, the time (x-) axis is not linear in terms of the conventional date of occurrence of the events, since the employed events appear equally spaced relative to x-axis as the natural time representation demands, although they are not equally spaced in conventional time.Date format: dd/mm.Figure 5. Natural time analysis of the ULF electromagnetic emission quantity F z at SGA for the time period from February 1 to April 12, 2013: (a-d) variations of the natural time analysis parameters for the different thresholds F z,thres 0.00000, 0.00250, 0.00375, and 0.00500, respectively.Note that the events employed depend on the considered threshold.Moreover, the time (x-) axis is not linear in terms of the conventional date of occurrence of the events, since the employed events appear equally spaced relative to x-axis as the natural time representation demands, although they are not equally spaced in conventional time.Date format: dd/mm.Figure 6.Natural time analysis of the ULF electromagnetic emission quantity D h at SGA for the time period from February 1 to April 12, 2013: (a -d) variations of the natural time analysis parameters for the different thresholds D h,thres 0, 50, 150, and 200, respectively.Note that the events employed depend on the considered threshold.Moreover, the time (x-) axis is not linear in terms of the conventional date of occurrence of the events, since the employed events appear equally spaced relative to x-axis as the natural time representation demands, although they are not equally spaced in conventional time.Date format: dd/mm.seismicity.The analysis starts from a specific date and all natural time analysis parameters (l 1 , S nt , S nt− , and 〈D〉, cf.Section 4) are calculated from the time-series rescaled in the natural time domain each time a new event is added.The analysis stops at the day of main shock on April 12, 2013.This way, calculation of the involved natural time analysis parameters is repeated as many times as the total number of the events revealed through the thresholding phase, while finally an equal number of sets of these parameters are obtained and their evolution in natural time is plotted and assessed.Although the selection of thresholds involved is arbitrary (usually more than 20 threshold values equispaced between zero and a maximum threshold value 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, then this is considered to be an indication of the validity of the performed analysis.This is because the underlying process is expected to be self-similar.
Results of natural time analysis for different ULF parameters at the station of SGA, the closest to the EQ epicenter, are presented here.Figure 4 illustrates the natural time analysis results for F h at SGA. Figures 4a to 4d indicate the corresponding results for increasing the threshold (F h,thres ) from smaller value (4a) to higher value (4d).The corresponding results for F z at SGA are shown in Figure 5, and those for D h at SGA are plotted in Figure 6.As already mentioned, the four critical parameters l 1 , S nt , S nt− , and 〈D〉 are calculated each time a new event is included in the analysis, so their evolution as a function of date of new event occurrence is presented.Criticality conditions are satisfied in some cases, as observed in these figures.Though our target EQ is the one on April 12, there were found two more possible effects in the ULF criticality study: one is another EQ (M = 6.3, depth ~10 km) on February 25 taken place close to KAK, and the other is a geomagnetic storm on March 17 (with Dst =−150 nT).Figures 4, 5, and 6 are the results only at SGA, but Table 1 summarizes those at other stations (KNY and KAK) for the April 12 EQ.
Because the true coincidence is of essential importance in the natural time analysis, we pay particular attention to the fact that the time of true coincidence should not vary when changing the threshold.First of all, we summarize the result for the April 12 EQ as follows.

Lithospheric ULF radiation
Criticality for F h at SGA is observed on April 2 and critical characteristics are kept for some days up to April 10.For higher thresholds this criticality seems to be identified in slightly earlier dates .Also the criticality for F h is reached at KNY, and further certain criticality is maybe related to this Kobe EQ even at KAK. Whereas, the criticality for F z is observed at SGA, but on a little bit dates of March 27 to April 1, as well as at KAK on April 3-8, or even a bit earlier (starting from March 29) for higher thresholds.No corresponding criticality is observed at KNY.

ULF depression
As related to ionospheric signature, the criticality for D h seems to be reached on April 3-7 at KNY, but it is questionable that the criticality is marginally reached on April 6 at KAK.The detection of ULF depression at KNY and KAK is likely to be consistent with our previous findings that this phenomenon can be detected at extremely far distances (up to few thousands km before earthquakes with magnitude M ~8-9 and up to one thousand km before earthquakes with magnitude about 6) [Schekotov et al. 2013a, Hayakawa et al. 2013c, Hayakawa and Schekotov 2014].Whereas, there is detected no criticality at SGA even though it is closest to the EQ epicenter.At this point, we have to mention that out of the three examined stations, the SGA station was the one carrying the higher noise contamination.Therefore, we consider that the fact that no criticality was revealed for D h at SGA may be a result of the high noise level there.
The fact that criticality has been revealed for F h at the same station is not inconsistent with this interpretation.Note that although D h is closely linked to (is the inverse of ) F h , they are used for the study of two different phenomena, seismo-ionospheric depression (corresponding to decrease of the average power of the horizontal component) and seismo-related lithospheric emissions (corresponding to increase of the average power of the horizontal component), respectively.Importantly, the "events" of these two ULF parameters considered during the natural time analysis are not linked through a simple relation as it holds for D h and F h values.First of all, there does not necessarily exist an event D h at the same (conventional) time when an F h event exists and vice versa, while their time series carry different dynamics.On the other hand, the influence of noise is not the same to the events of these two parameters.The D h events are more influenced since the percent increase of D h during a preseismic period (compared to the values of D h in a quiet period) is lower than the percent increase of F h , while the dynamic range available for the study of D h is limited by the electromagnetic background emissions level (natural or anthropogenic) and sensor noise.Note that the higher (maximum) observed depression, D h , values which are candidate "events" for natural time analysis correspond to the lower (minimum) levels of the average power of the horizontal component.
In order to avoid the confusion, we will here briefly indicate some criticality effects for two other events: the February 25 EQ and a geomagnetic storm on March 17.By looking at Figures 4, 5 and 6, it seems that criticality is reached on February 19 for F h and probably on February 21-24 for F z at SGA, but much more evident criticality (though not shown) is observed for F h at KAK on February 22-24 and also for F z at KAK on February 15-24.The conspicuous behavior at KAK can be reasonable by its relative location with respect to the EQ epicenter.Next, we move on to the criticality possibly associated with the March 17 geomagnetic storm.
Among various parameters at different stations, we have observed some criticality only at KNY on March 4-7 on F h and at KAK on March 15 on D h .This may indicate that criticality is not so clearly detected for any geomagnetic storm of this size in Dst.

Conclusion and discussion
Using exactly the same data as in the present work (Figure 3 and corresponding figures at other stations, though not shown), we have recently published a paper [Schekotov et al. 2015], in which the following conclusions have been obtained for the April 12 Kobe EQ based on conventional statistical analysis.
(1) As for the conventional topic of lithospheric ULF radiation, it was quite uncertain whether there may exist any definite signature on the presence of lithospheric ULF radiation.
(2) On the other hand, as related to the non-conventional ULF effect of depression of horizontal magnetic field, a clear anomaly was found on the depression at all ULF stations.
This work aimed at extending our study on the electromagnetic data recorded prior to the specific EQ by trying to find any significant phenomenon in ULF effects (both lithospheric radiation and the depression of horizontal component) by using the critical natural time analysis method.The main outcome of the present work is that the natural time analysis has suggested new information on criticality in the lithosphere, improving thus our insight into the physical processes which took place prior to the occurrence of the examined Kobe EQ.The here obtained result that the lithospheric ULF radiation (and accordingly the lithosphere) came to critical condition prior to the main seismic event provides a much clearer idea of the involved process, compared to our recently published study [Schekotov et al. 2015] in which we could not obtain any definite conclusion about it.Specifically, our results show that criticality was reached at the station of SGA close to the EQ epicenter from the end of March to April 10 for both F h and F z , but the clearest criticality was detected in the beginning of April (about 10-12 days before the This lead time seems to be consistent with previous findings by Hayakawa et al. [2007] and Hattori [2013].However, we have to note that the ULF parameters natural time analysis results for this EQ are not as striking as those for the huge 2011 Tohoku EQ [Hayakawa et al. 2015].This is likely to be due to that the magnitude of this EQ is much smaller than that of the 2011 Tohoku EQ, and correspondingly the possible lithospheric ULF radiation seems to be much smaller in amplitude for this April 12 Kobe EQ.In other words, the natural time analysis seems to be not so effective when the signal to noise ratio is not high enough.
Next we discuss the non-conventional ULF depression effect for the April 12 EQ.As already shown in our recent paper by Schekotov et al. [2015], the depression of ULF horizontal magnetic field is also confirmed by the natural time analysis of D h showing critical behavior on April 3-7 prior to the EQ, which seems to be consistent in temporal evolution with Schekotov et al. [2015].This lead time is also found to be in agreement with our previous works by Schekotov et al. [2006Schekotov et al. [ , 2013a]].Such a clear peak in the conventional statistical analysis [Schekotov et al. 2015] is found to be confirmed with the natural time method that it is definitely due to the criticality.
We should also mention at this point that different ULF quantities were also investigated.Specifically, we applied the natural time analysis on the quantities median (H 2 DT,f ) and median (Z 2 DT,f ), corresponding to the median of the instantaneous values of power of the considered one hour long narrowband filtered time series excerpts, as well as on the quantity 1/median (H 2 DT,f ), in an attempt to study the statistical characteristics of the instantaneous power values as in Regi et al. [2015].The natural time analysis of the above median-based quantities surprisingly yielded results compatible to the ones obtained for F h , F z , and D h , in the sense that criticality characteristics were also revealed within the same time frame of about two weeks before the EQ.This finding leads to the interesting conclusion that criticality signature was not only embedded in the energy content focusing ULF parameters F h , F z , and D h , but also in the new median-based quantities focusing on the distribution of the corresponding instantaneous power values.
So, we can conclude that such a critical analysis as the natural time method would be of great importance, providing additional new information on what is happening physically in the lithosphere, especially when the conventional statistical analysis is unsuccessful in giving us any definite answer on the presence of lithospheric ULF radiation.
Finally we will comment on the further findings for two events as a by-product of this paper: February 25 EQ and a geomagnetic storm on March 17.Some criticality was detected at the station of KAK, probably in its close proximity to the epicenter of the February 25 EQ.Concerning the geomagnetic storm, the only clear satisfaction of criticality conditions seems to be on March 15 on D h at KAK and on March 4-7 on F h at KNY.The effect of a geomagnetic storm of this size is not so evident in the critical analysis of ULF magnetic variations, though we have found before some self-organized criticality of the geomagnetic storm in ULF variations before the geomagnetic storm with the use of fractal analysis [Smirnova and Hayakawa 2007].

Figure 2 .
Figure 2. Indicative background electromagnetic activity at the ground-based station of Shigaraki: Daily spectrum (power spectral density) on February 2, 2013, for the H-component recorded at the specific station.
Dst index reflects possible "geomagnetic" interferences which are not connected with local magnetic field disturbance caused by seismic activity.It gives us possi-bility to separate one from another.The Kp index is obtained from a number of magnetometer stations at mid-latitudes.It reflects global geomagnetic activity.AE index is an auroral electrojet index obtained from a number of stations distributed in local time in the latitude region that is typical of the northern hemisphere auroral zone.It reflects local magnetic disturbances which is known as magnetospheric substorms and may have durations of tens of minutes to several hours.All indices correlate rather good during periods of noticeable disturbances as evident from top and second panels of Figure 3.The index of local seismicity, Kls, should not be confused with the geomagnetic K or Kp indices.

Figure 3 .
Figure 3. Temporal evolutions at SGA of geomagnetic (Dst in nT) and seismic (Kls) activities (top panel), Kp and AE indices (2nd panel), the average power of the horizontal magnetic field component (F h ) (3rd panel), the average power of the vertical magnetic field component (F z ) (4th panel), the depression of horizontal magnetic field component (5th panel) and the relative depression dDep (bottom panel).The frequency is 0.005-0.01Hz (5-10mHz).A red vertical line refers to the moment of the 2013 Kobe EQ.Date format: dd/mm/yyyy.

Figure 4 .
Figure 4. Natural time analysis of the ULF electromagnetic emission quantity F h at SGA for the time period of February 1 through April 12, 2013: (a-d) variations of the natural time analysis parameters (l 1 , S nt , S nt− , and 〈D〉) for the different thresholds F h,thres (in arbitrary unit) 0.0000, 0.0025, 0.0050, and 0.0075, respectively.Note that the events employed depend on the considered threshold.Moreover, the time (x-) axis is not linear in terms of the conventional date of occurrence of the events, since the employed events appear equally spaced relative to x-axis as the natural time representation demands, although they are not equally spaced in conventional time.Date format: dd/mm.

Table 1 .
Summary of natural time analysis results for F h , F z , and D h at three magnetic stations (SGA, KAK, and KNY) for the April 12 EQ.