SUMMARY

The Epidemic Type Aftershock Sequence (ETAS) model describes how an earthquake generates its own aftershocks. The regular ETAS model assumes that distribution F of the number of direct aftershocks is Poissonian, however there is evidence suggesting that a geometric distribution might be more adequate. Let |${{V}_M}({{m}_ \bullet })$| be the number of |$m > M$| aftershocks generated by |${{m}_ \bullet }$|event. In this study we consider the |${{V}_M}({{m}_ \bullet })$| distribution within Epidemic-type Seismicity models, ETAS(F). These models include the Gutenberg–Richter law for magnitude and Utsu law for average |${{m}_ \bullet }$|-productivity, but differ in the type of F distribution for the number |$v({{m}_ \bullet })$| of direct aftershocks. The class of F is quite broad and includes both the Poisson distribution, which is the basis for the regular ETAS model, and its possible alternative, the Geometric distribution. We replace the traditional |$M = {{m}_ \bullet } - \Delta $| threshold in |$\Delta $|-analysis with |$M = {{m}_a} - \Delta $| where |${{m}_a}$| is the distribution mode of the strongest aftershocks. Under these conditions we find the limit |${{V}_M}({{m}_ \bullet })$| distribution at |${{m}_ \bullet } > > 1$|⁠. In the subcritical case, the limit distribution is extremely simple and identical to the |$v({{m}_\Delta })$| distribution with a suitable magnitude |${{m}_\Delta }$|⁠. This result allows us to validate both the priority of the geometric distribution of F for direct aftershocks and the very concept of epidemic-type clustering on global seismicity data.

1 INTRODUCTION

We are interested in analysing the distribution |${{V}_M}({{m}_ \bullet })$| of the number of aftershocks with magnitude |$m \ge M$|⁠, caused by an |${{m}_ \bullet }$| earthquake. To weaken the dependence of the distribution on |${{m}_ \bullet }$|⁠, the so-called |$\Delta $|-analysis is traditionally used, assuming that M is equal to |${{M}_ \bullet } = {{m}_ \bullet } - \Delta \ge {{m}_c}$|⁠. Such approach, however, does not provide a clear answer about the |${{V}_{M \bullet }}({{m}_ \bullet })$| distribution to be preferred.

Soloviev & Solovieva (1962) were probably the first to postulate the Geometric distribution for |${{V}_{M \bullet }}({{m}_ \bullet })$|⁠. This conclusion was based on relatively scarce data (1954–1961) for the Pacific Belt (⁠|${{m}_ \bullet } \approx 7$|⁠; |$\Delta = 2$|⁠) and for the Kamchatka and Kuril Islands region (⁠|${{m}_ \bullet } \approx 6$|⁠; |$\Delta = 2$|⁠), Shebalin et al. (2018) confirmed the hypothesis using 850 aftershock sequences with main shocks |${{m}_ \bullet } \ge 6.5$| from the ANSS (1975–2018) global catalogue. Differently, Kagan (2010), using PDE (1977–2007) global data provided by USGS (United States Geological Survey), for |${{m}_ \bullet } = 7.1 - 7.2$| and |$\Delta = 2.5$|⁠, obtained a bimodal |${{V}_{M \bullet }}({{m}_ \bullet })$| distribution with a dominant peak at zero value of |${{V}_{M \bullet }}({{m}_ \bullet })$|⁠. Zaliapin & Ben-Zion (2016) used the global catalogue compiled at Northern California Earthquake Data Center, NCEDC (1975–2015, |${{m}_c} = 4$|⁠, |$\Delta = 2$|⁠, 2000 aftershock sequences) and clusters triggered not necessarily by the main (strongest) event. It was shown that the survival function |$\bar{P}(n) = P({{V}_{M \bullet }}({{m}_ \bullet }) > n)$| looks like a power-law |$C{{n}^{ - 2}}$| for locations with high heat flow, and as a parabola in the |$(\textit{logn},log\bar{P})$| representation otherwise.

Molchan et al. (2022) used epidemic-type cluster models to theoretically analyse the problem and found conditions for the existence of heavy and light tails in the |${{V}_{M \bullet }}({{m}_ \bullet })$| distribution. It was also shown that the distribution of |${{V}_{M \bullet }}({{m}_ \bullet })$| and the distribution of the number |$v({{m}_ \bullet })$| of direct aftershocks cannot be of the same type, say Poisson or Geometric (the exact definition of the type will be given later). This statement is inconsistent with the statement of Shebalin et al. (2018, 2020) on the Geometric distribution for both |${{V}_{M \bullet }}({{m}_ \bullet })$| and |${{v}_\Delta }({{m}_ \bullet })$|⁠, where |${{v}_\Delta }({{m}_ \bullet })$| is the number of direct aftershocks with |$m \ge {{m}_ \bullet } - \Delta $|⁠.

The uncertainty about the |${{V}_{M \bullet }}({{m}_ \bullet })$| distribution requires specific analysis. If the magnitude gap |$\delta {{m}_ \bullet }$| between the main event |${{m}_ \bullet }$| and its strongest aftershock, for instance, grows with |${{m}_ \bullet }$|⁠, the real magnitude range, |$\Delta - \delta {{m}_ \bullet }$|⁠, for aftershocks counting in |${{V}_{M \bullet }}({{m}_ \bullet })$| statistics decreases. In such case the choice |$M = {{m}_ \bullet } - \Delta $| excludes the independence of the |${{V}_{M \bullet }}({{m}_ \bullet })$| distribution on |${{m}_ \bullet }$|⁠. In particular, |${{V}_{M \bullet }}({{m}_ \bullet })$| statistics can be abnormally small for large |${{m}_ \bullet }$|⁠, if the gap |$\delta {{m}_ \bullet } > > \Delta $|⁠. A more flexible way is to use |$M = {{m}_a} - \Delta := {{M}_a}$|⁠, where |${{m}_a}$| is, for example, the peak value in the distribution of the magnitude of strongest aftershock in the clusters. We used this idea to find the distribution of |${{V}_M}_a({{m}_ \bullet })$| at |${{m}_ \bullet } > > 1$| in the framework of epidemic-type seismicity models, ETAS(F). These models include the Gutenberg–Richter magnitude distribution and the exponential Utsu productivity law, |$\lambda (m) = E\nu (m)$|⁠. They differ in the type of F distribution for the number |$\nu (m)$| of direct aftershocks. The F-distribution set includes both the traditional Poisson distribution and its possible alternative, the Geometric distribution. The magnitude distributions of the strongest aftershock for this F set have been studied previously in in Molchan & Varini (2024).

We are going to show that at the modified choice of M, the limit |${{V}_M}_a({{m}_ \bullet })$| distribution exists and is determined by the distribution of direct aftershocks. Moreover, in the subcritical regime, the limit |${{V}_M}_a({{m}_ \bullet })$| distribution is identical to the |$\nu ({{m}_\Delta })$| distribution with a suitable initial magnitude |${{m}_\Delta }$|⁠. This result agrees well with reality and corrects the hypothesis of Soloviev & Solovieva (1962) and Shebalin et al. (2018) about the geometric |${{V}_M}_ \bullet ({{m}_ \bullet })$| distribution type.

The basic concepts and seismostatistical results of the paper are contained in Sections 2.1, 2.2, 4 and 5. The remaining sections contain mathematical statements and proofs.

2 PRELIMINARY: KEY CONCEPTS AND FACTS

This paper is a continuation of Molchan et al. (2022) and Molchan & Varini (2024). Some key points and facts from these papers are summarized below.

2.1 ETAS(F) model

The regular epidemic-type aftershock sequence model (ETAS) is defined by means of the conditional intensity of an event at a point of the phase space S = (time, location, magnitude), given the past of the seismic process up to the current time (Hawkes & Oakes 1974; Ogata 1988). The model is not flexible enough to describe direct aftershocks (Zaliapin & Ben-Zion. 2013; Shebalin et al 2020). This shortcoming can be easily overcome by describing a cluster of seismic events in S as a Galton-Watson tree; in this way, the coordinates (time, location) become irrelevant. Accordingly, we recall the model in projection to the magnitude.

The initial event of magnitude m generates a random number |$\nu (m)$| of direct aftershocks according to the off-spring law F:

(1)

Each of these events is assigned an independent magnitude in accordance with the law F1

(2)

(Here and hereinafter the threshold magnitude |${{m}_c}$| is taken as a reference point, |${{m}_c} = 0$|⁠).

