Productivity within the epidemic-type seismicity model

SUMMARY

The productivity of a magnitude m event can be characterized in term of triggered events of magnitude above m − Δ: it is the number of direct ‘descendants’ |$\nu _\Delta$| and the number of all ‘descendants’ |$V_\Delta$|⁠. There is evidence in favour of the discrete exponential distribution for both |$\nu _\Delta$| and |$V_\Delta$| with a dominant initial magnitude m (the case of aftershock cluster). We consider the general Epidemic Type Aftershock Sequence model adapted to any distribution of |$\nu _\Delta$|⁠. Our first result shows that models with branching aftershock structure do not allow for the coincidence of distribution types of |$\nu _\Delta$| and |$V_\Delta$| (say, the discrete exponential, as in the scientific literature). The second problem is related to the tail behaviour of the |$V_\Delta$| distribution. We show the fundamental difference in tail behaviour of the |$V_\Delta$|-distributions for general-type clusters and clusters with a dominant initial magnitude: the tail is heavy in the former case and light in the latter. The real data demonstrate the possibilities of this kind. This result provides theoretical and practical constraints for distributional analysis of |$V_\Delta$|⁠.

Probability distributions, Earthquake interaction, forecasting, and prediction, Statistical seismology

1 INTRODUCTION

Earthquakes clustering, aftershocks occurrence in particular, is a prominent feature in observed seismicity. Different models have been proposed to describe the features of triggered events, especially in the case of aftershocks clusters with an initial dominant magnitude earthquake; yet, several theoretical aspects still need to be explored. In this study, we investigate the distributions of the number of triggered events as described by a branching process, which might provide useful practical constraints for empirical analysis of earthquake data.

Let |$V_\Delta$| denote the number of seismic events of magnitude m ≥ m₀ − Δ caused by an earthquake of magnitude m₀ ≥ m_c + Δ, where m_c is the lower magnitude threshold and Δ a positive default value. |$V_\Delta$| is further referred to as the total Δ-productivity.

Considering an earthquake cluster as the result of a branching process (Zaliapin et al. 2008), we can also deal with the directly triggered events of m₀, that is the first generation events of an m₀-event. The number |$\nu _\Delta$| of such events is in the same magnitude range m ≥ m₀ − Δ and is later named Δ-productivity of m₀.

The distributions of quantities |$V_\Delta$| and |$\nu _\Delta$| have recently become the object of active analysis in the literature (Baranov & Shebalin 2019; Shebalin et al. 2020; Shebalin & Baranov 2021; Baranov et al. 2022). According to these publications, the |$\nu _\Delta$|-distribution |$f(n)=P(\nu _\Delta =n)$| is universal in shape, namely, it is geometric (or discrete exponential):

$$\begin{eqnarray*} f_G(n)=p^n\left(1-p\right) n=0,1,\dots , \end{eqnarray*}$$

(1)

where p does not depend on the magnitude of triggering event. This result is relevant in statistical seismology because the popular Epidemic Type Aftershock Sequence (ETAS) model (Ogata 1998) is based on the Poisson distribution

$$\begin{eqnarray*} f_P(n)=e^{-\lambda }\lambda ^n / n! n=0,1,\dots \ . \end{eqnarray*}$$

(2)

Similarly, the geometric distribution (eq. 1) was adopted for statistics |$V_\Delta$|⁠, provided that it refers to clusters with a dominant initial magnitude m₀, that is, to aftershocks. For the first time this hypothesis appeared in the work by Soloviev & Solovieva (1962); it was based on relatively scarce data (1954–1961) for the Pacific Belt (m₀ ≈ 7; Δ = 2) and for the Kamchatka and Kuril Islands region (m₀ ≈ 6; Δ = 2). Shebalin et al. (2018) expanded the analysis and confirmed the hypothesis using 850 aftershock sequences with main shocks m₀ ≥ 6.5, Δ = 2 from the ANSS (1975–2018) catalogue. In the same time Kagan (2010), using aftershocks m ≥ 4.7 from PDE (1977–2007) catalogue for m₀ = 7.1−7.2, found that the empirical distribution of |$V_\Delta$| is bimodal with the dominant peak at n = 0.

Zaliapin & Ben-Zion (2016) provided a thorough analysis of earthquake clustering and showed its spatial dependence and connection with heat flow production. They showed that, on the global scale, the total Δ-productivity distribution (for clusters with max m ≥ 6 and Δ = 2) looks like a power-law for locations with high heat flow (H > 0.2 W m^–2), and as a log-normal distribution otherwise. These conclusions are qualitative because they are based on log–log plots. Nevertheless, the data for the regions of high heat flow differ significantly from the geometric hypothesis (eq. 1) for aftershocks. The same should be expected for any cluster with max m ≤ 5. In this case, according to Zaliapin & Ben-Zion (2016), the cluster size distribution in zones with both high and low heat flow is a power type distribution, f(n)∝n^{−κ − 1}, κ ≈ 2.5.

These observations rise the following questions: (i) is it possible to consider the geometric distribution of |$\nu _\Delta$| and |$V_\Delta$| (for aftershocks) as laws of seismicity and (ii) how to explain the observed difference in the tail behaviour of the |$V_\Delta$| distribution for arbitrary clusters and aftershocks. A rigorous analysis of these issues is possible within the framework of a model adequately related to the concepts of |$\nu _\Delta$| and |$V_\Delta$|⁠. A model suitable for this purpose may be the general ETAS model, namely the ETAS model with an arbitrary (not only Poisson) distribution of direct aftershocks associated to any event (Saichev et al. 2005). The purpose of this note is twofold: first, to show the incompatibility of the hypotheses about the geometric distribution of |$\nu _\Delta$| and |$V_\Delta$| (for aftershocks) within the general ETAS model; and second, to prove the fundamental difference between size distributions for general type clusters and aftershock clusters.

2 ETAS MODELS

Let us consider the seismic events x = (t, g, m) occurred in the time–location–magnitude space |$\mathcal {S}=(t\ge 0, g\in \boldsymbol {R}^2, m\ge 0)$|⁠. The lower magnitude threshold m_c is taken here as the zero reference point without loss of generality. According to the branching structure of the ETAS model, an earthquake cluster can be represented as follows (Saichev et al. 2005). The initial event x₀ = (t₀, g₀, m₀) generates a random number ν(m₀) of events {x_i: i = 1, …, ν(m₀)}, where ν(m₀) follows a special distribution F with an average λ(m₀) = Eν(m₀). The offspring events are distributed in |$\mathcal {S}$| as independent variables with the probability density

$$\begin{eqnarray*} f(x\mid x_0) = f_1(m) f_2(t-t_0) f_3(g-g_0\mid m_0) \quad . \end{eqnarray*}$$

(3)

Each new event independently generates its own family of events following the generation law described above, and so on. If the mean of the direct triggered events is Λ₁ = Eλ(m) < 1, then the mathematical expectation of the cluster size will be finite.

The independence of the magnitudes of direct ‘descendants’ between themselves and from the parent event x₀ is an important property of the model. This allows us to consider the distribution f₁(m) as a frequency–magnitude law for the cluster as a whole. This also is important for understanding the choice of the distribution models for |$\nu _\Delta$| and |$V_\Delta$|⁠.

The components of the regular ETAS model are specified as follows:

$$\begin{eqnarray*} \begin{array}{ll}\nu (m) {-distribution} & \mbox {- Poisson;} \\ \lambda (m) = \lambda e^{\alpha m}, \ \ m \ge 0 & \mbox {- a copy of the Utsu law for aftershocks;} \\ f_1(m)=\beta e^{-\beta m}, \ \ m\ge 0, \beta \gt \alpha & \mbox {- the Gutenberg-Richter law,} \\ f_2(t)=(p-1)(t/c+1)^{-p}/c & \mbox{- the time scattering (a copy of the Omori} \\ & \mbox {law for aftershocks);} \\ f_3(g\mid m)=(q-1)\left[|g/\psi _m|^2+1\right]^{-q}/(\pi \psi _m^2) & \mbox {- one of the laws of scattering in space;} \\ \psi _m=de^{\gamma m/2} & \mbox{- the scale parameter}. \end{array} \end{eqnarray*}$$

