A tentative model for the explanation of Ba˚th law using the order parameter of seismicity in natural time

Using the order parameter of seismicity deﬁned in natural time, we suggest a simple model for the explanation of Ba˚th law, according to which a mainshock differs in magnitude from its largest aftershock by approximately 1.2 regardless of the mainshock magnitude. In addition, the validity of Ba˚th law is studied in the Global Centroid Moment Tensor catalogue by using two different aftershock deﬁnitions. It is found that the mean of this difference, when considering all the pairs mainshock-largest aftershock, does not markedly differ from 1.2 and the corresponding distributions do not depend on the mainshock’s magnitude threshold in a statistically signiﬁcant manner. Finally, the analysis of the cumulative distribution functions provides evidence in favour of the proposed model.


Introduction
According to Båth law (Båth 1965;Lombardi 2002;Varotsos et al. 2011c), the magnitude difference between a mainshock and its largest aftershock is approximately 1.2 and this value is independent of the mainshock magnitude. That means, in a catalogue of seismic events, the mean hDMi of the difference between the magnitudes of each mainshock M ms and its respective largest aftershock M max as is given by the following relation: In this paper, we first propose (in Sect. 2) a simple model in order to understand the origin of Båth law. This model makes use of the concept that the occurrence of a mainshock can be considered as a critical point and employs the order parameter proposed for seismicity on the basis of the new time frame termed natural time (Varotsos et al. 2011c) to which we now turn.
For a time-series consisting of N events, the natural time can be defined (Varotsos et al. 2001(Varotsos et al. , 2002(Varotsos et al. , 2011c as v k ¼ k=N and it serves as an index for the occurrence of the k-th event of energy Q k . In natural time analysis, we study the evolution of the pair (v k ; Q k ) or equivalently of the pair (v k ; p k ) where p k ¼ Q k = P N n¼1 Q n . Focusing on the analysis of seismicity, the evolution of the pair (v k ; M 0k ) is considered (Varotsos et al. 2004;Sarlis et al. 2010a), where M 0k denotes the seismic moment of the k-th event since it is a measure of its energy. The normalized power spectrum was introduced (Varotsos et al. 2001(Varotsos et al. , 2002 as and p k ¼ M 0k = P N n¼1 M 0n . The Taylor expansion of PðxÞ leads to: and hf ðvÞi P N k¼1 p k f ðv k Þ, where f ðvÞ is a function of the real variable v (Varotsos et al. 2007, e.g. when f ðvÞ ¼ v 2 we have hf ðvÞi ¼ hv 2 i ¼ P N k¼1 v 2 k p k . It has been shown (Varotsos et al. 2011b) that j 1 becomes equal to 0.07 at the critical state for a variety of dynamical systems. In particular, for the seismicity j 1 converges to 0.07 a few days before a mainshock (Varotsos et al. 2001;Sarlis et al. 2008) when starting the computation of the j 1 values after the initiation of seismic electric signals (SES) activities (Varotsos and Lazaridou 1991), i.e. natural time is set to 0 at the initiation of the SES activity. SES are the transient changes of the electric field of the Earth generated by means of a mechanism of point defects (Varotsos et al. 1993;Varotsos 2008) and that have long been successfully used for short-term earthquake prediction (Varotsos et al. , 2009. In this computation, only the earthquakes that occur in the area of the forthcoming earthquake, which is determined on the basis of SES properties (Varotsos and Alexopoulos 1984a, b) are considered. At the time of the mainshock, j 1 tends to 0 and for this reason, j 1 was proposed as an order parameter for seismicity (Varotsos et al. 2011c;Sarlis et al. 2008Sarlis et al. , 2010bVarotsos et al. 2005Varotsos et al. , 2010Varotsos et al. , 2011aVarotsos et al. , 2013Sarlis et al. 2013;Flores-Márquez et al. 2014).
Additionally, we examine the validity of Båth law for the events in the Global Centroid Moment Tensor (CMT) catalogue (Dziewonski et al. 1981;Ekström et al. 2012) for the time period 1/1/1976-31/12/2013 (available online on http://www.ldeo.columbia.edu/*gcmt/projects/CMT/cata log/jan76_dec13.ndk). The seismic moment M 0 which is available in the Global CMT catalogue is used (In detail, we utilize the first two columns of the fourth line for each of the catalogue's records for the moment exponent and the columns 50-56 of every fifth line for the scalar seismic moment mantissa. For example, if the exponent has the value 24, the moment values that follow, expressed in dyn cm, are multiplied by 10 24 ). Two methods are used for the selection of the aftershocks as described in Sect. 3 and both of them show a behaviour which is compatible with the model proposed in Sect. 2. The latter discussion is made in Sect. 4. Finally, the conclusions are presented in Sect. 5.

A model for the explanation of Båth law
Here, we propose a simple model in order to understand the origin of Båth law by making use of the fact that, as already stated in Sect. 1, the variance of natural time designated as j 1 may serve as an order parameter for seismicity. In addition, we consider that once N events have been observed, the ratio Q k = P N n¼1 Q n ð¼ p k Þ can be considered as a probability for the occurrence (observation) of the k-th event, as it has been discussed in detail in Varotsos et al. (2016) in the light of authoritative aspects of Max Planck as far as the meaning of the energy is concerned. We note that upon the occurrence of an additional event, the value of v k changes from k=N to k=ðN þ 1Þ together with the change of p k from Q k = P N n¼1 Q n to Q k = P Nþ1 n¼1 Q n thus leading to a change of j 1 as well.
We employ the model, depicted in Fig. 1, in which most of the events in the sequence foreshocks-mainshockaftershocks have very small energies. Thus, the probabilities p k of occurring are almost zero apart from those in the positions k 1 ¼ l þ 1 considered as the mainshock with probability p 1 and k 2 ¼ l þ m þ 2 considered as the largest magnitude aftershock detected so far with probability p 2 ¼ 1 À p 1 so as p 1 þ p 2 ¼ 1.
By inserting in Eq. (4) all p k ¼ 0 apart from those of the v k 1 ¼ ðl þ 1Þ=N and v k 2 ¼ ðl þ m þ 2Þ=N events with probabilities p 1 , p 2 , respectively, we obtain where N is the length of the sequence, and ðm þ 1Þ=N ¼ ½ðl þ m þ 2Þ À ðl þ 1Þ=N ¼ Dv is the distance of the two prominent events in natural time, so that We now set p 1 ¼ Q 1 = P 2 n¼1 Q n and p 2 ¼ Q 2 = P 2 n¼1 Q n and Q n , which is proportional to the energy the earthquake produced (Kanamori 1978), is given by Fig. 1 The model proposed in order to understand the origin of Båth law using the order parameter of seismicity j 1 . Time is set to 0 at the initiation of the SES activity where M n is the event's moment magnitude. The following relation is finally obtained: Figure 2 depicts in three dimensions j 1 in relation to the difference in magnitudes DM of the two events mainshocklargest aftershock detected so far and their difference in natural time Dv. Setting j 1 to its critical value, i.e. j 1 ¼ 0:07, we plot Dv versus DM in Fig. 3. We observe that j 1 may reach the value of 0:07, and thus a largest aftershock than the one already detected is expected to occur, when jDMj DM 1 and for Dv values larger than 0:5. The value of DM 1 ,i.e. the maximum allowed jDMj in order to expect a stronger earthquake, as shown in Fig. 3 can be found from Eq. (8) for Dv ¼ 1 and j 1 ¼ 0:07 which results in DM 1 % 0:724. Moreover, Fig. 2 shows that j 1 is a monotonic function of jDMj, a fact which means that the condition jDMj 0:724 is unique. The mathematical proof is provided in Appendix 1. Now we come to the possible explanation this model provides for the origin of Båth law. If we consider the first event as the mainshock, as already stated, the result of the previous paragraph shows that the second event (largest magnitude aftershock) has to have DM [ 0:724 so that the fulfilment of the condition j 1 ¼ 0:07 is impossible. Thus, the picture in natural time is compatible with a mainshocklargest aftershock pair and no stronger earthquake is expected to take place ''close'' to the mainshock. On the other hand, if DM 0:724 there is a probability that j 1 may reach the value 0.07 and a stronger earthquake will occur. In this case, the original picture mainshock-largest aftershock pair is falsified by the true events and hence this pair of events is not included in the calculation of hDMi of Eq. (1). Of course, there is also a possibility that j 1 may not reach 0.07 as long as the activity ''close'' to the mainshock is in progress but at a later time, thus this pair of events contributes to the aforementioned average value with DM 0:724. Such cases, however, are expected to be fewer than those with DM [ 0:724 since in the latter it is impossible to have a stronger earthquake than the two events detected so far. Thus, we expect that the average value of DM should be larger than 0.724 as suggested by Båth law.
Specifically, if we assume a piecewise constant probability distribution for DM of the form where HðxÞ is the Heaviside unit step function, i.e. HðxÞ ¼ 1 for x ! 0 and 0 otherwise, which is drawn with the red solid line in Fig. 4, we practically assume that the probability to observe DM [ DM 1 is k times larger than that of observing DM DM 1 . Using Eq. (9), we obtain (see Appendix 2) that the average value of DM is given by which for k ! 1 results in hDMi % 1:086, i.e. a value very close to that of Båth law in Eq. (1). We note that Zaliapin and Ben-Zion (Zaliapin and Ben-Zion 2013a) have observed that for earthquake sequences in southern California DM follows an almost uniform distribution in the range [0,2] with an exponential tail for larger values of DM, thus supporting the value hDMi = 1.1 which has been found (Zaliapin and Ben-Zion 2013a) in the observed seismicity of southern California. They have also identified (Zaliapin and Ben-Zion 2013b) that there exist two basic types of earthquake clustering: the burst-like sequences which are consistent with highly brittle behaviour and the swarm-like sequences consistent with mixed brittle-ductile behaviour. The latter sequences result in smaller DM values than those of the former. The suggested model is not inconsistent with this finding since it differentiates between small (DM\DM 1 ) and large DM classes of aftershock sequences.
If we further extend our model to allow DM to be larger than 2DM 1 , which is closer to the real case as we will see in the next Section, we may consider the distribution which is drawn with the black-dashed line in Fig. 4.
If we substitute in the latter expression values of l which are close to the experimental ones, we can obtain a value for Båth law close to 1.2 as in Eq. (1). This will be shown in Sect. 4.

Båth law in the Global CMT catalogue: the procedures followed and the results
For the calculation of DM, we adopt as earthquake magnitude the moment magnitude M w which results from the seismic moment M o through the relation (Kanamori 1978) where M 0 is measured in dyn Á cm (10 À7 N Á m). The minimal reported magnitude results in 4.6 (cf. a single 4.3 earthquake is reported in the Global CMT catalogue). The completeness magnitude threshold M c lies within 5.3-5.4 for the period before 2004 and decreases to 5.0 during the period 2004-2010 (Ekström et al. 2012). Two methods for the identification of the mainshock and its largest aftershock were used.

First method
In this method, we follow the definitions suggested in Shcherbakov et al. (2004): For each mainshock, we consider as aftershocks all earthquakes that took place within a space-time window. The time span of this window is 1 year (Shcherbakov et al. 2013;Zakharova et al. 2013) after the mainshock and the spatial dimensions are those of a rectangular region of side ''length'' of dðL=L 1 Þ degrees centred around the mainshock. We consider (Shcherbakov et al. 2004) and L 1 ¼ 111 km is the length per degree. This value of L 1 corresponds to the length of one degree either in the North-South or East-West direction at the equator. For each earthquake in the catalogue a space-time window as described above was created around it. If this earthquake was not the largest event inside its respective window, then it was rejected as a mainshock, otherwise it was considered as a mainshock. The largest aftershock of each sequence was then designated as the second maximum magnitude inside each mainshock's window and a total of 3100 pairs of mainshock-largest aftershock were identified.
We now compute the mean value hDMi ¼ hM ms À M max as i and the standard deviation r of DM for various mainshock magnitude thresholds, i.e. first we compute r and hDMi for M ms ! 5ð M ms;thres Þ then for M ms ! 5:1 and so on. Figure 5 depicts hDMi ¼ hM ms À M max as i and r versus the M ms;thres and the number of mainshocks. We find that for M ms;thres ! 5:7 the value of hDMi does not differ significantly from 1.2 (see the yellow region in Fig. 5 which has been drawn as a guide to the eye in view of the large observed r values) at least for the cases where the statistics come from more than 30 mainshocks, i.e. for M ms;thres up to 7.9. The linear trend that appears is certainly related to the relatively high completeness magnitude. We verified this fact by discarding all earthquakes with M\M c ¼ 5:3 and repeating the computations. By comparing the new hDMi curve with the existing one, we observed almost the same linear trend and a horizontal displacement. The two curves result in the same values for Fig. 4 The piecewise constant probability distribution assumed for DM (red solid line for Eq. (9) and black-dashed line for Eq. (11)) where the probability to observe DM [ DM 1 is larger than that of observing DM DM 1 M ms;thres ! 7:9. For M ms;thres around 8.3 and above (where the number of events is less than 10) the results are not credible because only a few events with large magnitudes have occurred (see Fig. 5). . In order to examine whether the DM distributions corresponding to different M ms;thres can be statistically discriminated in practice, we generated 10 3 samples of five DM values coming from the distribution for M ms;thres ¼ 5:7 (assumed Gaussian) and compared them to another five DM values coming from the distribution for M ms;thres ¼ 7:7 (also assumed Gaussian) by using the t test for the difference of the means of the two samples. This study revealed that there is on average a 15% probability that the two samples come from distributions with the same mean. This reveals that although hDMi may differ, the corresponding distributions do not depend on the mainshock's magnitude threshold in a statistically significant manner.

Second method
In this case, we consider the definition of aftershocks suggested in Zaliapin et al. (2008) which uses the nearestneighbour distance: For each earthquake i in a catalogue, we consider its occurrence time t i , latitude u i (in radians), longitude h i (in radians) and magnitude m i . For a pair of earthquakes i and j, the time-space-magnitude distance is defined (Baiesi and Paczuski 2004) as where t ij ¼ t j À t i is the interoccurrence time in years, d is the fractal dimension of the epicentres, b the Gutenberg-Richter parameter and r ij ! 0 is the spatial distance between the epicentres (Batac and Kantz 2014; Baiesi and Paczuski 2004): where R E ¼ 6371 km, the Earth's radius. We studied the Global CMT catalogue following the method suggested in Zaliapin and Ben-Zion (2013a, b, 2015) for the identification of earthquake clusters: For each earthquake j, going backwards in time, we find its nearest-neighbour i and their corresponding distance n ij . The nearest-neighbour, which is called the parent, can be the parent of multiple events which are called offsprings. We should note that each possible parent-offspring pair is assigned a corresponding distance n but the nearestneighbour is defined as the pair with minimum n ij (Zaliapin et al. 2008;Zaliapin and Ben-Zion 2013a). Normalizing the time and space distances by the magnitude of the parent, we have (Zaliapin et al. 2008;Zaliapin and Ben-Zion 2015) where T is referred to (Zaliapin et al. 2008;Zaliapin and Ben-Zion 2015) as the rescaled time and R as the rescaled distance. We selected b ¼ 1, d ¼ 1:6 and p ¼ 0:5 (Zaliapin et al. 2008;Zaliapin and Ben-Zion 2015).
As it was shown in Zaliapin and Ben-Zion (2013a), Zaliapin et al. (2008) and Zaliapin and Ben-Zion (2015), earthquake catalogues show a bimodal joint distribution of ðlog 10 T; log 10 RÞ. One mode corresponds to the background events-like events in a Poisson process-while  the other corresponds to the clustered events with n ij \10 À5 , i.e. events closer in time and space to their parents than expected in a Poisson process Ben-Zion 2013a, 2015;Zaliapin et al. 2008). Figure 6 depicts this bimodal distribution for the Global CMT catalogue. By discarding all n ij [ 10 À5 , the events in a catalogue are divided into single events and families that contain multiple events connected to their parents with nearest-neighbour distances n ij \10 À5 , i.e. strong links (Zaliapin and Ben-Zion 2013a). The largest event in a family is considered to be the mainshock. The events before it are considered as foreshocks, while the events after it are considered to be the aftershocks (Zaliapin and Ben-Zion 2013a). Analysing Global CMT catalogue by employing this technique, 3144 pairs of mainshock-largest aftershock were extracted, i.e. 3144 families.
In order to study Båth law, we keep our computations to these families and Fig. 7 depicts hDMi ¼ hM ms À M max as i and r versus the mainshock threshold M ms,thres and the number of mainshocks used in the calculation. Table 2 shows the corresponding values of hDMi and r for different M ms;thres values alongside their standard errors (Sarlis and Christopoulos 2012;Ahn and Fessler 2003).
Here, we also find that for M ms;thres ! 5:8 the value of hDMi does not differ significantly from 1.2 at least for the cases where the statistics come from more than 10 mainshocks, i.e. for M ms;thres up to 8.2 (see the yellow region in Fig. 7). There is again the linear trend but we have to note that the results of this method lead to the existence of a clear plateau in hDMi (see Fig. 7 for M ms;thres 2 ½6:8; 8:1) where hDMi varies with a standard deviation 0.06 around the value 1.29, very close to that of the Båth law of Eq. (1) (cf. a similar calculation for the first method for M ms;thres 2 ½7:1; 8:2 leads to 1.53 with standard deviation 0.07). In our opinion this result is not fortuitous, but is due to the fact that in the second method the aftershocks are selected according to the clustering properties of seismicity, e.g. Eq. (15) and Fig. 6. Rescaled time, log 10 (T ) Fig. 6 The joint probability density of the rescaled components ðlog 10 T; log 10 RÞ of the earthquake nearest-neighbour distance for the Global CMT catalogue. The darker colour corresponds to larger density. The upper right mode is the background mode while the left bottom mode is the clustered mode (see text).The red solid line log 10 T þ log 10 R ¼ À5 is also drawn   (Båth 1965). As an example, we consider M ms ! 7:5 which is also larger than or equal to 3DM 1 þ M c for the whole period of study (cf. M c % 5:35, see Sect. 3). In Figs. 8 and 9, we plot the cumulative distribution function (CDF) of DM for the first and the second method of aftershock identification, respectively, and model these distributions according to Eq. (11). For this reason, we perform two linear fits for each CDF by dividing the data into two intervals, one before DM 1 and the other after DM 1 covering the remaining (''linear'') portion of the CDF up to approximately 90% so as to avoid the effect of the tails. The intervals we selected are namely ½0; DM 1 Þ and ðDM 1 ; 3DM 1 and the corresponding results are shown in Figs. 8 and 9. The ratio between the two slopes is l ¼ 3:08 AE 0:19 for the first method of aftershock identification and l ¼ 2:01 AE 0:06 for the second method. By inserting these values for l in Eq. (12), we obtain 1.30 and 1.23, respectively, for the value in Båth law. Thus, by inserting the experimental values into the proposed model the resulting hDMi values are in accordance with the corresponding value in Båth law.

Conclusions
We examined the validity of Båth law in the Global CMT catalogue. We found that the corresponding value in Båth law does not markedly differ from 1.2 and the corresponding distributions of DM do not depend on the mainshock's magnitude threshold in a statistically significant manner. We note that by using the second method of aftershock identification, which makes use of the notion of earthquake families, the experimental values of hDMi we derived are closer to 1.2 than those obtained by the first method that uses a space-time window around each mainshock, e.g. see Tables 1 and 2. A simple model, utilizing the order parameter j 1 of seismicity in natural time, to explain the value of 1.2 has been proposed. Briefly, we use the fact that after the initiation of the SES activity and when studying the seismicity in the area prone to suffer the large earthquake, j 1 approaches the value 0.07 before the mainshock. Considering that this condition should be valid before every mainshock the model suggests that the probability to observe DM [ 0:724 should be markedly larger than that of DM 0:724. Based on this, two simple probability distributions for DM have been considered and the related expressions for hDMi have been obtained. When these expressions are compared to the experimental data and combined with the modern aspects of earthquake clustering they may lead to values of hDMi close to 1.2 in accordance with the Båth law. where a [ 0. Since the model does not take under consideration whether the mainshock is larger than its largest aftershock or not, there are two possibilities: (1) DM [ 0 so that sinhð1:5DMln10Þ [ 0. In this case, oj 1 =oDM\0.
(2) DM\0 so that sinhð1:5DMln10Þ\0 since sinhðÀxÞ ¼ À sinhðxÞ. In this case, we have oj 1 =oDM [ 0. In either of the cases above, which are summarized in Fig. 2, j 1 is a monotonic function of jDMj, a fact which means that the condition jDMj 0:724 is unique.

Appendix 2
By normalizing the probability of observing DM, i.e. demanding the area under the red solid line in Fig. 4 to be unity, we have p 0 DM 1 þ kp 0 DM 1 ¼ 1 ð22Þ and solving for p 0 we derive The average value of DM is given by and by substituting Eq. (23) for p 0 we finally obtain For k ! 1, the last equation gives which for DM 1 ¼ 0:724 results in hDMi ¼ 1:086 as stated in the main text.

Appendix 3
By normalizing the probability of observing DM, i.e. demanding the area under the black-dashed line in Fig. 4 to be unity, we have and solving for p 0 we find The average value of DM is given by and by substituting Eq. (28) for p 0 we finally obtain which is just Eq. (12) of the main text.