Each event of the first generation independently generates new ones, following the described rule for its parent, etc. The process continues indefinitely. It stops with probability 1 if the criticality index (branching ratio)

(3)

where |$\lambda (m) = E\nu (m)$| is the average productivity of the m-event. The values |$n < 1$| and |$n = 1$| correspond to the subcritical and critical regimes, respectively.

We get the standard ETAS model if F has a Poisson distribution (P), that is

(4)

Shebalin et al. (2020) found that the Geometric distribution (G),

(5)

agrees better with empirical data.

2.2 Random thinning (RT) property

In the ETAS(F) model, F is a family of distributions parametrized by the magnitude m of the initial cluster event or its productivity |$\lambda (m)$|⁠, shortly |$F = \{ {{F}_\lambda }\} $|⁠, where |$\lambda $| is the mean value. Below we will consider a special class of F distributions that have the so-called Random Thinning (RT) property. To explain this, consider a random number |${{\nu }^{(\lambda )}} \ge 0$| of events with mean value |$\lambda $|⁠. Let |${{R}_p}\nu _{}^{(\lambda )}$| be the number of subevents that remain after randomly choosing each of the events with probability p, that is |${{R}_p}\nu _{}^{(\lambda )} = {{\varepsilon }_1} + ... + {{\varepsilon }_{{{\nu }^{(\lambda )}}}}$|⁠, where |${{\varepsilon }_i} \in \{ 0,1\} $| are independent and |$P({{\varepsilon }_i} = 1) = p$|⁠. For example, if the magnitudes of direct aftershocks are independent, the number of direct aftershocks |${{\nu }_\Delta }({{m}_ \bullet })$| associated with |$m > {{m}_ \bullet } - \Delta $| is stochastically equivalent to a random variable resulting from the thinning operation of |$\nu ({{m}_ \bullet })$| with probability |$p = 1 - {{F}_1}({{m}_ \bullet } - \Delta )$|⁠.

We will say that the family of random variables |${{\nu }^{(\lambda )}},\lambda \ge 0$| has RT property if |${{\nu }^{(\lambda p)}}$| and |${{R}_p}\nu _{}^{(\lambda )}$| are stochastically identical for any fixed |$\lambda > 0$| and 0 < p < 1:

that is both random variables have the same distributions. In this case, |${{\nu }^{(\lambda )}}$| and |${{R}_p}\nu _{}^{(\lambda )}$| differ only in their mean values, but do not change their type, which corresponds to the reference |${{\nu }^{(\lambda )}}$| distribution with mean |$\lambda = 1$|⁠.

Since |$E{{z}^{{{R}_p}{{\nu }^{(\lambda )}}}} = E{{(1 - p + pz)}^{{{\nu }^{(\lambda )}}}},$| the generating function of |${{\nu }^{(\lambda )}}$| with RT property has the following representation:

(6)

The |$\varphi (w),w < 0$| function uniquely determines the common distribution type for the family of random variables |${{\nu }^{(\lambda )}}$|⁠. Therefore, we will use the |${{\varphi }_F}$| symbol to denote the type of F-distribution. The class of random variables with the RT property is quite wide. The following statement shows that all RT distributions are limit points of distributions of |${{R}_{1/\lambda }}{{\nu }^{(\lambda )}}$| variables for which |$E{{\nu }^{(\lambda )}} = \lambda \gg 1$|⁠.

Statement 1. Let |$\{ {{\phi }_n}(z)\} $| be a sequence of generating functions of |${{\nu }^{({{\lambda }_n})}}$| with means |${{\lambda }_n} \to \infty $|⁠. Then there exists a subsequence |$\{ \tilde{n}\} $| such that for any |$\lambda > 0$| there exists the limit

where |$\varphi (z - 1)$| is the generating function of some distribution with unit mean.

The proof and properties of |$\varphi (w)$| are described in Appendix.

Examples. Let us note important examples of distributions with the RT property:

  1. The Negative Binomial distribution|$F = NB(\tau )$|⁠. This |$\tau - $|family of distributions includes the following RT types:
    (7)

The case |$\tau = 1$| corresponds to the Geometric distribution, while the limiting case |$\tau = \infty $| corresponds to the Poisson distribution. Decreasing |$NB(\tau )$| distributions were used in Kagan (2010; 2014) and Zaliapin & Ben-Zion (2013) to describe the number |$\nu (m)$| of direct aftershocks.

Among all |$NB(\tau )$| types, the geometric and Poisson distributions are distinguished by the following features: only in the Poisson case |$(\tau = \infty )$| the variance, |${{\sigma }^2}$|⁠, is minimal for any fixed mean |$\lambda > 0$|⁠; the same is true for the geometric distribution |$(\tau = 1)$| but with the additional condition (a) : the distribution density is decreasing function at any mean |$\lambda > 0$|⁠.

This follows from two facts: (i) |${{\sigma }^2} = {{\lambda }^2}(1/\tau + 1/\lambda )$| and (ii) condition (a) holds if and only if the coefficient of variation |$\sigma /\lambda \ge 1$| for any |$\lambda > 0$|⁠.

  1. The Kovchegov & Zaliapin (2021) distribution having the following RT type:
    (8)

Seismo-statistical applications of such distributions and their invariant properties are discussed in Kovchegov et al. (2022, 2023).

The distributions related to RT types eqs (7) and (8) have different tails: light |${{p}_n} \approx c{{n}^{\tau - 1}}{{(1 + \tau /\lambda )}^{ - n}},n > > 1$| in eq. (7), and heavy |${{p}_n} \approx c{{n}^{ - \tau - 1}},n > > 1$| in eq. (8).

2.3 Cluster thinning

In applications, the ETAS(F) model depends on a lower magnitude threshold |${{m}_0}$| which we have taken as a reference point. Natural candidates for |${{m}_0}$| may be a |${{m}_c}$| threshold for the representativeness of magnitudes or a fertility threshold, |${{m}_f}$|⁠, that determines the ability of the m event to generate aftershocks (Li et al. 2024).

The RT property for F distribution allows us to establish a simple relationship between models with different values of |${{m}_0}$|⁠.

Within an ETAS(F) cluster with |${{m}_0} = 0$|⁠, we can identify a subcluster of causally related events with |$m \ge {{m}_f} > 0$|⁠, starting from the initial event. Following Duquesne & Winkel (2007), we remove from the original cluster all events with |$m \le {{m}_f}$| along with their descendants of all generations. Then, we shift the magnitude to use the scale |$\tilde{m} = m - {{m}_f}$|⁠. The resulting sub cluster we will call ETAS(F, |${{m}_f}$|⁠). The following statement is a simple consequence of the RT property of F.

Statement 2. If the ETAS(F) model is such that F distribution has the RT property, and

then the ETAS(F, |${{m}_f}$|⁠) model is stochastically equivalent to the ETAS|$(\tilde{F})$|⁠, where |$\tilde{F}$| has the same RT type as F,

In the regular case: |$\lambda (m) = {{\lambda }_0}{{e}^{\alpha m}},{{\bar{F}}_1} = {{e}^{ - \beta m}}$|⁠, the corresponding ETAS|$(\tilde{F})$| characteristics look similar except for |${{\lambda }_0}$|⁠, which needs to be replaced with |${{\tilde{\lambda }}_0} = {{\lambda }_0}{{e}^{ - (\beta - \alpha ){{m}_f}}}$|⁠. The RT type of F distribution will remain unchanged.

As we can see, cluster modelling with different lower magnitude thresholds does not take us beyond the chosen ETAS(F) model with fixed RT type of F.

2.4 Cluster types

The initial cluster magnitude |${{m}_ \bullet }$| may be the strongest (dominant) or arbitrary. Therefore, we will consider clusters of both types: |$DM({{m}_ \bullet })$| with a Dominant initial Magnitude and |$AM({{m}_ \bullet })$| with an Arbitrary initial Magnitude. Statement 6 in the Appendix show that any |$DM({{m}_ \bullet })$| cluster in ETAS(F) model can be viewed as an |$AM({{m}_ \bullet })$| cluster in the ETAS|$(\hat{F})$| model with a suitable off-spring distribution |$\hat{F}$|⁠.

2.5 The strongest cluster event

The following statement describes the limit distributions of strongest aftershock magnitude in a cluster caused by large event |${{m}_ \bullet } > > 1$|⁠.

Statement 3 (Molchan & Varini 2024).