Note: The exponents (α, β) in the laws of Utsu and Gutenberg–Richter are more convenient (Saichev et al. 2005) to use with the base ‘e’ instead of the usual 10.

A theoretical analysis of various aspects of the general ETAS model with the regular components λ(m) and f₁(m) is provided in the works of Saichev and co-authors (see Saichev & Sornette 2017, for further details). In particular, they proved the power law behaviour of the cluster size distribution f_cl(n)∝n^{−1 − κ} with κ ∈ [0.5, 1] (Saichev et al. 2005; Baró 2020). Note that, according to Zaliapin & Ben-Zion (2016), κ ∈ [1, 2.5] for real data.

In the general case we need to specify the family of ν(m) distributions: F = {f(n∣m)}. This choice will be reflected in the designation of the model as ETAS(F). In particular, we will set F = P for Poisson, and F = G for Geometric distribution. The component λ(m) = Eν(m) will be common to all models. The relation of λ(m) with the parameter p(m) of the geometric distribution is determined by the equality λ = p/(1 − p), that is

$$\begin{eqnarray*} p(m)=\lambda (m)/(1+\lambda (m))\ . \end{eqnarray*}$$

(4)

Recall that the regular ETAS model can be defined in terms of the conditional intensity of events (Hawkes & Oakes 1974):

$$\begin{eqnarray*} \lambda ^{*}\left( x \mid \lbrace x_i : t_i\lt t\rbrace \right) = \mu (g) f_1(m) + \sum _i \lambda (m_i) f(x\mid x_i) \end{eqnarray*}$$

(5)

where x = (t, g, m) and μ(g) is the rate of background seismicity. The general definition is broader because it allows for any ν(m) distribution.

To provide the G-distribution of Δ-productivity for any Δ, it is natural to consider the G-distribution of ν(m). It follows from the following seemingly obvious but useful statement.

Statement. Let the distribution type F of ν(m₀) be F = G or P. Then |$\nu _\Delta (m)$|-distribution will have the same type for m₀ ≥ Δ. In both cases the statistics |$\nu _\Delta (m_0)$| and ν(m₀) differ only in average values:

$$\begin{eqnarray*} \lambda _\Delta (m_0) = \lambda (m_0) \bar{F}_1(m_0-\Delta ) / \bar{F}_1(\Delta ) \end{eqnarray*}$$

(6)

where |$\bar{F}_1(m) = \int _m^\infty f_1(u){\rm d}u$|⁠, |$\bar{F}_1(0) = 1$|⁠. Let a = f₁(0) > 0 and |$\lambda _\Delta (m_0)$| be independent of m₀ in the area 0 < Δ ≤ m₀ < M. Then λ(m) = ce^am and |$\bar{F}_1(m)= e^{-am}$|⁠, m ∈ (0, M).

The proof of the statement and all subsequent ones are contained in the  Appendix.

 

Remark 1.

The conditions for the independence of the |$\nu _\Delta (m_0)$|-distribution from m₀ are very strict, and therefore, in contrast with Shebalin et al. (2020), this property should be unstable.

 

Remark 2.

The Random Thinning (RT) property. Statement 1 allows the following generalization. Consider the RT(π) operation of a random cluster thinning: RT(π) saves each cluster event independently with probability π, and deletes it with probability 1 − π; the transformed cluster size V is RT(π)V. Given the magnitude independence in the ETAS(F) model, the productivity |$\nu _\Delta$| has the same distribution as RT(π)ν provided that π = P(m ≥ m₀ − Δ).

In the case F = (P or G), E[RT(π)ν] = πEν and the distribution type of ν and RT(π)ν is identical. We will call these two properties of the discrete distribution the RT property. We can formalize the RT property as follows.

Let V_λ and V_λ(π) = RT(π)V_λ be the cluster sizes, where λ = EV_λ. If |$Ez^{V_\lambda }=\varphi (z,\lambda )$|⁠, then by eq. (A1) from Appendix, |$Ez^{V_\lambda (\pi )}=\varphi (1-\pi +\pi z, \lambda )$|⁠, EV_λ(π) = λπ. The RT property of V_λ implies

$$\begin{eqnarray*} \varphi (1-\pi +\pi z, \lambda )=\varphi (z,\lambda \pi )\ , \ \ \ \lambda \lt \lambda _0\le \infty \end{eqnarray*}$$

Let for simplicity λ₀ ≥ 1. Setting π = Λ/n < 1, λ = n and ψ⁽ⁿ⁾(z) = φ(1 + z/n, n), we get:

$$\begin{eqnarray*} \psi ^{(n)}(\Lambda (z-1)) = \varphi (z,\Lambda )\ ,\ \ \ \Lambda \le n\lt \lambda _0\ . \end{eqnarray*}$$

(7)

Due to the analyticity of φ(z, λ) in the disk ∣z∣ < 1, ψ⁽ⁿ⁾(z) is analytical in K_n = {z: ∣1 + z/n∣ < 1}. Moreover, (7) implies that ψ⁽ⁿ⁾(z) and ψ^{(n + 1)}(z) are equal in K_n∩K_{n + 1} = K_n. Hence ψ⁽ⁿ⁾(z) is an analytical continuation of ψ⁽¹⁾(z) = φ(1 + z, 1). The function ψ(z) = φ(1 + z, 1) specifies the type of F in the ETAS(F) model for which Statement 1 remains true. For example, ψ(z) = e^z and ψ(z) = (1 − z)⁻¹ correspond to the types of Poisson and Geometric distributions, respectively.

Examples of the RT property.

• The Negative Binomial Distribution (NBD) for whichφ(z, λ) = (1 − λ(z − 1))^−τ, τ > 0, holds the RT property. In the case τ = 1 it corresponds to the Geometric distribution.

• Combining Nindependent clusters withRT property also preserves this property. Formally this follows from eq. (7) and the relation for the cluster characteristics φ_i(z, λ_i):

  $$\begin{eqnarray*} \varphi (z,\lambda ) = \prod _i \varphi _i(z,\lambda _i) = \prod _i \psi _i(p_i\lambda (z-1)), \quad \quad p_i\lambda =\lambda _i \quad {\rm and} \quad \sum p_i=1 \end{eqnarray*}$$

  (8)

  The NBD model was used by Kagan (2014) to describe the distribution of annual earthquake numbers for the CIT catalogue 1932–2001, m ≥ 5.0. Below we will show that the RT property, which NBD possesses, is poorly compatible with the branching cluster structure. However, in the Kagan’s case, the statistics of events occurring worldwide in one year consists mainly of independent weakly grouped events. Such case is closer to the model (eq. 8) with components: φ_i(z, λ_i) = 1 + λ_i(z − 1), λ_i < 1.
• A cluster selected randomly with probabilities {ρ_i > 0} from N clusters withRT property retains RT property. In this case

  $$\begin{eqnarray*} \varphi (z,\lambda )=\sum _i \varphi _i(z,\lambda _i)\rho _i =\sum _i \psi _i(p_i\lambda (z-1))\rho _i , \quad \quad \lambda p_i=\lambda _i \quad {\rm and} \quad \sum _i \lambda _i \rho _i = \lambda . \quad \end{eqnarray*}$$

  (9)

  However, if the cluster sizes distributions are of the same type, that isφ_i(z, λ_i) = ψ(λ_i(z − 1)), then the random cluster usually loses this type provided that λ_i ≠ const.

  For example, assuming a geometric distribution for both the size of a random cluster and the sizes of its components, (9) gives (1 − w)⁻¹ = ∑_i(1 − wp_i)⁻¹ρ_i, w = λ(z − 1), which is possible only when p_i = 1. This is also true in general case if φ(z, λ_i) is N-smooth at z = 1. The considered case is typical for the observed clusters taken as one statistical population. Consequently, empirical distributions of |$\nu _\Delta$| or |$V_\Delta$| obtained over large territories generally do not represent the true types of their distributions. This circumstance makes it difficult to prove the universality of the type of |$\nu _\Delta (m_0)$| distribution.

 

Remark 3.