Consider the ETAS(F) model with the distribution F having the RT property, that is

Suppose that

(9)
(10)
(11)

where |$\alpha \le \beta $|⁠, |$1 \le {{K}_1} < \infty $|⁠; in addition, |${{K}_1} \le \infty $| if |$\alpha < \beta $|⁠.

Let |${{\mu }_a}$| be a maximal magnitude in the AM/DM cluster and n is the criticality index.

A. Under |${{m}_ \bullet } > > 1$| condition, the following regression relationships are valid for AM/DM clusters:

(12)

where

(13)
(14)

[since n is fixed, |${{\lambda }_0}$| in eq. (14) is a function of |${{m}_ \bullet }$|⁠:|${{\lambda }_0} = n \cdot {{(\beta {{K}_1})}^{ - 1}}/{{m}_ \bullet }$|];

(15)

and |$A = 2{{(\beta {{\varphi ^{\prime\prime}}_F}(0))}^{ - 1}}(\beta - 2\alpha + {{1}_{\beta = 2\alpha }})$|⁠.

The random component |$\varsigma $| in eq. (12) has the limit distribution:

(16)

B. In the case |$(\alpha < \beta < 2\alpha ,n = 1)$|⁠, the regression relationship (eq. 12) depends on the cluster type.

Let |$\psi (x),x > 0$| be defined by the following equation

(17)

Then regression (eq. 12) is valid for the AM- cluster with the parameters:

(18)

In the case of a DM cluster, the limit |${{\mu }_a}$| distribution is as follows

(19)

For small |${{\lambda }_0} \approx 1 - \alpha /\beta $|⁠, eq. (19) admits the following simplification:

Remark. The non-random component, |${{m}_a}$|⁠, in the regression relation (eq. 12), will be conventionally called the peak value of |${{\mu }_a}$|⁠. In the most interesting case, when F is the Negative-Binomial distribution |$NB(\tau )$|⁠, this conventional name becomes precise (Molchan & Varini 2024).

3 MAIN RESULTS

3.1 The main equation

From the structure of the ETAS(F) model, it follows that an m-event generates |$\{ {{\mu }_i},i = 1,...,\nu (m)\} $| direct aftershocks, which in turn independently generate clusters of events with magnitude m > M and sizes |$\{ {{V}_M}({{\mu }_i})\} $|⁠. All random variables |${{V}_M}({{\mu }_i})$| are independent and equally distributed, given the random nature of the initial event |${{\mu }_i}$|⁠. As a result we get the following stochastic equation:

(20)

where the |${{1}_{{{\mu }_i} > M}}$| term accounts for the |${{\mu }_i}$| aftershock if |${{\mu }_i} > M$|⁠.

3.2 Choice of M

If |$n < 1$|⁠, eq. (20) entails

(21)

For the exponential |$\lambda (m)$| and |${{f}_1}(m)$| characteristics presented in eqs (10) and (11),

Given |$M = {{m}_ \bullet } - \Delta $| we have |$E{{V}_M}({{m}_ \bullet }) \approx C{{e}^{ - (\beta - \alpha ){{m}_ \bullet }}}{{e}^{\beta \Delta }} = o(1)$|⁠. Hence, to have a non-trivial limit |${{V}_M}({{m}_ \bullet })$| distribution at |${{m}_ \bullet } > > 1$|⁠, the choice of M should be such that |$0 < c < E{{V}_M}({{m}_ \bullet }) < C,{{m}_ \bullet } \to \infty $|⁠, that is |$M = (\alpha /\beta ){{m}_ \bullet } + \textit{const}$|⁠. The non-random component |${{m}_a}$| of the strongest aftershock (see Statement 3) has exactly the same form in the subcritical regime.

The choice |$M = {{m}_a}$| is quite probable and in the general case. Indeed the events |$\{ {{V}_M}({{m}_ \bullet }) = 0\} $| and |$\{ {{\mu }_a} < M\} $| are identical for any M, and hence

(22)

Assuming |$M = {{m}_a}$| the right-hand part of eq. (22) will have a non-trivial limit at |${{m}_ \bullet } > > 1$|⁠.

Relation (22) shows that the problem under consideration includes the problem of distribution of the strongest aftershock.

3.3 The limit VMa(m) distribution

We will consider the ETAS(F) model described in Statement 3 and use the notation adopted there.

Statement 4. Under the assumptions of Statement 3, the (⁠|${{V}_M}_a({{m}_ \bullet })$|⁠,|${{M}_a} = {{m}_a} - \Delta $|⁠) statistics has a non-trivial limit distribution at |${{m}_ \bullet } > > 1$|⁠.

I. In |$(\alpha < \beta ,n < 1)$| and |$(\alpha = \beta ,n \le 1)$| regimes, the limit |${{V}_M}_a({{m}_ \bullet })$| distribution is independent of the AM/DM cluster type and has generation function |${{\varphi }_F}({{e}^{\beta \Delta }}(z - 1))$| corresponding to the |$\nu ({{m}_\Delta })$| statistics with mean |$\lambda ({{m}_\Delta }) = {{e}^{\beta \Delta }}$|⁠, that is

(23)

In the case |$(2\alpha \le \beta ,n = 1)$|⁠, the limit distribution is also independent of the AM/DM cluster type but has an infinite mean and generation function

(24)

II. The case |$(\alpha < \beta < 2\alpha ,n = 1)$|⁠, AM cluster

Let |$M = {{m}_ \bullet } - \Delta $|⁠, then the limit |${{V}_M}({{m}_ \bullet })$| distribution has a generation function |${{\varphi }_F}( - {{e}^{\alpha \Delta }}{{\tilde{A}}_M}(z))$|⁠, where |${{\tilde{A}}_M}(z)$| is defined by the equation

(25)

At small |${{\lambda }_0}$|⁠, that is |$\alpha /\beta \approx 1$|⁠,

(26)

III. The case |$(\alpha < \beta < 2\alpha ,n = 1)$|⁠, DM cluster.

If |$M = {{m}_ \bullet } - \Delta $|⁠, then the limit |${{V}_M}({{m}_ \bullet })$| distribution has generation function |${{\varphi }_F}({{\hat{A}}_M}(z))$|⁠, where |${{\hat{A}}_M}(z)$| is defined by the equation

(27)

|$\mu (dy) = \beta /\alpha {{y}^{ - \beta /\alpha - 1}}dy$| and |${{\lambda }_0} = (1 - \alpha /\beta )$|⁠.

At small |${{\lambda }_0}$|⁠, that is |$\alpha /\beta \approx 1$|⁠,

(28)

that is

Remark 1. The generating function |${{\varphi }_F}( - {{(1 - z)}^\gamma }{{\lambda }_\Delta }),0 < \gamma < 1$| in eqs (24) and (26) corresponds to the random variable |$V \sum\nolimits_{i = 1}^{\nu ({{m}_\Delta })} {{{\varsigma }_i}} $| with independent components |$\{ \nu ({{m}_\Delta }),{{\varsigma }_1},...{{\varsigma }_i},..\} $| such that |$E{{z}^{{{\varsigma }_i}}} = 1 - {{(1 - z)}^\gamma }$| and |$E\nu ({{m}_\Delta }) = {{\lambda }_\Delta }$|⁠.

Applying to V the general theory of subexponential distributions (e.g. Denisov et al. 2010), we get the following consequence:

if |$\nu ({{m}_\Delta })$| has all moments (it is the case of F = P/G), then for large n

(29)

where Г(x) is the Euler Gamma function.

For the regular ETAS model, Saichev et al. (2005) obtained a similar result for the total number of cluster events V. Namely, at the critical regime |$P(V = n) = O({{n}^{ - 1 - \gamma }})$|⁠, where |$\gamma = \max (1/2,\alpha /\beta )$| is the same as in eqs (24) and (26). Importantly, this result contains no constraints on the proximity of |$\alpha $| and |$\beta $|⁠.

Remark 2 Relationship to Statement 3. According to eq. (22), the function |$M \to P\{ {{V}_M}({{m}_ \bullet }) = 0\} $| determine the magnitude distribution of the strongest aftershock |${{\mu }_a}$|⁠: |$P\{ {{\mu }_a} < M\} = E{{z}^{{{V}_M}({{m}_ \bullet })}}| {_{z = 0}} $|⁠. In the limit case this connection requires additional verification. In regimes (n < 1) and |$(\alpha < \beta < 2\alpha ,n = 1)$|⁠, it does not cause difficulties because of the simplicity of the limit generating function of |${{V}_M}_a({{m}_ \bullet })$|⁠. The case |$(2\alpha > \beta ,n = 1)$| requires additional analytics, which are presented in the Appendix (see A39, A55–A57).

Remark 3 AM/DM clusters. The probabilistic structure of AM/DM clusters is similar and differs only in |$({{f}_1}(m),\lambda (m))$| characteristics that become close at |${{m}_ \bullet } > > 1$|⁠. Therefore, it is quite expected that the limit |${{V}_{{{M}_a}}}({{m}_ \bullet })$| distributions for both cluster cases will be the same. These intuitive considerations proved themselves in the subcritical regime and even in the critical one, if |$\alpha = \beta $| or |$2\alpha \le \beta $|⁠. Radical differences in the |${{V}_{{{M}_a}}}({{m}_ \bullet })$| distributions for AM/DM clusters arise when (⁠|$2\alpha > \beta ,n = 1)$|⁠, that is when the variance of the number of direct aftershocks |$\nu (\mu )$| for random |$\mu $| is unbounded and the regime is critical.

Remark 4|${{V}_{M \bullet }}({{m}_ \bullet })$|distribution. In the ETAS(F) model, statistics |$\{ {{{V}_{M \bullet }}({{m}_ \bullet }),{{M}_ \bullet } = {{m}_ \bullet } - \Delta } \}$| and |$\nu ({{m}_ \bullet })$| cannot have the same RT type (Molchan et al. 2022). Replacing threshold |$M = {{m}_ \bullet } - \Delta $| with |${{M}_a} = {{m}_a} - \Delta $| corrects the situation, at least in part. According to (23), in subcritical regime, the limit |${{V}_M}_a({{m}_ \bullet })$| distribution coincides with the |$\nu ({{m}_\Delta })$| distribution having the following mean and variance:

(30)

where |${{\varphi ^{\prime\prime}}_F}(0) = 2$| and 1 for Geometric, F = G and Poisson, F = P, models respectively.

The existence of a non-trivial limit for |${{V}_{M \bullet }}({{m}_ \bullet }),$| it is possible only if |${{m}_a} - {{m}_ \bullet } = C$| is a constant. This is only true in the critical regime under the condition |$(\alpha < \beta < 2\alpha )$|⁠. However, in the practically interesting case then |$\alpha /\beta \sim 1 $|⁠, the difference |${{m}_a} - {{m}_ \bullet }$| depends weakly on |${{m}_ \bullet }$| independently of the criticality index, |$n \le 1$|⁠. Therefore, in this case too, the |${{V}_{M \bullet }}({{m}_ \bullet })$| distributions may be close to the |$\nu ({{m}_{\tilde{\Delta }}})$| distribution with |$\tilde{\Delta } = \Delta - C$|.

4 TESTING OF VM (m) DISTRIBUTIONS

Comparison of the limit |${{V}_{{{M}_a}}}({{m}_ \bullet })$| distribution with reality provides an opportunity to test the type of F-distribution, as well as the very concept of epidemic-type clustering.

4.1 Main events and clusters identification

Since the limit distribution is associated with strong |${{m}_ \bullet }$| events, it is difficult to find a large enough homogeneous sample of clusters for testing purposes. Following Zaliapin & Ben-Zion (2016), we consider the main events |${{m}_ \bullet } \ge 6$| reported in the ANSS (2024) global catalogue for the time period 1.1.1980–1.4.2024, using the lower magnitude threshold Mc = 4.0 and depth |$\le $|70 km.

We identify earthquake clusters using the nearest neighbour (NN) method (Zaliapin & Ben-Zion 2016). The parameters required by the NN method to compute the generalized distance |$\eta $| between two earthquakes, namely the slope of Gutenberg-Richter law (b-value) and the fractal dimension of epicentres (df), are set for the whole global catalogue to the following values: b-value = 1.0, df = 1.3. The threshold distance |${{\eta }_0}$|⁠, used to separate clusters from background seismicity, is set automatically by the Gaussian mixture model approach (see Zaliapin & Ben-Zion 2016 for further details). This resulted in the identification of 1427 clusters with |${{m}_ \bullet } \ge 6$| and 189 clusters with |${{m}_ \bullet } \ge 7$|⁠.

4.2 Strongest aftershocks |${{\mu }_a}({{m}_ \bullet })$| and their ‘peak values’ |${{m}_a}({{m}_ \bullet })$|

The |${{V}_M}({{m}_ \bullet })$| statistics with |$M = {{m}_a} - \Delta $| threshold contains unknown |${{m}_a}$| value. According to Statement 3, |${{m}_a}$| is a non-random component in the linear regression |${{\mu }_a}$| versus |${{m}_ \bullet }$|⁠, that is

(31)

The corresponding |${{m}_a} = A{{m}_ \bullet } + B$| estimates are presented in Fig. 1. In the interval |$6 \le {{m}_ \bullet } < 6.5$| the |${{\mu }_a}$| data are censored from below by the threshold |${{m}_c} = 4$|⁠. Therefore, we display two versions of the linear regression:

(32)

In the ETAS(F) model |$A = \alpha /\beta \le 1$|⁠, therefore the A parameter in eq. (32) is weakly defined.

Relationship ${{\mu }_a}$ versus ${{m}_ \bullet }$ for ${{m}_ \bullet } \ge 6$ and linear regressions ${{m}_a} = A{{m}_ \bullet } + B$ based on: Båth's law ${{m}_a} = {{m}_ \bullet } - 1.2$ (red line); ${{m}_ \bullet } \ge 6.5$ data, A = 1.03, B = −1.51 (dashed line); ${{m}_ \bullet } \ge 6$ data, A = 0.84, B = −0.13 (dotted line).
Figure 1.

Relationship |${{\mu }_a}$| versus |${{m}_ \bullet }$| for |${{m}_ \bullet } \ge 6$| and linear regressions |${{m}_a} = A{{m}_ \bullet } + B$| based on: Båth's law |${{m}_a} = {{m}_ \bullet } - 1.2$| (red line); |${{m}_ \bullet } \ge 6.5$| data, A = 1.03, B = −1.51 (dashed line); |${{m}_ \bullet } \ge 6$| data, A = 0.84, B = −0.13 (dotted line).

Another way of estimating |${{m}_a}$| is based on a version of Båth's law, namely, the assumption that the distribution of the Båth's gap, |$\delta = {{m}_ \bullet } - {{\mu }_a}({{m}_ \bullet })$|⁠, is independent of |${{m}_ \bullet }$|⁠. According to Zaliapin & Ben-Zion (2016) this assumption implies homogeneity of the data. It follows from the assumption that

(33)

where |$\hat{\delta }$| is the mode of the |$\delta $|-distribution.

Fig. 2 shows |$\delta $| distributions for |${{m}_ \bullet } \ge 6$| and |${{m}_ \bullet } \ge 7$|⁠; the distributions are similar and have common classical mode |$\hat{\delta } \approx 1.2$|⁠. This justifies the use of the estimate (eq. 33) suitable in the regime |$\alpha /\beta \sim 1 $| For completeness of the aftershock statistics, we use the main shocks |${{m}_ \bullet }$| for which |${{m}_a}({{m}_ \bullet }) - \Delta > {{m}_c}$|⁠.

Distribution of Båth's gap, $\delta = {{m}_ \bullet } - {{\mu }_a}({{m}_ \bullet })$ based on main shocks: ${{m}_ \bullet } \ge 6$ (blue) and ${{m}_ \bullet } \ge 7$ (red). The curves correspond to Gaussian approximations with parameters (mean, sigma): (1.18, 0.56) and (1.27, 0.59), respectively.
Figure 2.

Distribution of Båth's gap, |$\delta = {{m}_ \bullet } - {{\mu }_a}({{m}_ \bullet })$| based on main shocks: |${{m}_ \bullet } \ge 6$| (blue) and |${{m}_ \bullet } \ge 7$| (red). The curves correspond to Gaussian approximations with parameters (mean, sigma): (1.18, 0.56) and (1.27, 0.59), respectively.

4.3 Main characteristics of the limit |${{V}_{{{M}_a}}}({{m}_ \bullet })$| distribution