ETAS* model. To get the geometric distribution of |$\nu _\Delta$| in framework of the ETAS model, Shebalin et al. (2020) used randomization of the parameter λ(m₀). They used the fact that a Poisson distribution with an exponentially distributed random parameter is equivalent to a geometric distribution. The resulting model was named ETAS*. Its analogue is the ETAS(G) model in which the average productivity does not depend on the initial magnitude, namely λ(m₀) = λ. Because of (6), in this case |$\lambda _\Delta (m_0)$| becomes dependent on m₀ for any f₁(m) > 0.

The analogue of the ETAS* model may be interesting due to its connection with the classical branching Galton–Watson model (Harris 1963). Both models are identical if we omit the space–time component in the model. This gives us an exact formula for the distribution of the total productivity of any initial event (cluster size not including the initial event):

$$\begin{eqnarray*} P(V=n) &=& \Gamma (1/2+n)/\Gamma (2+n) \times (4pq)^n (q/\sqrt{\pi })\ , \quad q=1-p , \end{eqnarray*}$$

(10)

where Γ(x) is the Gamma function, and p < 1/2 is the parameter of the geometric distribution (Harris 1963). Formula (10) remains valid in the critical regime when λ = p/(1 − p) = 1 or p = 1/2. The distribution (eq. 10) has the following asymptotic behaviour:

$$\begin{eqnarray*} P(V=n) &\propto & n^{-3/2} \exp(-C_1 n)\ , \quad n \gg 1 \end{eqnarray*}$$

(11)

where C₁ = −ln (4pq)  ≥  0. Such asymptotic with the factor n^−3/2 is typical for the Galton–Watson models (Harris 1963) and, consequently, for any ETAS(F) model with λ(m) = const.

The generating function of the distribution (eq. 10) is

$$\begin{eqnarray*} Ez^{V} = (1+2\Lambda -\sqrt{1-4(1+\Lambda )\Lambda (z-1)}) / (2\Lambda z), \end{eqnarray*}$$

(12)

where EV = Λ = λ/(1 − λ) = p/(1 − 2p) (see eq. 4). It is obvious that eq. (12) cannot be represented as eq. (7). Therefore, the analogue of ETAS* model has the geometric distribution for the Δ-productivity |$\nu _\Delta$|⁠, but does not have any RT property for the total Δ-productivity |$V_\Delta$|⁠.

3 THE TOTAL Δ -PRODUCTIVITY

We will consider the ETAS(F) model in a subcritical regime: Λ₁ = Eλ(m) < 1. Eν(m) = λ(m) is a non-decreasing, finite, positive function; the other components, namely f₁(m), f₂(t) and f₃(g∣m), are free from additional restrictions.

As above m₀ > 0 is the magnitude of an initial cluster event and |$V_\Delta (m_0)$| is the total Δ-productivity of m₀. Any initial event can be fixed or random, namely it has the f₁(m) distribution. A further distinction concerns the dominant initial events, that is those initial events having maximum magnitude within their cluster. The case of a random m₀ is important because identifying the initial cluster event is an unstable practical procedure.

Statement 2. Consider the specified ETAS(F) model. Assume that all moments of F-distribution are finite and F has the RT property.

• Let |$V_\Delta (m_0)$|⁠, m₀ > Δ > 0, be the size of the cluster with any condition C: the initial magnitude m₀ is fixed or random. Then the number |$N_{V_\Delta }$| of finite moments of |$V_\Delta (m_0)$| is independent of C. Moreover,

  $$\begin{eqnarray*} N_{V_\Delta } = N_\lambda := max\left\lbrace n: I_n = \int _0^\infty \lambda ^n(m) f_1(m) dm \lt \infty \right\rbrace \ . \end{eqnarray*}$$
• Let F be Poisson or Geometric distribution and m₀ is fixed and dominant. Then |$N_{V_\Delta }=\infty$| regardless of the N_λ-value. This is also true for the random dominant m₀-event if additionally

  $$\begin{eqnarray*} \lambda (m) \int _0^\Delta f_1(m-u) {\rm d}u \le {\rm const}\lt \infty \ . \end{eqnarray*}$$

  (13)

 

Remark 4.

Statement 2 allows us to judge on the tail attenuation of the |$V_\Delta$| distribution because, by the Chebyshev–Markov inequality, P(ξ ≥ u) ≤ Eξⁿ/uⁿ for any random variable ξ ≥ 0 and u > 0. Distributions having all moments will be called distributions with light tails. Otherwise, we will talk about heavy tails.

Let’s consider clusters with fixed or random initial magnitude. The conditions of the Statement 2 are easily verified for the regular ETAS components: λ(m) = λe^αm, α > 0 and f₁(m) = βe^−βm, β > α. In this case we have N_λ < β/α. This means that the total Δ-productivity distribution has only a heavy tail for small β/α. The data of Zaliapin & Ben-Zion (2016) demonstrate the possibilities of this kind. In the considered case Saichev et al. (2005) give more information on the tail behaviour of |$V_\Delta$| distribution: if β/α is large, then the tail is subexponential and looks like eq. (11).

The case α = 0 corresponds to the analogue of the ETAS* model; here N_λ = ∞ in accordance with eq. (10).

Given α ≤ β, condition (13) is satisfied. This guarantees the light tails in the distribution of the total Δ-productivity of main shocks. All the observations mentioned above support this fact.

 

Remark 5.

The difference in tail behaviour of |$V_\Delta$| statistics can be partially explained as follows. The total Δ-productivity of the main shock with a fixed magnitude m₀ is based on events in the finite range of magnitudes (0, m₀). In this case, the frequency–magnitude law f₁(m) is replaced by its truncated counterpart in the interval (0, m₀), where all moments of λ(m) are finite. Since |$N_{V_\Delta }=N_\lambda$|⁠, the tail of the |$V_\Delta$| distribution becomes typical for the case N_λ = ∞ that is quasi exponential. The similar effect for a random dominant magnitude is not so obvious and is a non-trivial analytical fact.

Statement 3. (1) Let m₀ be a dominant initial magnitude in the ETAS(G) cluster and let f₁(m) be strictly positive. Then the statistics |$V_\Delta (m_0)$| can have a geometric distribution for no more than one value Δ ∈ (0, m₀).

(2) Let ν and V be the sizes of clusters of direct descendants and all descendants, respectively, in the ETAS(F) cluster with the random initial event. Assume that both distributions of ν and V have the RT property, as well as Eν² < ∞ and ∫λ²(m)f₁(m)dm < ∞. Then the distribution types of ν and V cannot be the same.

 

Remark 6.

According to Statement 3 we can say that, in the framework of generalized ETAS model, the geometric laws for |$\nu _\Delta$| and |$V_\Delta$|⁠, including the case of aftershocks, are incompatible. In the ETAS(G) model case with λ(m) = const., the geometric law of |$\nu _\Delta$| completely excludes the RT property for |$V_\Delta$| (see eq. 12). The reason for this lies in the branching structure of the cluster.

Numerical examples

To support the theoretical conclusions, we simulated a synthetic earthquake catalogue from the regular ETAS(P) model. To mimic the features of real seismicity, the set of ETAS parameters was mainly derived from the analysis of the earthquake catalogue of Northeastern Italy (Benali et al. 2020). The simulation region also corresponds to that area, namely: Lon 11.5–14.0 and Lat 45.5–47.0.

The model was based on the following parameters:

1. The uniform background seismicity with the rate μ = 0.04 events (degree²·day)^–1, which roughly corresponds to the rate of background events identified in the earthquake catalogue of Northeastern Italy. The uniformity makes it easier the declustering process, as it allows to better identify the clustered component.

2. Magnitude range m ≥ 2; the Gutenberg–Richter parameter β = 2.07, which corresponds to the b-value 0.9.

3. The Utsu Law parameters: λ = 0.22 (event/day), α = 1.54. The estimated λ = 0.67 for Northeastern Italy was reduced to 0.22 to satisfy the stability conditions Λ₁ = Eλ(m) < 1.