By virtue of eq. (30), the interval |$I(\Delta )$|⁠: |$\lambda ({{m}_\Delta }) \pm 2\sigma (\Delta )$| contains the principal values of the limit |${{V}_{{{M}_a}}}({{m}_ \bullet })$| statistics. Therefore, reasonable intervals for testing the |${{V}_{{{M}_a}}}({{m}_ \bullet })$| distribution are

(34)

For a geometric distribution with mean |${{e}^{\beta \Delta }}$|⁠, the value |$3{{e}^{\beta \Delta }}$| corresponds to approximately 95 per cent quantile. The exact 95 per cent interval for F = G is

(35)

Graphical representation of the|${{V}_{{{M}_a}}}({{m}_ \bullet })$|distribution. In frame of the ETAS(F) model, the limit |${{V}_{{{M}_a}}}({{m}_ \bullet })$| distribution preserves the F type. Considering F = G as the main hypothesis, it is natural to analyse the size distributions using the survival function |$\bar{P}(n)$|=Prob(cluster size > n) on the plane (X = n, Y = log|$\bar{P}$|⁠).

In this representation, the limit distribution is a linear function with slope

(36)

For F = P, the limiting distribution is Poissonian, which is close to Gaussian when the mean, |${{e}^{\beta \Delta }}$|⁠, is large. Therefore, the limit distribution on the plane X = n, Y = ln |$\bar{P}$| is concave with the vertex at the |${{e}^{\beta \Delta }}$| point.

Zaliapin & Ben-Zion (2016) consider the |${{V}_{{{m}_ \bullet } - 2}}({{m}_ \bullet })$| statistics for |${{m}_ \bullet } \ge 6$| from the global catalogue ANSS (1975–2015) and use the plane (X = log n, Y = log |$\bar{P}$|⁠). Such representation is useful to see the power type of the |${{V}_{{{m}_ \bullet } - 2}}({{m}_ \bullet })$| distribution. For control, we presented our |${{V}_{{{m}_ \bullet } - 2}}({{m}_ \bullet })$| distribution on the same (X,Y) plane (Fig. 3). The distribution looks like a parabola, similar to the data of Zaliapin & Ben-Zion (2016) for locations with low heat flux. Because these data dominate in the global catalogue, |${{V}_{{{m}_ \bullet } - 2}}({{m}_ \bullet })$| distributions in our and in Zaliapin & Ben-Zion (2016) cases are consistent. This should be kept in mind because different graphical representations emphasize different features of the distribution.

Survival function $P({{V}_{{{m}_ \bullet } - 2}}({{m}_ \bullet }) \ge n)$ for main shocks with ${{m}_ \bullet } \ge 6,{{m}_c} = 4$
Figure 3.

Survival function |$P({{V}_{{{m}_ \bullet } - 2}}({{m}_ \bullet }) \ge n)$| for main shocks with |${{m}_ \bullet } \ge 6,{{m}_c} = 4$|

|${{V}_{{{M}_a}}}({{m}_ \bullet })$|empirical distributions. We have taken all of the above into account to present the empirical |${{V}_{{{M}_a}}}({{m}_ \bullet })$| distributions for different |$\Delta \in (0.2,1.5)$| and |${{m}_a}$| models: (i) |${{m}_a} = {{m}_ \bullet } - 1.2$|⁠, (ii) |${{m}_a} = 1.03{{m}_ \bullet } - 1.51$|⁠. |${{V}_{{{M}_a}}}({{m}_ \bullet })$| accounts only for clusters for which it holds |${{m}_a}({{m}_ \bullet }) - \Delta \ge {{m}_c}$|⁠; the number of clusters identified for different models and different |$\Delta $| values are given in Table 1. The results for (i) model are shown in Fig. 4 and Table 2. Since Fig. 4 is similar for both models, model (i) is only presented in the final Fig. 4.

${{V}_{{{M}_a}}}({{m}_ \bullet })$ distributions for different $\Delta $ and the model ${{m}_a} = {{m}_ \bullet } - 1.2$. The Y-axis shows the survival function $\bar{P}(n) = \textit{prob}({{V}_{{{M}_a}}}({{m}_ \bullet }) \ge n)$. The dotted horizontal line corresponds to the 5 per cent level of $\bar{P}(n)$. ${{V}_{{{M}_a}}}({{m}_ \bullet })$ refers to main shocks with ${{m}_ \bullet } \ge 6$ for which it holds ${{m}_a}({{m}_ \bullet }) \ge {{m}_c} + \Delta $.
Figure 4.

|${{V}_{{{M}_a}}}({{m}_ \bullet })$| distributions for different |$\Delta $| and the model |${{m}_a} = {{m}_ \bullet } - 1.2$|⁠. The Y-axis shows the survival function |$\bar{P}(n) = \textit{prob}({{V}_{{{M}_a}}}({{m}_ \bullet }) \ge n)$|⁠. The dotted horizontal line corresponds to the 5 per cent level of |$\bar{P}(n)$|⁠. |${{V}_{{{M}_a}}}({{m}_ \bullet })$| refers to main shocks with |${{m}_ \bullet } \ge 6$| for which it holds |${{m}_a}({{m}_ \bullet }) \ge {{m}_c} + \Delta $|⁠.

Table 1.

Number of clusters, identified for the different models and values of |$\Delta $|⁠, used to obtain |${{V}_{{{M}_a}}}({{m}_ \bullet })$| distribution.

Model|$\Delta $|
0.20.50.81.01.5
|${{m}_a} = {{m}_ \bullet } - 1.2$|142714271427967367
|${{m}_a} = 1.03{{m}_ \bullet } - 1.51$|14271427967655243
Model|$\Delta $|
0.20.50.81.01.5
|${{m}_a} = {{m}_ \bullet } - 1.2$|142714271427967367
|${{m}_a} = 1.03{{m}_ \bullet } - 1.51$|14271427967655243
Table 1.

Number of clusters, identified for the different models and values of |$\Delta $|⁠, used to obtain |${{V}_{{{M}_a}}}({{m}_ \bullet })$| distribution.

Model|$\Delta $|
0.20.50.81.01.5
|${{m}_a} = {{m}_ \bullet } - 1.2$|142714271427967367
|${{m}_a} = 1.03{{m}_ \bullet } - 1.51$|14271427967655243
Model|$\Delta $|
0.20.50.81.01.5
|${{m}_a} = {{m}_ \bullet } - 1.2$|142714271427967367
|${{m}_a} = 1.03{{m}_ \bullet } - 1.51$|14271427967655243
Table 2.

Slopes |$\gamma (\Delta )$| of the |$- \log P({{V}_{{{M}_a}}}({{m}_ \bullet }) > n)$| curves in the 95 per cent intervals, eq. (35), of principal |${{V}_{{{M}_a}}}({{m}_ \bullet })$| values. Notation: (Th) theoretical slopes eq. (36) with b-value = 1; (i) |$\gamma (\Delta )$| for the model |${{m}_a} = {{m}_ \bullet } - 1.2$| (Fig. 4); (ii) |$\gamma (\Delta )$| for the model |${{m}_a} = 1.03{{m}_ \bullet } - 1.51$|⁠; (Sh) slopes for |${{m}_ \bullet } \ge 6.5$| based on Shebalin et al. (2018).

|$\gamma (\Delta )$||$\Delta $|
0.20.50.81.01.5
(Th)0.210.120.060.040.015
(i)0.220.100.050.030.015
(ii)0.220.100.040.030.013
(Sh, i)0.200.110.06
|$\gamma (\Delta )$||$\Delta $|
0.20.50.81.01.5
(Th)0.210.120.060.040.015
(i)0.220.100.050.030.015
(ii)0.220.100.040.030.013
(Sh, i)0.200.110.06
Table 2.

Slopes |$\gamma (\Delta )$| of the |$- \log P({{V}_{{{M}_a}}}({{m}_ \bullet }) > n)$| curves in the 95 per cent intervals, eq. (35), of principal |${{V}_{{{M}_a}}}({{m}_ \bullet })$| values. Notation: (Th) theoretical slopes eq. (36) with b-value = 1; (i) |$\gamma (\Delta )$| for the model |${{m}_a} = {{m}_ \bullet } - 1.2$| (Fig. 4); (ii) |$\gamma (\Delta )$| for the model |${{m}_a} = 1.03{{m}_ \bullet } - 1.51$|⁠; (Sh) slopes for |${{m}_ \bullet } \ge 6.5$| based on Shebalin et al. (2018).

|$\gamma (\Delta )$||$\Delta $|
0.20.50.81.01.5
(Th)0.210.120.060.040.015
(i)0.220.100.050.030.015
(ii)0.220.100.040.030.013
(Sh, i)0.200.110.06
|$\gamma (\Delta )$||$\Delta $|
0.20.50.81.01.5
(Th)0.210.120.060.040.015
(i)0.220.100.050.030.015
(ii)0.220.100.040.030.013
(Sh, i)0.200.110.06

Fig. 4 represents the |${{V}_{{{M}_a}}}({{m}_ \bullet })$| distribution by means of the survival function. Everything above the 5 per cent level of this function (i.e. above the black dashed line) is of interest for comparing with the theory. At all |$\Delta \in (0.2,1.5)$| the log-survival functions here are linear. This is consistent with the geometric distribution of |${{V}_{{{M}_a}}}({{m}_ \bullet })$| in the zone of principal values and is in agreement with the theoretical conclusion for the ETAS(F = G) model. This agreement automatically excludes the hypothesis F = P, for which the slope of the graph near 0 should be positive.

Table 2 allows us to compare the slopes of the curves with the theoretical ones in eq. (36). Here we have added the slopes from Shebalin et al. (2018) obtained for the main events |${{m}_ \bullet } \ge 6.5$|⁠. These data are related to model (i): |${{m}_a} = {{m}_ \bullet } - 1.2$| and are particularly important because the authors used the alternative method of Molchan & Dmitrieva (1992) to determine aftershocks.

Table 2 shows surprisingly good agreement of all empirical |$\gamma (\Delta )$| data with theory.

As a result, the |${{V}_{{{M}_a}}}({{m}_ \bullet })$| distributions for |${{m}_ \bullet } \ge 6$| agree with the theoretically predicted geometric distributions in zone of the principal |${{V}_{{{M}_a}}}({{m}_ \bullet })$| values and thus independently support the F = G hypothesis for direct aftershocks within the ETAS (F) model with |$\alpha /\beta $| close to 1.

In addition, Figs 3 and 4 show that the global data analyses from different authors are not inconsistent regarding the |${{V}_{{{M}_a}}}({{m}_ \bullet })$| distribution.

5 CONCLUSIONS

An important feature of a seismic process is the clustering of its events. In the ETAS model and its generalization ETAS(F), the dynamics of cluster formation is similar to a branching Galton–Watson type process. One way to test this mechanism can be based on the distributions of two integral characteristics: the number |$\nu ({{m}_ \bullet })$| of direct aftershocks (the input) and the total number |${{V}_M}({{m}_ \bullet })$| of cluster events with |$m \ge M({{m}_ \bullet })$| (the output).

The works of Shebalin et al. (2020, 2022) strongly suggests that the geometric F-distribution of |$\nu ({{m}_ \bullet })$| is preferable to the basic Poisson distribution in the ETAS model.

The empirical data on the |${{V}_M}({{m}_ \bullet })$| distribution are ambiguous, and hence the theoretical results are important. For the ETAS(F) model, they are obtained here for the first time even for the regular F = P case.

The key elements of the theoretical results are as follows:

  1. the choice of a special class of F distributions possessing the RT property. This choice simplifies the use of Utsu's law of productivity. Moreover, the ETAS(F) model with the RT property preserves the probabilistic structure when the lower magnitude threshold |${{m}_c}$| is changed. In this case only rescaling of the regular productivity is required (Proposition 2);

  2. the traditional selection of strong cluster events was changed: the lower |${{m}_ \bullet } - \Delta $| magnitude threshold was replaced by |$M({{m}_ \bullet }) = $||${{m}_a} - \Delta $|⁠, where |${{m}_a}$| is the distribution mode of the strongest aftershock.

  3. It is assumed that the initial event |${{m}_ \bullet }$| is sufficiently large.

In the subcritical regime, |${{V}_m}_{_a - \Delta }({{m}_ \bullet })$|-statistics has a quite simple limit distribution identical to the |$\nu ({{m}_\Delta })$| distribution with the mean value|$E\nu ({{m}_\Delta }) = {{e}^{\beta \Delta }}$|⁠.

Based on the hypothesis of geometric distribution of the number of direct aftershocks and using the global catalogue, we confirmed the theoretical result about the closeness of the |${{V}_m}_{_a - \Delta }({{m}_ \bullet })$| distribution to the geometric one in the region of the main values of this statistics.

This can be seen as a confirmation of the basic concept of the ETAS model about the epidemic structure of clusters and as an additional argument in favour of non standard ETAS(F) model with F = G.

AUTHOR CONTRIBUTIONS

George Molchan (Conceptualization, Formal analysis, Investigation, Supervision, Validation, Visualization, Writing – original draft, Writing – review & editing) and Antonella Peresan (Data curation, Investigation, Software, Visualization, Writing – original draft, Writing – review & editing).

Acknowledgement

Prof D. Sornette and Prof J. Zhuang stimulated the presented analysis of RT properties in seismicity models; the computational aspects of the work were fruitfully discussed with Dr E. Varini. We are indebted to Prof Y. Kovchegov for non-easy review and support of the work. The code for the nearest neighbor analysis belongs to Prof I. Zaliapin. This research received no external funding. The authors declare no conflict of interest.

DATA AVAILABILITY

The data analysed in this study, namely the global ANSS Comprehensive Earthquake Catalogue (ComCat) and related description, are available online via the USGS website (United States Geological Survey), at the following link: https://earthquake.usgs.gov/data/comcat (Last accessed on 1 april 2024). Additional information and the catalog with identified clusters are available upon request from the corresponding author (Antonella Peresan; [email protected]).

References

Denisov
D.
,
Foss
S.
,
Korshunov
D.
,
2010
.
Asymptotics of randomly stopped sums in the presence of heavy tails
,
Bernoulli
,
16
(
4
),
971
994
.

Duquesne
T.
,
Winkel
M.
,
2007
.
Growth of Levy trees
,
Prob. Theory Relat. Fields
,
139
(
3-4
),
313
371
.

Hawkes
A.G.
,
Oakes
D.
,
1974
,
A cluster process representation of a self-exciting process
,
J. Appl. Probab.
,
11
(
3
),
493
503
.

Kagan
Y.
,
2010
.
Statistical distributions of EQ numbers: consequence of brunching process
,
Geophys. J. Int.
,
180
(
3
),
1313
1328
.

Kagan
Y.
,
2014
.
Earthquakes: Models, Statistics, Testable Forecasts
,
283pp
,
AGU, Wiley
.

Kovchegov
Y.
,
Xu
G.
,
Zaliapin
I.
,
2023
.
Invariant Galton-Watson trees: metric properties and attraction with respect to generalized dynamical pruning
,
Adv. Appl. Probab.
,
55
(
2
),
643
671
.

Kovchegov
Y.
,
Zaliapin
I.
,
2021
.
Invariance and attraction properties of Galton-Watson trees
,
Bernoulli
,
27
(
3
),
1789
1823
.

Kovchegov
Y.
,
Zaliapin
I.
,
Ben-Zion
Y.
,
2022
.
Invariant Galton–Watson branching process for earthquake occurrence
,
Geophys. J. Int.
,
231
(
1
),
567
583
.

Li
J.
,
Sornette
D.
,
Wu
Z.
,
Zhuang
J.
,
Ch
J.
,
2024
.
Revisiting seismicity criticality: a new framework for bias correction of statistical seismology model calibrations
,
arXiv:2404.16374

Molchan
G.
,
Varini
E.
,
2024
.
The strongest aftershock in seismic models of epidemic type
,
Geophys. J. Int.
,
236
,
1440
1454
.

Molchan
G.
,
Varini
E.
,
Peresan
A.
,
2022
.
Productivity within the ETAS seismicity model
,
Geophys. J. Int.
,
231
,
1545
1557
.

Molchan
G.M.
,
Dmitrieva
O.E.
,
1992
.
Aftershock identification: methods and new approaches
,
Geophys. J. Int.
,
109
,
3
,
501
516
.

Ogata
Y.
,
1988
.
Statistical models for earthquake occurrence and residual analysis for point processes
,
J. Am. Statist. Assoc.
,
83
,
9
27
.

Saichev
A.
,
Helmstetter
A.
,
Sornette
D.
,
2005
.
Power law distributions of offspring and generation numbers in branching models of earthquake triggering
,
Pure. appl. Geophys.
,
162
,
1113
1134
.