4. The Omori law parameters: c = 0.015 (day), p = 1.037. The estimated p parameter very close to 1 implies an extremely long tail of f₂(t), so that offspring events may span several thousands years (Harte 2013); this practically hinders the analysis of cluster size for synthetic catalogues. However, we recall that, from the theoretical viewpoint, the temporal component of the ETAS model does not affect cluster size. Therefore, in order to prevent this inconvenience, instead of increasing the p parameter, we considered a truncated distribution f₂(t) on the support t < 0.7 yr, which corresponds to the 95 per cent quantile of the cluster lifetime for Northeastern Italy.

5. The space scattering law parameters: d = 0.0085 (degree), q = 2.25, γ = 0.62.

The synthetic earthquake catalogue, consisting of 143 568 events, spans about 500 yr and includes information about the links between each event and its direct descendants. The synthetic time-series were processed by the nearest-neighbour method (Zaliapin & Ben-Zion 2013) to identify clusters of events. This method allows partitioning earthquakes into background and clustered components, based on nearest-neighbour distances between earthquakes in the space–time–magnitude domain (Baiesi & Paczuski 2004; Zaliapin et al. 2008). In this study, the parameters necessary for the computation of nearest-neighbour distances, namely the b-value and the fractal dimension of epicentres d_f, were set according to the values estimated from the analysis of Northeastern Italy catalogue (Benali et al. 2020, and references therein): b-value = 0.9 and d_f = 1.1. The threshold of the nearest-neighbour distance, which separates the clustered and background components, was automatically set (log₁₀η₀ = −4.2) following a criterion based on a 1-D Gaussian mixture model with two modes, where the threshold is the maximum likelihood boundary between the two modes (Zaliapin & Ben-Zion 2013). The total Δ-productivity distribution was examined considering the clusters extracted from the synthetic earthquake catalogue.

Fig. 1 represents the total Δ-productivity distributions (Δ = 2) (a) for arbitrary initial events and (b) for main shocks, as identified by nearest-neighbour method applied to the synthetic ETAS(P) seismicity. In the first case, the distribution |$V_\Delta$| demonstrates a heavy tail (requiring the use of logarithmic values in the figure), which is also confirmed by the coefficient of variation |$\sigma (V_\Delta )/EV_\Delta \approx 10.8$|⁠. On the contrary, the distribution |$V_\Delta$| for main shocks (Fig. 1b) demonstrates the light tail behaviour; for comparison with the previous case |$\sigma (V_\Delta )/EV_\Delta \approx 0.7$|⁠.

Figure 1.

Clusters identified by nearest-neighbour method applied to the synthetic ETAS(P) seismicity: The |$V_\Delta$| distributions, Δ = 2, for (a) arbitrary initial events, m₀ ≥ 4 and (b) main shocks, m₀ ≥ 4.

Fig. 2 shows similar results based on the clusters originally defined by the branching structure of the synthetic ETAS(P) catalogue. The coefficients of variation are also consistent with previous results: (a) |$\sigma (V_\Delta )/EV_\Delta \approx 9.7$| for arbitrary initial events and (b) |$\sigma (V_\Delta )/EV_\Delta \approx 0.8$| for main shocks.

Figure 2.

Clusters in the synthetic ETAS(P) seismicity: The |$V_\Delta$| distributions, Δ = 2, for (a) arbitrary initial events, m₀ ≥ 4 and (b) main shocks, m₀ ≥ 4.

A good agreement of the empirical data Figs 1 and 2 demonstrates the robustness of conclusions about the tail behaviour of |$V_\Delta$| distribution.

4 CONCLUSION

The concept of the Δ-productivity |$\nu _\Delta$| is associated with a priori ideas about the structure of a cluster of seismic events in the form of a random branching tree. Such a representation is not the only possible one. Therefore, the Geometric distribution of |$\nu _\Delta$| can be considered as a useful alternative within the framework of the generalized ETAS(F) model. The other characteristic of the Δ-productivity, namely |$V_\Delta$| for main shocks, is more objective and stable.

For large class of the ETAS(F) models, we answered the question when the total Δ-productivity distribution has light tails for main shocks and heavy ones for arbitrary events. The real data demonstrate the possibilities of this kind.

Note the difficulties of substantiating the mentioned geometric laws for |$\nu _\Delta$| and aftershock’s |$V_\Delta$|⁠. They are associated both with the declustering of seismicity and with the effect of averaging the distributions of |$\nu _\Delta$| and |$V_\Delta$| when grouping data.

The first difficulty can be judged by experiments of Zaliapin & Ben-Zion (2013) with the nearest-neighbour method. This method adapted to the ETAS structure has 40 per cent error of incorrectly determining the parent of an event. The stability of the conclusions in this case may depend on the strict similarity in the real hierarchical structure of the seismic cluster.

The second difficulty is expressed in the fact that the mixture of geometric distributions retains the original property of monotonicity, but ceases to be geometric, that is loses the original G-type. This is unavoidable if the distribution parameter is dependent on location or the magnitude of the initial event (see eq. 6). Therefore, we can speak more confidently about the light tail of the |$V_\Delta$| distribution for the main shocks.

Finally, the theoretical result for ETAS(F) model shows that the types of distributions of |$\nu _\Delta$| and |$V_\Delta$| cannot be the same due to the branching structures of the clusters.

ACKNOWLEDGEMENTS

We are grateful to Prof J. Zhuang and another anonymous reviewer for their careful reading of the manuscript and comments that allowed us to significantly revise and improve the text.

DATA AVAILABILITY

The data (synthetic earthquake catalogue) underlying this paper will be shared on reasonable request to the corresponding author.

APPENDIX

A1 Proof of Statement 1

Let |$E\left(z^{\nu (m_0)}\mid m_0\right)=\varphi (z, \lambda (m_0))$| be the generating function of the ν(m₀) with the characteristic λ(m₀) = Eν(m₀). If {m_i: i = 1, …, ν(m₀)} are magnitudes of the direct descendants of the m₀ event, then the |$\nu _\Delta (m_0)$| has the same distribution as |$\sum _{i=1}^{\nu (m_0)} \varepsilon (m_i)$|⁠, where |$\varepsilon (m)=\mathbf {1}_{m\ge m_0-\Delta }$| is the logical 1-0 function. Since {m_i} are independent,

$$\begin{eqnarray*} \varphi _\Delta (z,\lambda (m_0)) &:=& E\left(z^{\nu _\Delta (m_0)}\mid m_0\right)= \sum _{n\ge 0} \left(Ez^{\varepsilon (m)}\right)^n P(\nu (m_0)=n) = \varphi \left(\pi _{m_0-\Delta }+\bar{\pi }_{m_0-\Delta }z,\lambda (m_0)\right)\quad , \end{eqnarray*}$$

(A1)

where |$\bar{\pi }_{m-\Delta }=P(m\ge m_0-\Delta \mid m_0\ge \Delta )$| and |$\pi _m=1-\bar{\pi }_m$|⁠.

Recall that φ(z, λ(m)) = (1 − p(m))/(1 − p(m)z) = [1 − λ(m)(z − 1)]⁻¹ for F = G and φ(z, λ(m)) = exp(λ(m)(z − 1)) for F = P. Substituting any of the functions φ(z, λ(m)) in the right-hand part of eq. (A1), we get the same function with the replacement of λ(m₀) by

$$\begin{eqnarray*} \lambda _\Delta (m_0)=\lambda (m_0) \bar{F_1}(m_0-\Delta )/\bar{F_1}(\Delta ) \quad . \end{eqnarray*}$$

(A2)

Assume that the right-hand part of eq. (A2) is independent of m₀ in the interval Δ ≤ m₀ < M. Then λ(m) and F₁(m) are differentiable functions simultaneously. Since |$\lambda _\Delta (m_0)={\rm const}$| for m₀ ≥ Δ, we have

$$\begin{eqnarray*} 0 = \left. \frac{d\lambda _\Delta (m_0)}{dm_0} \right|_{m_0=\Delta } = \left[ \dot{\lambda }(\Delta ) -\lambda (\Delta ) f_1(0) \right] /\bar{F}_1(\Delta ) \quad . \end{eqnarray*}$$

Here we used condition |$\bar{F}_1(0)=1$|⁠. Solving the equation |$\dot{\lambda }(\Delta )-\lambda (\Delta )f_1(0)=0$| in the interval 0 < Δ ≤ m₀, we get λ(m) = λe^am, m ∈ (0, M), where a = f₁(0). Substitute in eq. (A2) λ(m) = λe^am, |$\lambda _\Delta (m_0)=C(\Delta )$| and Δ = 0, to get |$\bar{F}_1(m)=C(0)\lambda ^{-1}e^{-am}$|⁠. Since |$\bar{F}_1(0)=1$|⁠, we have |$\bar{F}_1(m)=e^{-am}$|⁠.

A2 Proof of the statement 2

A2.1 General remarks.

We will need the following notation: N_M(m) is the total number of events of magnitude ≥M triggered by an event of magnitude m. In particular, the total Δ-productivity |$V_\Delta (m_0)$| of the initial event with magnitude m₀ is |$N_{M(m_0)}(m_0)$|⁠, M(m) = (m − Δ)₊, where a₊ = max(a, 0). The notation E(*∣m) denotes the conditional mean given m.

The spatio-temporal components of the ETAS model do not affect cluster size, since the law of generating new events depends only on the parent magnitude, and the problem of declustering does not arise in our theoretical analysis. Therefore, the spatio-temporal components in theoretical analysis can be any and the distribution of N_M(m) is the same for any x = (t, g, m).

By definition, if {m_i:  i = 1, …, ν(m₀)} are offspring magnitudes of m₀, then {N_M(m_i): i = 1, …, ν(m₀)} are independent for a given m₀. Moreover, they are identically distributed together with N_M(m₀) since the magnitudes m_i, i ≥ 0, have a common distribution f₁(m). Let us consider a new variable

$$\begin{eqnarray*} \tilde{N}_M(m_0) = \sum _{i=1}^{\nu (m_0)} \left\lbrace N_M(m_i) +\mathbf{1}_{m_i\ge M}\right\rbrace \end{eqnarray*}$$

(A3)

where the symbol |$\mathbf {1}_{a\ge b}$| denotes the logical function 1-0. Obviously, eq. (A3) is nothing more than counting events with magnitude ≥M in the cluster through first-generation events and its triggered events of magnitude ≥M. The component |$\mathbf {1}_{m_i\ge M}$| takes into account the ith event of the first generation if its value is greater than M. Therefore N_M(m₀) and |$\tilde{N}_M(m_0)$| are identically distributed, that is |$N_M(m_0)\overset{law}{=}\tilde{N}_M(m_0)$|⁠.

The moments |$\mathbf {EV_\Delta ^n}$|⁠.

We will use the notation:

$$\begin{equation*} \Lambda _n=\int \lambda ^n(m) f_1(m)dm\ , \quad \bar{\Lambda }_1=1-\Lambda _1, \quad \lambda (m)=p(m)/(1-p(m))\ \mbox{in the case}\ F=G\ , \end{equation*}$$

$$\begin{equation*} b_M^{(n)}(m)=E(N_M^n(m)\mid m),\quad b_M^{(n)}=EN_M^n(m),\quad \pi _M=P(m\ge M). \end{equation*}$$

Mean value: |$\mathbf {EV_\Delta }$|⁠.

Averaging eq. (A3) given m₀ = m, we have

$$\begin{eqnarray*} b_M^{(1)}(m)=\lambda (m)(b_M^{(1)}+\pi _M)\ . \end{eqnarray*}$$

After averaging over m, we obtain the equation for |$b_M^{(1)}$|⁠:

$$\begin{eqnarray*} b_M^{(1)}=\Lambda _1(b_M^{(1)}+\pi _M)\ . \end{eqnarray*}$$

For finite |$b_M^{(1)}$|⁠, this equality can be true if Λ₁ < 1. Hence

$$\begin{eqnarray*} b_M^{(1)}=\pi _M\cdot \Lambda _1/\bar{\Lambda }_1\ , b_M^{(1)}(m)=\lambda (m)\pi _M/\bar{\Lambda }_1\ . \end{eqnarray*}$$

(A4)

Since |$V_\Delta =N_{M(m)}(m)$|⁠, M(m) = (m − Δ)₊, we have

$$\begin{eqnarray*} \begin{array}{l}E(V_\Delta \mid m) =\lambda (m) \bar{F}_1(m-\Delta )/\bar{\Lambda }_1 \ , E(V_\Delta )\le \Lambda _1/\bar{\Lambda }_1\ , \\ E(V_\Delta \mid m\ge \Delta ) =E(V_\Delta \mathbf {1}_{m\ge \Delta })/\bar{F}_1(\Delta ) \le E(V_\Delta )/\bar{F}_1(\Delta ) \le (\Lambda _1/\bar{\Lambda }_1)/\bar{F}_1(\Delta )\ , \\ \end{array} \end{eqnarray*}$$

that is, |$V_\Delta$| has finite mean value in conditional and unconditional situations if Λ₁ < 1.

The case |$\mathbf {n\gt 1}$|⁠.

 

(a) Fixed/random initial |$\mathbf {m_0}$|

It is easy to see that for any λ(m) ≥ 0 there exists K ≤ ∞ such that Λ_k = Eλ^k(m) is finite for k ≤ K and is infinite for k > K. Let us show that there are finite constants B_k such that

$$\begin{eqnarray*} b_M^{(k)}\le B_k \,\,\,\, b_M^{(k)}(m)\lt \infty &\mbox{if {\rm and} only if} \Lambda _k\lt \infty \,\,\,\, \Lambda _1\lt 1\ . \end{eqnarray*}$$

(A5)

We have shown (see above) that eq. (A5) is true for k = 1. We will prove eq. (A3) by the method of induction on k. Assume that all values in eq. (A5) are finite for k ≤ n − 1. Putting |$\xi _i=N_M(m_i)+\mathbf {1}_{m_i\ge M}$| in eq. (A3), one has

$$\begin{eqnarray*} N_M^n(m)=\sum _{i_1\ne i_2\ne \dots \ne i_n}^{\nu (m)} \xi _{i_1}\dots \xi _{i_n} + R_M^{(n)} + \sum _{i=1}^{\nu (m)} \xi _i^n \ , \end{eqnarray*}$$

(A6)

where |$R_M^{(n)}$| represents the secondary component of |$N_M^{n}(m)$|⁠. The number of terms in the first sum is η_n = ν(m)(ν(m) − 1)…(ν(m) − n + 1). By the RT property of the F-distribution, Ez^ν(m) = φ(z, λ(m)) = ψ(λ(m)(z − 1)), where φ(z, λ(m)) has all derivatives at z = 1 because all moments of ν(m) are finite. Hence

$$\begin{eqnarray*} E\left[\eta _n \mid m\right] = \left.\left(\frac{\rm d}{{\rm d}z}\right)^n \varphi (z,\lambda (m)) \right|_{z=1} =a_n(F)\lambda ^n(m) \ , \end{eqnarray*}$$

(A7)

where |$a_n(F)=\left(\frac{\rm d}{{\rm d}z}\right)^n \left.\psi (z)\right|_{z=0}\lt \infty$|⁠. Therefore

$$\begin{eqnarray*} E\left( \sum _{i_1\ne i_2\ne \dots \ne i_n}^{\nu (m)} \xi _{i_1}\dots \xi _{i_n} \mid m \right) = a_n(F) \lambda ^n(m)[b_M^{(1)}+\pi _M]^n \ . \end{eqnarray*}$$

Similarly

$$\begin{eqnarray*} E\left( R_M^{(n)} \mid m \right) = \sum _\Omega c_n(\omega ) a_k(F) \lambda ^k(m) \tilde{b}_M^{(r_1)}\dots \tilde{b}_M^{(r_k)} \ , \end{eqnarray*}$$

(A8)

where

$$\begin{eqnarray*} \tilde{b}_M^{(r)} = E\left\lbrace N_M(m_i) + \mathbf {1}_{m_i\ge M} \right\rbrace ^r \ , \end{eqnarray*}$$

and

$$\begin{eqnarray*} c_n(\omega ) = \frac{n!}{r_1!\dots r_k!k!}\ , \Omega =\left\lbrace \omega =(k,r_1,\dots ,r_k): 1\lt k\lt n,\ r_i\ge 1,\ r_1+\dots +r_k=n \right\rbrace . \end{eqnarray*}$$