Shebalin
P.
,
Baranov
S.
,
Dzeboev
B.
,
2018
.
The law of the repeatability of the number of aftershocks
,
Dokl. Akad. Nauk
,
481
(
3
),
963
966
. DOI:

Shebalin
P.
,
Baranov
S.
,
Vorobieva
I.
,
2022
.
Earthquake productivity law in a wide magnitude range
,
Front. Earth Sci.
,
10
,
881425
.

Shebalin
P.
,
Narteau
C.
,
Baranov
S.
,
2020
.
Earthquake productivity law
,
Geophys. J. Int.
,
222
,
1264
1269
.

Soloviev
S.
,
Solovieva
O.
1962
.
Exponential distribution of the total number of subsequent earthquakes and decreasing with the depth of its average value
,
Izv. AN USSR (Geophys.)
,
12
,
1685
1694
.

Stoilov
S.
,
1962
.
Theory of Functions of Complex Variable
.
Foreign Literature Publishing House
.

Zaliapin
I.
,
Ben-Zion
Y.
,
2013
.
Earthquake clusters in southern California I: identification and stability
,
J. geophys. Res..
,
118
(
6
),
2847
2864
.

Zaliapin
I.
,
Ben-Zion
Y.
,
2016
.
A global classification and characterization of earthquake clusters
,
Geophys. J. Int.
,
207
(
1
),
608
634
.

APPENDIX A: AUXILIARY STATEMENTS

Statement 5 (Molchan & Varini 2024).

Let |${{\nu }^{(\lambda )}},\lambda > 0$| be a family of random variables with the distribution: |$P({{\nu }^{(\lambda )}} = n) = {{p}_n}(\lambda ),n \ge 0$| such that |$E{{\nu }^{(\lambda )}} = \lambda $|⁠, |$E{{[{{\nu }^{(\lambda )}}]}^2} < \infty $| and |${{p}_0}(\lambda ) \to 0,\lambda \to \infty $|⁠.

Assume that for any |$\lambda > 0$|

Then

(a) |$\varphi (w),w < 0$| function is analytic, increasing from |$\varphi ( - \infty ) = 0$| to |$\varphi (0) = 1$|⁠;

(A1)
(A2)

Statement 6 (Molchan & Varini 2024).

Consider ETAS(F) model with magnitude and productivity characteristics |$({{f}_1}(m),\lambda (m))$|⁠. Let |${{\phi }_F}(z| m ) = E{{z}^{\nu (m)}}$| be the generating function for direct aftershocks of an m-event. Then |$DM({{m}_ \bullet })$| cluster in the ETAS(F) model can be viewed as |$AM({{m}_ \bullet })$| cluster in the ETAS(⁠|$\hat{F}$|⁠) model such that the magnitude law has the form

(A3)

and the off-spring law |$\hat{F}$| is contiguous with F:

(A4)

The |$\hat{F}$|generating function looks as follows

(A5)

APPENDIX B: PROOFS

Proof of Statement 1

Let |${{\phi }_n}(z)$| be generating functions оf integer variables |${{\nu }_n}$| with means |${{\lambda }_n} \uparrow \infty $| and |${{R}_p}$| is the random thinning operation with the parameter |$p \in (0,1)$|⁠. Then the generating function of |${{R}_{\lambda /{{\lambda }_n}}}{{\nu }_n}$| is |${{\psi }_n}(z| \lambda ) = {{\phi }_n}(1 + \lambda (z - 1)/{{\lambda }_n})$|⁠.

Since |${{\phi }_n}(w)$| is analytic and |$| {{{\phi }_n}(w)} | < 1$| in the disk |$| w | < 1$|⁠, the same properties hold for the function |${{\psi }_n}(z| \lambda )$| in the |${{S}_n}(\lambda )$| disk with the diameter |$(1 - 2{{\lambda }_n}/\lambda ,1)$|⁠.

Under these conditions, there exists a subsequence |${{\psi }_{\tilde{n}}}(z| \lambda )$| converging to an analytic function|${{\Psi }_N}(z| \lambda )$| in any disk with diameter (−N,1) (Stoilov 1962). These limits analytically continue each other and represent one analytic function|$\Psi (z| \lambda )$| in the half-plane |${\mathop{\rm{Re}}\nolimits} z < 1$|⁠.

Since |${{\psi }_n}(z| \lambda ) = {{\psi }_n}(1 + \lambda /\mu (z - 1)| {\mu )} $| on |${{S}_n}(\lambda \vee \mu )$|⁠, the same is true for the limit function|$\Psi (z| \lambda )$|. Hence |$\Psi (z| \lambda ) = \Psi (1 + \lambda (z - 1)| {1)} $|⁠.

The convergence of the generating functions |${{\psi }_{\tilde{n}}}(z| \lambda )$| in the disk |$| z | < 1$| entails convergence

|${{R}_{\lambda /{{\lambda }_n}}}{{\nu }_{\tilde{n}}}$|in distribution to an integer variable |$\nu (\lambda )$|⁠. By virtue of Fatou's theorem, |$E\nu (\lambda ) = \lambda $| since |$E{{R}_{\lambda /{{\lambda }_{\tilde{n}}}}}{{\nu }_{\tilde{n}}} = \lambda $|⁠. Hence, the |$\nu (\lambda )$|set has the RT property, since |$E{{z}^{\nu (\lambda )}} = \varphi (\lambda (z - 1))$| and |$E\nu (\lambda ) = \lambda $|⁠.

Proof of Statement 4

A.The basic equations forAM clusters

Let |$E{{z}^{\nu (m)}} = \varphi (\lambda (m)(z - 1))$|⁠, M is a fixed magnitude threshold,

(A6)

Using the main eq. (20), we have

(A7)

Now we fix the initial magnitude |${{m}_ \bullet }$|and |$M = M({{m}_ \bullet })$|⁠. Consider

(A8)

Note that |${{a}_u}(z) \ge {{a}_0}(z)$| since |${{U}_M}(z| {m)} \le 1$|⁠. We can rewrite eq. (A7) as follows:

(A9)

where

(A10)

Let us subtract 1 from both parts of equality eq. (A8) and integrate them over interval |$(u,\infty )$| with weight |$\lambda ({{m}_0}){{f}_1}(m)$|⁠. We obtain a closed system of two equations with respect to |${{a}_u}(z),u = 0,M:$|

(A11)

By eq. (A9), the solution of our problem is the limit of

(A12)

If a non-trivial limit of |${{U}_M}(z| {{{m}_ \bullet })} $| exists, the limit|${{A}_M}(z)$| function is separated from 0 and |$\{ - \infty \} $|⁠, where |$\varphi ( \cdot )$|is 1 and 0, respectively. So we will search for a solution under the assumption that at fixed |$z \in (0,1)$|

(A13)

Let's substitute

into eq. (A11) and change the variables under the integrals: |$m \to x = \exp (\beta ({{m}_ \bullet } - m))$|⁠. Omitting the |$(1 - {{e}^{ - \beta {{M}_1}}}) = 1 + o(1)$| term we have

(A14)

where

We have |$y(x) = {{A}_M}(z){{x}^{ - \alpha /\beta }} \to 0$|as |$x \to \infty $| by virtue of eq. (A13), and |$\varphi (y) - 1 = y(1 + o(1)),y \to 0$| by virtue of Statement 5. Applying the L'Hôpital's rule, we have for large R

(A15)

In addition, since|$0 \le {{\varphi }_F} \le 1$|⁠, we have

(A16)

Relations (A15A16) are key to simplify eq. (A14) when |${{m}_ \bullet }$| and |$M = {{m}_a} - \Delta $| are large.

A1. The case |$(\alpha < \beta ,n < 1)$|

In this case

(A17)

and |$R = {{e}^{\beta ({{m}_ \bullet } - u)}} > > 1$| both for u = 0 and u = M. Therefore eq. (A15) give

(A18)
(A19)

Substituting eqs (A18A19) in eq. (A10), we have

That is |${{A}_M}(z) \approx (z - 1){{e}^{\beta \Delta }}$|⁠. As a result we get

(A20)

A2. The case|$(\alpha = \beta ,n < 1)$|

In this case |$n = {{\lambda }_0}\beta {{M}_1} = {{\lambda }_0}\beta {{m}_ \bullet }{{K}_1}$| and

(A21)

Therefore |$R = {{e}^{\beta ({{m}_ \bullet } - u)}} > > 1$| both for u = 0 and u = M since |${{\lambda }_0}{{m}_ \bullet } = \textit{const}$|⁠. Using eq. (A10), we have

(A22)
(A23)

and by eq. (A10),

Hence, |${{U}_M}(z| {{{m}_ \bullet })} $| limit is again eq. (A20).

A3.The case |$(2\alpha \le \beta ,n = 1)$|

In this case

(A24)

where |$B = 2{{(\beta {{\varphi ^{\prime\prime}}_F}(0))}^{ - 1}}(\beta - 2\alpha + {{1}_{\beta = 2\alpha }})$|⁠. Since |$R = \beta ({{m}_ \bullet } - M) > > 1$|⁠, eq. (A14) for |${{a}_M}(z)$| is analysed as above using eq. (A15) asymptotics. We have

(A25)

Now we consider the equation for|${{a}_0}(z)$|⁠. Since

we have

(A26)

To simplify eq. (A26), we will need the properties of |${{\varphi }_F}(w)$| from Statement 5:

(A27)
(A28)

Acting in the same way as in (A16A17) we have

(A29)
(A30)

Substituting |$R = {{e}^{\beta {{m}_ \bullet }}}$|⁠, eq. (A26) is transformed by eq. (A29) as follows:

(A31)

Substituting here M from eq. (A24), we obtain |$A_M^2(z) = (1 - z){{e}^{\beta \Delta }}$|⁠. Therefore,

(A32)

A4.The case |$(\alpha < \beta < 2\alpha ,n = 1)$|

We have

(A33)

where |$\psi ({{\lambda }_0})$| is given in eq. (17). Hence, by eqs (A14) and (A27A28)

(A34)

For |${{a}_0}(z)$|we use eq. (A26) to have

(A35)

Due to eqs (A27A28) the following integral

(A36)

is finite if|$(\alpha < \beta < 2\alpha $|⁠. Therefore, in the limit the integration interval in eq. (A28) can be replaced by a semi-straight line. Then, eqs (A34) and (A35) give

Finally,

(A37)

Here |$\exp ( - \alpha \Lambda ) = \psi ({{\lambda }_0}){{e}^{ - \alpha \Delta }}$|⁠. The obtained equation can be simplified at small |${{\lambda }_0} = (1 - \alpha /\beta )n$|⁠. In this case |$\psi ({{\lambda }_0}) = {{\lambda }_0}(1 + o(1))$| and |$\int_{{_\varepsilon }}^{\infty }{{{{\varphi }_F}}}( - u){{u}^{ - \beta /\alpha - 1}}du = \alpha {{\beta }^{ - 1}}{{\varepsilon }^{ - \beta /\alpha }}(1 + o(1)$|as |$\varepsilon \to 0$|⁠.

Hence at small |${{\lambda }_0}$|⁠,

(A38)

that is

The relation eq. (A37) at z = 0. This case is equivalent to (17, 18), that is

(A39)

where |$\psi = | {{{A}_M}} (z = 1 ) |{{e}^{ - \alpha \Delta }}$|⁠. To verify this, it is enough to substitute |${{I}_F}$| from eq. (A36) into eq. (A37). After simple algebra we obtain at z = 0

To get eq. (A39) we should multiply the obtained ratio by |${{\psi }^{\beta /\alpha }}(\beta /\alpha - 1)$| and take into account that |${{\lambda }_0} = 1 - \alpha /\beta $|⁠.

B. DM clusters

Now we consider |$DM({{m}_ \bullet })$| cluster in ETAS(F) model with RT property of F. According to Statement 6, statistics |$\nu (m)$| in such cluster has the following generating function

(A40)

Following eqs (A7)–(A11) we can obtain the equation for

where |${{\hat{f}}_1}(m) = {{f}_1}(m)/{{F}_1}({{m}_ \bullet }),0 \le m \le {{m}_ \bullet }$|⁠.

Using notation

one has

(A41)

where

(A42)

Let's subtract 1 from both parts of equality eq. (A41) and integrate over interval |$(u,{{m}_ \bullet })$| with weight |$\lambda ({{m}_ \bullet }){{f}_1}(m)$|⁠. Then we obtain a closed system of two equations with respect to |${{\hat{a}}_u}(z),u = 0,M$|⁠:

(A43)

After replacement of variables |$m \to x = \exp (\beta ({{m}_ \bullet } - m))$| we obtain the analogue of relation eq. (A14)

(A44)

where |$\varepsilon = {{\lambda }_0}{{e}^{ - (\beta - \alpha ){{m}_ \bullet }}}$|⁠.

In the case |$\alpha = \beta $|⁠,|${{\lambda }_0}{{m}_ \bullet } = \textit{const}$|⁠. Therefore|$\varepsilon = o(1)$|⁠, when |$\alpha \le \beta $| and |${{m}_ \bullet } > > 1$|⁠.

Our goal is the limit function for |${{U}_M}(z| {{{m}_ \bullet })} = $||${{\varphi }_F}({{\hat{A}}_M}(z))$||$/{{\varphi }_F}( - \varepsilon )$|⁠.

Since |$1 - x < {{\varphi }_F}( - x) \le 1$|⁠,

(A45)

Hence, as above, we can assume that at fixed |$z \in (0,1)$|

To simplify eq. (A43), note that.

From here

(A46)

where |$u \le {{m}_ \bullet }$|and

Therefore

(A47)

Relations (eqs A46A47) allow us to transfer the proofs of the results for AM clusters to DM clusters in the following regimes: |$(\alpha \le \beta ,n \le 1) \cup (2\alpha \le \beta ,n = 1)$|⁠, where M is the same for these clusters.

B2. The case |$(\alpha < \beta < 2\alpha ,n = 1)$|⁠, DM cluster

Let's remember that|$M = {{m}_ \bullet } - \Delta $|⁠, |${{\lambda }_0} = 1 - \alpha /\beta $|⁠, |$\varepsilon = {{\lambda }_0}{{e}^{ - (\beta - \alpha ){{m}_ \bullet }}}$|⁠.

Replacing in eq. (A44) the variable x by |$y = {{x}^{ - \alpha /\beta }}$| we have

(A48)

where |$\mu (dy) = \beta /\alpha {{y}^{ - \beta /\alpha - 1}}dy$|⁠.

By eq. (A42),

(A49)

By eq. (A48),

(A50)

where

(A51)
(A52)

The estimations in eqs (A51A52) follow from eq. (A27). Each of the relations eqs (A49), (A50) allows us to find an expression for |$Q ({{\hat{A}}_M}(z) - {{\hat{a}}_0}(z))/\varepsilon + 1$|⁠.

According to eq. (A48),

Given eq. (A50), we have

Since |${{I}_F} = \int_{0}^{\infty }{{[{{\varphi }_F}}}( - u) - 1 + u]{{u}^{ - \beta /\alpha - 1}}dx$| is bounded and |${{\varphi }_F}( - y\varepsilon ) - 1 = O(\varepsilon )$| we can go to the limit to get the equation for the limit |${{\hat{A}}_M}(z)$| function:

(A53)

Finally, the desired generating function of |${{V}_M}({{m}_ \bullet })$| is |${{\varphi }_F}({{\hat{A}}_M}(z))$|⁠.

The case of small|${{\lambda }_0}$|⁠. When |${{\lambda }_0}$| is small, eq. (A53) can be simplified. Let's substitute |${{\hat{A}}_M}(z)/{{\lambda }_0} = B(z)$| into eq. (A53), then we obtain

For |${{\lambda }_0} < < 1$|

that is

Hence

(A54)

The relation (eq. A53) at z = 0. This case is equivalent to eq. (17), that is

(A55)

where |$\psi = | {{{{\hat{A}}}_M}} (z = 0 ) |{{e}^{ - \alpha \Delta }}$|⁠. To see this, let us split the integration interval (0,1) in eq. (A53) into two |$(0,\rho ) \cup (\rho ,1)$|⁠, where |$\rho = {{e}^{ - \alpha \Delta }}$| Then eq. (A53) at z = 0 is as follows:

(A56)

where |$A = {{\hat{A}}_M}(z = 1)$| and

(A57)

Let us substitute eq. (A57) into eq. (A56) and multiply both parts of the equality (eq. A56) by |${{\lambda }_0}{{\rho }^{\beta /\alpha }}$|⁠. As a result, we obtain eq. (A55). From here, for small |${{\lambda }_0} = 1 - \alpha /\beta $|

which is consistent with eq. (A54) at z = 0.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.