Since for any random variable ξ ≥ 0:

$$\begin{eqnarray*} E\xi ^n E\xi ^m \le E\xi ^{n+m} \ , \end{eqnarray*}$$

we get from eq. (A8):

$$\begin{eqnarray*} E\left( R_M^{(n)} \mid m \right) \le C_n (1 \vee \lambda (m))^{n-1} \underset{1\le k\lt n}{\max} \left( \tilde{b}_M^{(k)} \tilde{b}_M^{(n-k)} \right) \ , \end{eqnarray*}$$

(A9)

where C_n is some finite constant. The expectation of eq. (A6), given m, is

$$\begin{eqnarray*} b_M^{(n)}(m) = a_n(F) \lambda ^n(m) \left(\tilde{b}_M^{(1)}\right)^n + E\left( R_M^{(n)} \mid m \right) + \lambda (m) \left[ b_M^{(n)} +\delta _M^{(n)} \right] \ , \end{eqnarray*}$$

(A10)

where |$\delta _M^{(n)}=\tilde{b}_M^{(n)}-b_M^{(n)}$| is a secondary term of |$\tilde{b}_M^{(n)}$|⁠. By eq. (A8),

$$\begin{eqnarray*} \delta _M^{(n)} &\lt & E\left[\left( N_M(m)+1 \right)^n -N_M^{(n)}(m) \right] \lt n E\left[\left( N_M(m)+1 \right)^{n-1}\right] \lt n(2^{n-1} b_M^{(n-1)} +1) \ . \end{eqnarray*}$$

(A11)

Averaging eq. (A10) over m, we get

$$\begin{eqnarray*} b_M^{(n)}(1-\Lambda _1) = a_n(F) \Lambda _n \left(\pi _M / \bar{\Lambda }_1 \right)^n + ER_M^{(n)} + \Lambda _1 \delta _M^{(n)} \ . \end{eqnarray*}$$

(A12)

The relations (A11)–(A12) show that |$b_M^{(n)}$| is finite if and only if Λ_n < ∞. Similarly, eqs (A10)–(A11) show that |$b_M^{(n)}(m)\lt \infty$| if and only if |$b_M^{(n)}\lt \infty$|⁠. Hence, at the nth step of induction, these three quantities are finite simultaneously, that is eq. (A5) holds. By setting M = (m − Δ)₊, we get |$E(V_\Delta ^n \mid m)=b_{(m-\Delta )_+}^{(n)}(m)\lt \infty$| only when Λ_n < ∞.It remains to consider |$E(V_\Delta ^n)$|⁠. By eqs (A10) and(A12),

$$\begin{eqnarray*} b_M^{(n)}(1-\Lambda _1) = a_n(F) \Lambda _n \left( \pi _M /\bar{\Lambda }_1 \right)^n + ER_M^{(n)} +\Lambda _1 \delta _M^{(n)} \ , \\ E(V_\Delta ^n) = Eb_{(m-\Delta )_+}^{(n)}(m) \ge E\lambda (m) b_{(m-\Delta )_+}^{(n)} \ge E\lambda (m) a_n(F) \Lambda _n \pi _{(m-\Delta )_+}^{n} / \bar{\Lambda }_1^{n+1} \ . \end{eqnarray*}$$

Therefore |$E(V_\Delta ^n)=\infty$| if Λ_n is infinite.Let Λ_n < ∞. Using (A10) again, we have

$$\begin{eqnarray*} E(V_\Delta ^n) \le a_n(F) \Lambda _n (\bar{\Lambda }_1)^{-n} + E\left(R_{M(m)}^{(n)}\mid m\right) + E\lambda (m) \left[b_{M(m)}^{(n)}+\delta _{M(m)}^{(n)}\right] \ . \end{eqnarray*}$$

Due to eq. (A9), we have |$E\left(R_{M(m)}^{(n)}\mid m\right)\le \tilde{C}_n~\max(1,\Lambda _1,\dots ,\Lambda _{n-1})$|⁠; by eq. (A9), |$E\lambda (m)\delta _{M(m)}^{(n)} \le D_n\Lambda _1$|⁠. Finally, |$E\lambda (m)b_{M(m)}^{(n)}\le B_n\Lambda _1$|⁠, because we have proved (A3) for index n. Hence, |$E(V_\Delta ^n)\lt \infty$| if Λ_n < ∞. The proof is complete.

(b) Fixed dominant |$\mathbf {m_0}$|

Let us consider the ETAS cluster with a dominant initial magnitude m₀. Let |$\hat{N}_{0M}(m)$| be the number of events of magnitude >M in the cluster with the initial event m ≤ m₀ in which all events have magnitude <m₀. In such cluster, direct descendants of any initial magnitude m are realized under the condition Ω(m) = {m_i < m₀: i = 1, …, ν(m₀)}.

In term of the generating function Ez^ν(m) = φ(z, λ(m)), the probability of Ω(m) is

$$\begin{eqnarray*} P(\Omega (m)) = \varphi (F_1(m_0),\lambda (m)) = \psi (-\lambda (m)\bar{F}_1(m_0)). \end{eqnarray*}$$

(A13)

Therefore, the conditional probability density of the direct descendants is

$$\begin{eqnarray*} P\left\lbrace m_i=\mu _i: i=1,\dots ,n;\ \nu (m_0)=n \mid \Omega (m) \right\rbrace = \hat{f}_{01}(\mu _1) \dots \hat{f}_{01}(\mu _n)P(\hat{\nu }(m)=n). \end{eqnarray*}$$

where

$$\begin{eqnarray*} && \hat{f}_{01}(m) = f_1(m)/F_1(m_0) \ , 0\le m \le m_0 \ ,\ \ \ F_1(m)=\int _0^m f_1(u){\rm d}u \end{eqnarray*}$$

(A14)

$$\begin{eqnarray*} && P(\hat{\nu }(m)=n) = P(\nu (m)=n \mid \Omega (m)) = P(\nu (m)=n) F_1^n(m_0) / \psi (-\lambda (m) \bar{F}_1(m_0)) \end{eqnarray*}$$

(A15)

Substituting Geometric (eq. 1) and Poisson (eq. 2) distributions here, we get the distributions of the same type

$$\begin{eqnarray*} P(\nu (m_0)=n \mid \Omega (m)) = \left\lbrace \begin{array}{cc}\hat{\lambda }_0^n(m) e^{-\hat{\lambda }_0(m)} / n! & [P{\text {-model}}] \\\hat{p}_0^n(m) (1-\hat{p}_0(m)) & [G{\text {-model}}] \\\end{array}\right. , \end{eqnarray*}$$

(A16)

where

$$\begin{eqnarray*} \hat{p}_0(m) &=& p(m) F_1(m_0) = \lambda (m) F_1(m_0) / (1+\lambda (m)) [F=G] \nonumber \\ \hat{\lambda }_0(m) &=& \lambda (m) F_1(m_0) \left\lbrace \begin{array}{cc}1 & [F=P] \\ (1+\lambda (m)\bar{F}_1(m_0))^{-1} & [F=G] \\ \end{array}\right. \end{eqnarray*}$$

(A17)

Therefore, given the dominance of the initial magnitude m₀, we are again dealing with an ETAS model in which the magnitude m is transformed into |$\hat{m}$| with a distribution (eq. A14), and the productivity ν(m) into |$\hat{\nu }(\hat{m})$| with the same distribution (Poisson/Geometric), but with new parameters |$\hat{\lambda }_0(\hat{m})$| as in eq. (A17). Therefore, as above, we will have the following relations:

$$\begin{eqnarray*} \hat{N}_{0M}(\hat{m}) \overset{law}{=} \sum _{i=1}^{\hat{\nu }(\hat{m})} \left( \hat{N}_{0M}(\hat{m}_i) + \mathbf {1}_{\hat{m}_i\ge M} \right) \ , \end{eqnarray*}$$

(A18)

where |$\left\lbrace \hat{N}_{0M}(\hat{m}_i): i=1,\dots , \hat{\nu }(\hat{m})\right\rbrace$| are independent and identical distributed given |$\hat{\nu }(\hat{m})$| and |$\left\lbrace \hat{m}_i: i=1,\dots , \hat{\nu }(\hat{m}) \right\rbrace$|⁠.

As a result

$$\begin{eqnarray*} V_\Delta (m_0) = \hat{N}_{0(m_0-\Delta )_+}(m_0) \ ,\ \ \ m_0\ge \Delta \ . \end{eqnarray*}$$

(A19)

Since m₀ is a parameter, in the random case of initial event, m₀ has distribution |$F_1(m)/\bar{F}_1(\Delta )$|⁠.As we can see, eq. (A18) is a special case of eq. (A3), in which the magnitude range is finite. By eq. (A17), |$\hat{\lambda }_0(\hat{m})\le \lambda (m_0)$| and therefore |$\hat{\lambda }_0(\hat{m})$| has all moments. Consequently, the total Δ-productivity will have all the moments at any fixed m₀.

(c) Random dominant |$\mathbf {m_0}$|

We need additional notation:

$$\begin{eqnarray*} \begin{array}{l}\hat{b}_{0M}^{(n)}(\hat{m})=E\left[ \hat{N}_{0M}^{n}(\hat{m}) \mid \hat{m} \right] \ , \hat{b}_{0M}^{(n)} = E\hat{N}_{0M}^{n}(\hat{m}) \ , \\ {}\\ \tilde{b}_{0M}^{(n)} = E\left( \hat{N}_{0M}(\hat{m}) +\mathbf {1}_{\hat{m}\ge M} \right)^n \ , \\ {}\\ \hat{\pi }_{0M} = P(\hat{m}\ge M) = [F_1(m_0)-F_1(M)]_+/F_1(m_0) \ , \\ {}\\ \hat{\Lambda }_{0n} = E\hat{\lambda }_0^n(\hat{m}) \ . \end{array} \end{eqnarray*}$$

(A20)

Using eq. (A17) and non-decreasing function λ(m), we have for the P-model

$$\begin{eqnarray*} \hat{\Lambda }_{0n} = \int _0^{m_0} \lambda ^n(x) F_1^{n-1}(m_0) f_1(x) dx \le \hat{\lambda }_0^{n-1}(m_0) \Lambda _1 \ . \end{eqnarray*}$$

(A21)

The same is true for the G-model. Really, by eq. (A17) we have the following representation:

$$\begin{eqnarray*} \hat{\Lambda }_{0n} = \hat{\lambda }_0^{n-1}(m_0) \int _0^{m_0} [\psi ^{n-1}(\lambda (x))/(1+\lambda (m) \bar{F}(m_0))] \lambda (x) f_1(x) dx \ , \nonumber \end{eqnarray*}$$

where

$$\begin{eqnarray*} \psi (m) = (\lambda ^{-1}(m_0) + \bar{F}_1(m_0)) / (\lambda ^{-1}(m) + \bar{F}_1(m_0)) \ . \nonumber \end{eqnarray*}$$

Since λ(m) is a non-decreasing function, ψ(m) ≤ 1 on the interval (0, m₀). Therefore, the integrand in square brackets is less than 1, which proves eq. (A21).As above, the relation eq. (A18) gives,

$$\begin{eqnarray*} \hat{b}_{0M}^{(1)} &=& \hat{\Lambda }_{01}\hat{\pi }_{0M}/(1-\hat{\Lambda }_{01}) \le \Lambda _1\hat{\pi }_{0M}/(1-\Lambda _1) \ , \\ EV_\Delta & \le & \Lambda _1/(1-\Lambda _1) \ . \nonumber \end{eqnarray*}$$

(A22)

Setting |$\xi _i=\hat{N}_{0M}(\hat{m}_i)+\mathbf {1}_{\hat{m}_i\ge M}$| in eq. (A18) and using eq. (A6), we get an analogue of eq. (A10):

$$\begin{eqnarray*} \hat{b}_{0M}^{(n)}(\hat{m}) = \sum _{k=1}^n \sum _{\mathbf {r^{(k)}}} c_k(r_1,\dots ,r_k) \hat{\lambda }_0^k(\hat{m}) \tilde{b}_{0M}^{(r_1)}\dots \tilde{b}_{0M}^{(r_k)} \ , \end{eqnarray*}$$

(A23)

where the summation is over all integer vectors |$\mathbf {r^{(k)}}=(r_1,\dots ,r_k)$|⁠: r₁ + … + r_k = n, r_i ≥ 1. Averaging eq. (A23) over |$\hat{m}$| we will also have

$$\begin{eqnarray*} \hat{b}_{0M}^{(n)} = \sum _{k=2}^n \sum _{\mathbf {r^{(k)}}} c_k(r_1,\dots ,r_k) \hat{\Lambda }_{0k} \tilde{b}_{0M}^{(r_1)}\dots \tilde{b}_{0M}^{(r_k)} + \hat{\Lambda }_{01} \tilde{b}_{0M}^{(n)} \ . \end{eqnarray*}$$

(A24)

Let’s extract the term |$\hat{\Lambda }_{01} \hat{b}_{0M}^{(n)}$| from |$\hat{\Lambda }_{01} \tilde{b}_{0M}^{(n)}$| and move it to the left-hand part of eq. (A24); then apply eq. (A21) to |$\hat{\Lambda }_{0k}$| and take into account the relation |$\hat{\Lambda }_{01}/(1-\hat{\Lambda }_{01}) \le \Lambda _1/(1-\Lambda _1)$|⁠. Then multiplying eq. (A24) by |$\hat{\lambda }_0(m_0)$|⁠, we get

$$\begin{eqnarray*} \hat{\lambda }_0(m_0)\hat{b}_{0M}^{(n)} &\le & \displaystyle \frac{\Lambda _1}{1-\Lambda _1} \sum _{k=2}^n \sum _{\mathbf {r^{(k)}}} \left\lbrace c_k(r_1,\dots ,r_k) (\hat{\lambda }_0(m_0)\tilde{b}_{0M}^{(r_1)}) \dots (\hat{\lambda }_0(m_0)\tilde{b}_{0M}^{(r_k)}) + \hat{\lambda }_0(m_0) (\tilde{b}_{0M}^{(n)}-\hat{b}_{0M}^{(n)}) \right\rbrace . \end{eqnarray*}$$

(A25)

Since |$\tilde{b}_{0M}^{(n)} = E\left( \hat{N}_{0M}(\hat{m}_i) +\mathbf {1}_{\hat{m}_i\ge M}\right)^n$|⁠, the elements of the right-hand part of eq. (A25) can be obtained by summing and multiplying elements of the form: |$\hat{\lambda }_0(m_0)\hat{\pi }_{0M}$| and

$$\begin{eqnarray*} D_{k,\varepsilon } =\hat{\lambda }_0(m_0) E\hat{N}_{0M}^k(\hat{m}) (\mathbf {1}_{\hat{m}\ge M})^\varepsilon\, k=1,\dots ,n-1;\ \ \varepsilon =0,1 \ , \end{eqnarray*}$$

(A26)

where |$D_{k,\varepsilon } \le \hat{\lambda }_0(m_0) \hat{b}_{0M}^{(k)}$|⁠.We have assumed that |$\lambda (m)\displaystyle \int _0^\Delta f_1(m-x)dx \le C$|⁠. Hence

$$\begin{eqnarray*} \hat{\lambda }_0(m_0) \hat{\pi }_{0}(m_0-\Delta ) = \lambda (m_0) F_1(m_0) \left[ F_1(m_0)-F_1(m_0-\Delta ) \right] /F_1(m_0) \le C \end{eqnarray*}$$

and by virtue of eq. (A22)

$$\begin{eqnarray*} \hat{\lambda }_0(m_0) \hat{b}_{0M(m_0)}^{(1)} \le C_1 \, M(m_0) =(m_0-\Delta )_+ \ . \end{eqnarray*}$$

Using eq. (A25) and doing induction over n, we can conclude that |$\hat{\lambda }_0(m_0) \hat{b} _{0M(m_0)}^{(n)}\le C_n$| for any n. But then, due to eq. (A19), |$EV_\Delta ^n = E\hat{b}_{0M(m_0)}^{(n)}(m_0) \le K_n\lt \infty$| for any n.

The proof is complete.

A3 Proof of the statement 3.1

Let ξ be a random variable with a geometric distribution. Then

$$\begin{eqnarray*} E\xi (\xi -1)/(E\xi )^2=2\ . \end{eqnarray*}$$

(A27)

Consider |$\xi =V_\Delta (m_0)$|⁠, where m₀ is the initial dominant magnitude in the ETAS(G) cluster. We will show that the relation (A27) is possible for at most one parameter Δ ∈ (0, m₀). For simplicity of notation, we will temporarily proceed from eq. (A3). Recall that

$$\begin{eqnarray*} \begin{array}{ll}\Lambda _n = \displaystyle \int \lambda ^n(m) f_1(m) dm \ , & \pi _M = P(m\ge M)\ , \\ b_M^{(n)}(m) = E(N_M^n(m)\mid m) \ , & b_M^{(n)} =EN_M^n(m)\ . \end{array} \end{eqnarray*}$$

Above we show that

$$\begin{eqnarray*} \begin{array}{lll}b_M^{(1)}(m) = \lambda (m) \pi _M / \bar{\Lambda }_1 \ , & b_M^{(1)} = \pi _M\Lambda _1 / \bar{\Lambda }_1 \ , & b_M^{(1)} + \pi _M = \pi _M/\bar{\Lambda }_1 \ . \end{array} \end{eqnarray*}$$

(A28)

According to eq. (A10) for the ETAS(G) model, we have

$$\begin{eqnarray*} b_M^{(2)}(m) = 2(\lambda (m) \pi _M/\bar{\Lambda }_1)^2 +\lambda (m)(b_M^{(2)}+2b_M^{(1)}\pi _M+\pi _M^2) \ . \end{eqnarray*}$$

(A29)

By eq. (A28),

$$\begin{eqnarray*} 2b_M^{(1)}\pi _M+\pi _M^2 = (1+\Lambda _1)\pi _M^2/\bar{\Lambda }_1 \ . \end{eqnarray*}$$

(A30)

Averaging eq. (A29) over m, we can find |$b_M^{(2)}$|⁠; using eq. (A30), we will have

$$\begin{eqnarray*} b_M^{(2)} = 2\Lambda _2\pi _M^2/\bar{\Lambda }_1^3 + \Lambda _1(1+\Lambda _1)\pi _M^2/\bar{\Lambda }_1^2 \ . \end{eqnarray*}$$

(A31)

If we substitute eq. (A31) into eq. (A29), we get

$$\begin{eqnarray*} b_M^{(2)}(m) = 2(b_M^{(1)})^2 + \lambda (m) (\pi _M/\bar{\Lambda }_1)^2 (2\Lambda _2+1-\Lambda _1^2) /\bar{\Lambda }_1 \ . \end{eqnarray*}$$

As a result,

$$\begin{eqnarray*} (b_M^{(2)}(m)-b_M^{(1)}(m)) / (b_M^{(1)}(m))^2 = 2+\lambda ^{-1}(m) \left[ (2\Lambda _2+1-\Lambda _1^2) / \bar{\Lambda }_1 - \bar{\Lambda }_1/\pi _M \right] \ . \end{eqnarray*}$$

(A32)

The condition (A27) is satisfied only if

$$\begin{eqnarray*} \pi _M = (1-\Lambda _1)^2 / (2\Lambda _2+1-\Lambda _1^2) \ , \end{eqnarray*}$$

(A33)

where |$\Lambda _2-\Lambda _1^2 = \sigma ^2(\lambda )\gt 0$|⁠.The eq. (A33) with respect to M has unique solution if π_M is strictly monotone.Now we can return to the ETAS(G) cluster with dominant initial magnitude m₀. In this case the main equation is (A18). It is identical to eq. (A10), but requires the following replacement of the characteristics used (see A14 and A17):

$$\begin{eqnarray*} f_1(m) & \Longrightarrow & f_1(m)/F_1(m_0), 0\le m\le m_0 \ ,\nonumber \\ \lambda (m) & \Longrightarrow & \lambda (m)F_1(m_0)/(1+\lambda (m) \bar{F}_1(m_0)) \ , \nonumber \\ \pi _M & \Longrightarrow & [F_1(m_0)-F_1(M)]_+ / F_1(m_0) \ , \end{eqnarray*}$$

(A34)

$$\begin{eqnarray*} \Lambda _n & \Longrightarrow & \int _0^{m_0} \lambda ^n(m) F_1^{n-1}(m_0) (1+\lambda (m)\bar{F}_1(m_0))^{-n} f_1(m) {\rm d}m \ . \end{eqnarray*}$$

(A35)

Since |$V_\Delta (m_0)=N_{m_0-\Delta }(m_0)$|⁠, we have to consider eq. (A33) given the substitutions (A34) and (A35) and M = m₀ − Δ. Given f₁(m) > 0, the function (A34) will be strictly monotonic with respect to Δ. Hence, the relation (A27) is possible for at most one value of Δ = Δ(m₀).

The proof is complete.

A4 Proof of the statement 3.2

Let ν(m) and V(m) be the sizes of clusters of direct descendants and all descendants, respectively, associated to an initial event with magnitude m. Let ϕ_V(z, λ(m)) = Ez^V(m) and |$\phi _V(z)=\displaystyle \int _0^\infty Ez^{V(m)} f_1(m) dm$| be the conditional and unconditional generating functions of V(m). Due to the RT property of ν(m) distribution, we have

$$\begin{eqnarray*} Ez^{\nu (m)} = \psi (\lambda (m)(z-1)), \quad \quad \psi (0)=\psi ^{\prime }(0)=1, \quad \quad \psi ^{\prime \prime }(0)\lt \infty . \end{eqnarray*}$$

Then by eq. (A3)

$$\begin{eqnarray*} \phi _V(z, \lambda (m)) = \psi (\lambda (m)(z\phi _V(z)-1)), \end{eqnarray*}$$

and after integration over m we get

$$\begin{eqnarray*} \phi _V(z) = \int _0^\infty \psi (\lambda (m)(z\phi _V(z)-1)) f_1(m) dm. \end{eqnarray*}$$

(A36)

Since it is assumed that ∫λ²(m)f₁(m)dm < ∞ and Eν² < ∞, it follows EV² < ∞ from Statement 2. This means that |$\phi _V^{\prime \prime }(1)\lt \infty$|⁠. Set |$\mu =EV=\phi _V^{\prime }(1)$|⁠, |$\Lambda _k=\int _0^\infty \lambda ^k(m) f_1(m) dm$|⁠. Differentiating eq. (A36) at point 1, we get

$$\begin{eqnarray*} \mu = \Lambda _1 (1+\mu ) \quad \mbox{or} \quad \mu = \Lambda _1/(1-\Lambda _1). \end{eqnarray*}$$

(A37)

Differentiating eq. (A36) twice at point 1, we get

$$\begin{eqnarray*} \phi _V^{\prime \prime }(1) = \psi ^{\prime \prime }(0) \Lambda _2(1+\mu )^2 + \Lambda _1(2\mu +\phi _V^{\prime \prime }(1)). \end{eqnarray*}$$

Given eq. (A37), we have

$$\begin{eqnarray*} \phi _V^{\prime \prime }(1) = \psi ^{\prime \prime }(0) \Lambda _2(1-\Lambda _1)^{-3} +2\mu ^2. \end{eqnarray*}$$

(A38)

Assume that the distribution types of ν and V are the same. Then ϕ_V(z) = ψ(μ(z − 1)) and |$\phi _V^{\prime \prime }(1)=\psi ^{\prime \prime }(0)\mu ^2$|⁠. Combining this equality with eq. (A38) we have

$$\begin{eqnarray*} \psi ^{\prime \prime }(0) = \psi ^{\prime \prime }(0) \Lambda _2\Lambda _1^{-2}(1-\Lambda _1)^{-1} +2. \end{eqnarray*}$$

(A39)

By definition, all terms in this equality are non-negative, |$\Lambda _2\Lambda _1^{-2}\ge 1$|⁠, and (1 − Λ₁)⁻¹ > 1. Therefore, eq. (A39) is contradictory.