# A positive answer on the existence of correlations between positive
earthquake magnitude differences

Eugenio Lippiello  Department of Mathematics and Physics, University of
Campania “Luigi Vanvitelli”, Caserta, Italy  Lucilla de Arcangelis  Department
of Mathematics and Physics, University of Campania “Luigi Vanvitelli”,
Caserta, Italy  Cataldo Godano  Department of Mathematics and Physics,
University of Campania “Luigi Vanvitelli”, Caserta, Italy

###### Abstract

The identification of patterns in space, time, and magnitude, which could
potentially encode the subsequent earthquake magnitude, represents a
significant challenge in earthquake forecasting. A pivotal aspect of this
endeavor involves the search for correlations between earthquake magnitudes, a
task greatly hindered by the incompleteness of instrumental catalogs. A novel
strategy to address this challenge is provided by the groundbreaking
observation by Van der Elst (Journal of Geophysical Research: Solid Earth,
2021), that positive magnitude differences, under certain conditions, remain
unaffected by incompleteness. In this letter we adopt this strategy which
provides a clear and unambiguous proof regarding the existence of correlations
between subsequent positive magnitude differences. Our results are consistent
with a time-dependent $b$-value in the Gutenberg-Richter law, significantly
enhancing existing models for seismic forecasting.

Seismic catalogs suffer from significant incompleteness. Small earthquakes
often go unreported due to their occurrence at distances beyond the detection
range of monitoring stations [1, 2, 3], or they may be overshadowed by the
coda-wave of preceding larger earthquakes [4, 5, 6, 7, 8, 9, 10, 11, 12]. This
latter mechanism, in particular, renders the catalog inherently incomplete [7,
13] and has left fundamental questions about seismic activity unanswered. One
such crucial inquiry involves seismic forecasting and the potential dependence
of earthquake magnitudes on the organization of seismic events in time, space,
and magnitude preceding its occurrence. Studies have convincingly demonstrated
correlations between earthquake magnitudes, with these correlations depending
on both temporal and spatial distances [14, 15, 16, 17, 18, 19, 20, 21, 22].
However, the challenge lies in discerning whether these observed correlations
are genuine features or spurious artifacts of catalog incompleteness [23, 22,
24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34]. Consequently, the most commonly
adopted statistical forecasting models, such as the ETAS model [35, 36, 37,
38], assume earthquake magnitudes to be independent and identically
distributed (i.i.d.) variables. Within this framework, forecasting is severely
limited, relying solely on the information that large earthquakes have a
higher probability of occurring when and where the seismic rate is higher.
Modifications of the ETAS model, incorporating magnitude correlations have
been proposed in the literature [14, 39, 29, 30]

Recently, Van der Elst [40] introduced a novel and remarkable tool designed to
address issues associated with catalog incompleteness in a simple yet highly
effective manner. This tool is grounded in the observation that, under
specific conditions, positive magnitude differences remain largely unaffected
by incompleteness [40, 41]. Leveraging this approach, researchers have been
able to derive the true magnitude distribution, the true occurrence rate [42],
and the correct parameters of the ETAS model [43], even when working with
incomplete catalogs. In this study, we employ this tool to unequivocally
demonstrate the presence of genuine correlations between earthquake
magnitudes.

A seismic catalog is a multi-variate dataset ${\cal C}$ consisting of $N$
elements, where the $i$-th element is the vector ${m_{i},t_{i},\vec{r}_{i}}$
representing the earthquake magnitude, occurrence time, and epicentral
coordinates, respectively, of the $i$-th earthquake. Specifically, we examine
the relocated catalog for Southern California (SC) [44], spanning from January
1981 to March 2022, encompassing $N=800499$ earthquakes with $m\geq-1$.
Additionally, we consider the relocated catalog for Northern California (NC)
[45], covering the period from January 1984 to December 2021, including
$N=879547$ earthquakes with $m\geq-0.8$. From the original catalog, we also
construct a catalog of magnitude differences, denoted as ${\cal D}$, where
each element $i$ is represented by the vector ${\delta
m_{i},t_{i},\vec{r}_{i}}$, with $\delta m_{i}=m_{i+1}-m_{i}$ representing the
magnitude difference between two successive earthquakes in the original
catalog.

To investigate correlations between earthquake magnitudes, we extend the
method developed in [16]. We define the conditional probability as:

$P(\delta m_{i}>M|\delta y_{i}<Y)\equiv\frac{N(M,Y)}{N(Y)}.$ (1)

Here, $\delta y_{i}$ denotes either the epicentral distance $\delta
r_{i}=|\vec{r}_{i+1}-\vec{r}_{i}|$ or the temporal distance $\delta
t_{i}=t_{i+1}-t_{i}$, depending on the context. In Eq.(1), $N(M,Y)$ represents
the number of pairs of subsequent events where both $\delta m_{i}>M$ and
$\delta y_{i}<Y$, while $N(Y)$ is the total number of pairs with $\delta
y_{i}<Y$. Starting from the initial catalog $\mathcal{C}$, we generate a
catalog with reshuffled magnitude differences $\mathcal{D}_{\text{ran}}$,
where the $i$-th element is the vector ${\delta m^{*}_{i},t_{i},\vec{r}_{i}}$,
with $\delta m^{*}_{i}=m_{k}-m_{i}$, and $k$ represents the index of an
earthquake randomly chosen within the catalog. Specifically, the catalog
$\mathcal{D}_{\text{ran}}$ retains the same occurrence time and epicentral
coordinates as catalog $\mathcal{D}$, but the magnitude differences are
randomized. From a single instrumental catalog $\mathcal{C}$, we can derive
multiple reshuffled catalogs $\mathcal{D}_{\text{ran}}$. In each reshuffled
catalog, we evaluate the quantity $P(\delta m^{*}_{i}>M|\Delta y_{i}<Y)$,
which exhibits a different value in each reshuffled catalog
$\mathcal{D}_{\text{ran}}$. This quantity follows a Gaussian distribution with
a mean value $\overline{P(M,Y)}$ and a standard deviation $\Sigma(M,Y)$.

Figure 1: (a) The quantity $E\left[\Delta P(M,Y)\right]$ for SC, where $\Delta
P(M,Y)$ is defined in Eq.(2), with $\delta y_{i}=\delta r_{i}$ and $Y=10$ km.
Different curves correspond to different thresholds $M_{th}$. Error bars
represent $2\Sigma(M,Y)$. (b) The magnitude distribution $P(m_{i})$ and the
positive magnitude difference distribution $P(q_{i})$ for SC and NC. The
dashed orange line represents the GR law for SC with $\beta=2.23$.

The key quantity of the study is defined as:

$\Delta P(M,Y)=P(\delta m_{i}>M|\delta y_{i}<Y)-P(\delta m^{*}_{i}>M|\Delta
y_{i}<Y).$ (2)

We then calculate the mean value $E\left[\Delta P(M,Y)\right]$, obtained by
averaging over 1000 realizations of the reshuffled catalog
$\mathcal{D}_{\text{ran}}$. $E\left[\Delta P(M,Y)\right]$ for SC is plotted in
Fig.1a for the case where $\delta y_{i}=\delta r_{i}$ and $Y=10$ km. The error
bars represent twice the standard deviation $\Sigma(\delta M,Y)$. Values such
as $|E\left[\Delta P(M,Y)\right]|>2\Sigma(M,Y)$ indicate that it is highly
improbable for a catalog with uncorrelated magnitudes, such as
$\mathcal{D}_{\text{ran}}$, to produce the same probability of observing a
given magnitude difference $P(\delta m_{i}>M|\Delta y_{i}<Y)$ of the real
catalog. Therefore, Fig.1a demonstrates that magnitudes in the instrumental
catalog $\mathcal{C}$ are correlated. In Fig.1a, we also examine the influence
of a magnitude threshold $M_{\text{th}}$ on $E\left[\Delta P(M,Y)\right]$. By
considering only earthquakes with magnitudes larger than $M_{\text{th}}$ and
increasing $M_{\text{th}}$, we observe a decrease in $E\left[\Delta
P(M,Y)\right]$. This decrease can be interpreted as a signature that magnitude
correlations are primarily a spurious effect of incompleteness, as increasing
$M_{\text{th}}$ focuses on a catalog less affected by incompleteness.
Nevertheless, the observation that even for $M_{\text{th}}=3$, which exceeds
the completeness magnitude $M_{c}$ estimated from the magnitude distribution
(Fig.1b), $|E\left[\Delta P(M,Y)\right]|>2\Sigma(M,Y)$ suggests the
persistence of genuine magnitude correlations beyond incompleteness. This
scenario finds further support in the resilience of $E\left[\Delta
P(M,Y)\right]$ to variations in the quality of the seismic network [22] and is
reinforced by laboratory experiments on rock fracture [33]. These experimental
findings demonstrate the existence of magnitude clustering irrespective of
loading protocols, rock types, and observable magnitude ranges.

Next, we will review these results in light of recent findings concerning
positive magnitude difference statistics. We assume as null hypothesis that
magnitudes are i.i.d variables which obey the Gutenberg-Richter (GR) law,
$p(m)\propto\exp(-\beta(m-M_{L}))$ for magnitudes larger than a lower limit
$M_{L}$. Although a precise estimate of $M_{L}$ is unavailable, it is
reasonable to assume that $M_{L}$ takes on a very negative value [25].
Additionally, we assume that not all magnitudes are reported in the catalog
due to detection issues. Therefore, the observed magnitude distribution is
given by [41]

$p\left(m_{i}\right)=\beta
e^{-\beta(m_{i}-M_{L})}\Phi\left(m_{i}-M_{c}\right),$ (3)

where $M_{c}>M_{L}$ is the completeness magnitude and $\Phi(x)$ is a detection
function which is a monotonic increasing function of $x$ ranging from
$\Phi(x)=0$ for $x\lesssim-2\sigma$ to $\Phi(x)=1$ for $x\gtrsim 2\sigma$.
Essentially, earthquakes with magnitudes $m_{i}<M_{c}-2\sigma$ are undetected,
while those with $m_{i}>M_{c}+2\sigma$ are always detected. A realistic
functional form for $\Phi(x)$ is
$\Phi(x)=\left(1+\text{Erf}(x/\sigma)\right)/2$, with typical estimates
placing $\sigma$ in the range $\sigma\in(0.2,0.3)$ [46, 40, 41]. Eq.(3) is
consistent with the magnitude distribution of instrumental catalogs (Fig.1b),
which conforms to the GR law only for magnitudes $m_{i}>M_{c}$, determined
using the maximum curvature method [47], with $M_{c}=2.6$ for SC and
$M_{c}=2.7$ for NC (Fig.1b).

Using Eq.(3) in Eq.(1), we obtain

$\displaystyle P(\delta m_{i}>M|\Delta
y_{i}<Y)=\beta^{2}\int_{m_{L}}^{\infty}dm_{i}\int_{m_{i}+M}^{\infty}dm_{i+1}$
(4) $\displaystyle
e^{-\beta(m_{i+1}+m_{i}-2m_{L})}\Phi\left(m_{i}-M_{c}\right)\Phi\left(m_{i+1}-M_{c}|m_{i}\right).$

In the above equation, we explicitly use the notation
$\Phi\left(m_{i+1}-M_{c}|m_{i}\right)$ to specify that the detection function
must be evaluated under conditions where the previous earthquake $m_{i}$ has
been identified and reported in the catalog. It is worth noticing that it is
exactly this conditioned detection function that introduces correlations
between the magnitudes $m_{i}$ and $m_{i+1}$. Indeed, the same expression
Eq.(4) holds for the probability $P(\delta m^{*}_{i}>M|\Delta y_{i}<Y)$, with
$m_{k}$ instead of $m_{i+1}$. This makes a fundamental difference since, by
construction, $k$ is a random index, and therefore
$\Phi\left(m_{k}-M_{c}|m_{i}\right)=\Phi\left(m_{k}-M_{c}\right)$,
independently of $m_{i}$.

Accordingly, $\Delta P(M,Y)$ can be rewritten as

$\displaystyle\Delta
P(M|Y)=\beta^{2}\int_{m_{L}}^{\infty}dm_{i}\int_{m_{i}+M}^{\infty}dm_{i+1}e^{-\beta(m_{i+1}+m_{i}-2m_{L})}$
(5)
$\displaystyle\Phi\left(m_{i}-M_{c}\right)\left[\Phi\left(m_{j}-M_{c}|m_{i}\right)-\Phi\left(m_{j}-M_{c}\right)\right].$

The origin of magnitude correlations in the null hypothesis scenario
originates from the difference of the two detection functions in the term
between square brackets in Eq.(5). It is evident that, if $M>0$, $m_{j}>m_{i}$
and $\Phi\left(m_{j}-M_{c}|m_{i}\right)$ is on average larger than
$\Phi\left(m_{j}-M_{c}\right)$. Indeed, if a previous earthquake has been
detected, the probability to detect a subsequent larger one is larger compared
to the case where no information about previously detected earthquakes is
available. At the same time, imposing $m_{j}\geq M_{\text{th}}$ and by
considering increasing values of $M_{\text{th}}$, both detection functions
$\Phi\left(m_{j}-M_{c}|m_{i}\right)$ and $\Phi\left(m_{j}-M_{c}\right)$
approach $1$, and $\Delta P(M,Y)$ approaches zero. Accordingly, the behavior
of $E\left[\Delta P(M,Y)\right]$ in Fig.1a cannot exclude the null hypothesis
scenario where magnitudes are i.i.d variables in the presence of detection
problems.

We now explore the prediction of the null hypothesis on the magnitude
difference $q_{i}=m_{i+1}-m_{i}$, restricting to $m_{i+1}>m_{i}$ and $\delta
y_{i}=\delta r_{i}<Y=d_{0}$. The parameter $d_{0}$ is chosen sufficiently
small to ensure that earthquake pairs are close enough in space, guaranteeing
similarity in the distance to the stations necessary for their identification.
Under this condition, since the earthquake $m_{i}$ has been detected, it
implies that $m_{i}>M_{c}-2\sigma$, and therefore
$\Phi\left(m_{j}+M-M_{c}|m_{i}\right)\simeq 1$, for $M>2\sigma$. Accordingly,
from Eq.(4), for $M>2\sigma$, $P(q_{i}>M|\Delta y_{i}<Y)=e^{-\beta M}K$, where
$K$ is a constant given by $K=\int_{m_{L}}^{\infty}dm_{i}e^{-2\beta
m_{i}}\Phi\left(m_{i}-M_{c}\right)$. After derivation, we therefore obtain,
for $M>2\sigma$

$P(q_{i}=M|\Delta y_{i}<Y)=\beta Ke^{-\beta M},$ (6)

showing that, even in the presence of detection problems, the distribution of
magnitude differences $P(q_{i})$ follows an exponential GR decay for
$q_{i}>2\sigma$, independently of the previous seismic history. This behavior
is observed in instrumental catalogs (Fig.1b), where we find that $P(q_{i})$
follows a pure exponential decay when $q_{i}$ is larger than a completeness
threshold $\delta M_{c}=0.4$, for both SC and NC, consistent with Eq.(6), with
the same value of the exponent ($\beta=2.23\pm 0.05$ in SC, $\beta=2.20\pm
0.05$ in NC) extrapolated from the magnitude distribution $P(m_{i})$ (Eq.(3))
when $m_{i}>M_{c}$ (Fig.1b).

Figure 2: The quantity $E\left[\Delta P_{+}(Q,Y)\right]$ for the SC, where
$\Delta P_{+}(Q,Y)$ is defined in Eq.(7), with $\delta y_{i}=\delta t_{i}$ and
$Y=1$ h. Different curves correspond to different thresholds $\delta M_{th}$.
Error bars represent $2\Sigma(Q,Y)$ and are not shown for the curve with
$\delta M_{th}=1$ to enhance readability.

It is crucial to note that while detection issues introduce spurious
correlations between $m_{i}$ and $m_{i+1}$, positive magnitude differences
remain unaffected by such detection problems. Consequently, when
$q_{i}>2\sigma$, the null hypothesis predicts no correlations between $q_{i}$
and $q_{i+1}$. To delve deeper into this aspect, we select the subset of the
catalog $\mathcal{D}$, denoted as $\mathcal{D}^{\prime}$, which contains only
those earthquakes satisfying both constraints $\delta r_{i}<d_{0}$ and $\delta
m_{i}>0$, and therefore the $i$-th element of the catalog
$\mathcal{D}^{\prime}$ will be the vector ${q_{i},t_{i},x_{i},y_{i}}$. From
the catalog $\mathcal{D}^{\prime}$, we next obtain the catalog $\mathcal{DD}$,
whose elements are the vectors $\delta q_{i},t_{i},\vec{r}_{i}$ with $\delta
q_{i}=q_{i+1}-q_{i}$. It is important to note that, in general, $\delta
q_{i}\neq m_{i+2}-m_{i}$, and this equality only holds when
$m_{i+2}>m_{i+1}>m_{i}$ and $dr_{i}$ and $dr_{i+1}$ are both smaller than
$d_{0}$. We also generate several reshuffled catalogs
$\mathcal{DD}_{\text{ran}}$, where the $i$-th element is the vector $\delta
q^{*}_{i},t_{i},\vec{r}_{i}$ with $\delta q^{*}_{i}=q_{k}-q_{i}$, and $k$ is
the index of a random earthquake in the catalog $\mathcal{D}^{\prime}$. We can
define, $P(\delta q_{i}>Q|\delta y_{i}<Y)$, $\overline{P(Q,Y)}$, and
$\Sigma(Q,Y)$ by simply replacing $m_{i}$ with $q_{i}$ and $M$ with $Q$ in all
the definitions from Eq.(1) to Eq.(4). Accordingly, from Eq.(5), under the
null hypothesis Eq.(3), we obtain

$\displaystyle\Delta P_{+}(Q,Y)=P(\delta q_{i}>Q|\Delta y_{i}<Y)-P(\delta
q^{*}_{i}>Q|\Delta y_{i}<Y)$ (7)
$\displaystyle=\int_{0}^{\infty}dq_{i}\int_{q_{i}+Q}^{\infty}dq_{j}P(q_{i})P(q_{j})\Phi\left(q_{i}\right)\left[\Phi_{q}\left(q_{j}|q_{i}\right)-\Phi_{q}\left(q_{j}\right)\right],$

where $\Phi_{q}\left(q_{j}|q_{i}\right)$ and $\Phi_{q}(q_{j})$ represent the
conditioned and unconditioned probabilities, respectively, of detecting a
magnitude difference $q_{j}$. Correlations between $q_{i}$ and $q_{i+1}$ are
induced by the term $\Phi_{q}\left(q_{j}|q_{i}\right)$. Nevertheless, in
contrast to $\Phi\left(m_{j}-M_{c}|m_{i}\right)$ where the information that
event $m_{i}$ has been detected strongly constrains the detectability of
$m_{i+1}>m_{i}$, in this case, the information that the previous magnitude
difference $q_{i}$ has been detected only weakly affects the probability
$\Phi_{q}\left(q_{j}|q_{i}\right)$ to detect the next magnitude difference
$q_{i+1}$. In fact, we expect that
$\Phi_{q}\left(q_{j}|q_{i}\right)\simeq\Phi_{q}\left(q_{j}\right)$, and
$\Delta P_{+}(Q,Y)$ is expected to be very small under the null-hypothesis
scenario. In particular, by considering only magnitude differences
$q_{i}\geq\delta M_{th}$, as $\delta M_{th}$ increases, the detection
probability $\Phi_{q}$ also increases and approaches $1$ for $\delta
M_{th}\geq 2\sigma$. Consequently, we expect that $\Delta P_{+}(Q,Y)$ presents
small values when $\delta M_{th}=0$, and decreases as $\delta M_{th}$
increases, vanishing when magnitude differences are detected with probability
$1$ for $\delta M_{th}\geq 2\sigma$.

In the following we investigate the behavior of the mean value $E\left[\Delta
P_{+}(Q,Y)\right]$ for the SC catalog. The main result is that $E\left[\Delta
P_{+}(Q,Y)\right]$ exhibits a markedly different pattern (Fig.2) compared to
the one predicted by the null hypothesis (Eq.(7)). While data for $\delta
M_{th}=0$ could be consitent with the null hypothesis Eq.(7), as evidenced by
$E\left[\Delta P_{+}(Q,Y)\right]$ being only slightly larger than
$2\Sigma(Q,Y)$, the remarkable finding is that for larger values of $\delta
M_{th}$, $E\left[\Delta P_{+}(Q,Y)\right]$ not only fails to decrease but
rather increases. Particularly notable is that for $\delta M_{th}\geq 0.5$,
significantly surpassing $\delta M_{c}\simeq 0.4$, when the variables $q_{i}$
remain unaffected by detection issues (Fig.1b), the correlation between
$q_{i}$ and $q_{i+1}$ becomes even more pronounced. The significant growth of
$E\left[\Delta P_{+}(Q,Y)\right]$ as a function of $\delta M_{th}$ is notable
until $\delta M_{th}=0.75$, after which no substantial increase is observed up
to $\delta M_{th}=1$. Results plotted in Fig.2 unequivocally demonstrate the
presence of non-trivial correlations among magnitude differences.
Consequently, the result derived from the null hypothesis of i.i.d.
magnitudes, as expressed in Eq.(7), can be confidently dismissed.

We note that the results presented in Fig.2 were derived by considering only
earthquakes with epicentral distances smaller than $d_{0}=20$ km as elements
of the catalog ${\cal D}^{\prime}$. Similar results were obtained for other
choices of $d_{0}\leq 50$ km. Additionally, comparable outcomes were observed
when analyzing the NC catalog (Fig.3a). This suggests that correlations
between magnitude differences appear to be a consistent characteristic of
seismic catalogs, with patterns showing weak dependency on the geographic
region.

In order to gain further insights into the mechanisms responsible for
correlations between $q_{i+1}$ and $q_{i}$, we examine the influence of $Y$ on
$E\left[\Delta P(Q,Y)\right]$ in Fig.3. Specifically, in Fig.3d, we consider
$\delta y_{i}=\delta t_{i}$ and investigate decreasing values of $Y$, ranging
from $Y=\infty$ to $Y=0.1$ hours. The results show that as $Y$ decreases,
imposing that earthquake pairs $(q_{i},q_{i+1})$ are closer in time, their
correlation strengthens. In Fig.3c, we conduct a similar analysis but with
$\delta y_{i}=\delta r_{i}$, with $Y$ decreasing from $Y=\infty$ to $1$ km.
Here, we observe that magnitude correlations become more pronounced as we
transition from $Y=\infty$ to $Y=10$ km, but then stabilize, remaining roughly
constant regardless of the specific value of $Y$. Consequently, magnitude
correlations are stronger for earthquakes occurring closer in time, whereas
they are only weakly affected by the spatial distance between earthquakes.

Further insights are provided by examining the derivative $p^{\prime}(\delta
Q,Y)=-\left.\frac{\partial E\left[\Delta P(Q,Y)\right]}{\partial
Q}\right|_{\delta Q}$, which measures the difference in probability between
the instrumental and reshuffled catalogs to observe $q_{i+1}=q_{i}+\delta Q$.
The results (Fig.3b) reveal that it is more probable for a magnitude
difference to be followed by a larger one. Given that the inverse of the
average value of $q_{i}$ coincides with $\beta$ [40], the origin of magnitude
correlations could be linked to a decreasing $\beta$ value over time during
aftershock sequences. This observation aligns with the pattern identified by
Gulia & Wiemer [27] from their analysis of $58$ main-aftershock sequences. On
the contrary, the scarcity of earthquakes preceding significant shocks poses a
challenge in deriving conclusive insights from our analysis regarding the
behavior of $\beta$ for foreshocks, and in comparing them with the predictive
framework proposed in references [48, 49].

Figure 3: (a) The quantity $E\left[\Delta P_{+}(Q,Y)\right]$ for the SC and
NC, with $\delta y_{i}=\delta t_{i}$ and $Y=1$h and $\delta M_{th}=0.75$. (b)
The quantity $p^{\prime}(\delta Q,Y)$ for the SC, with $\delta y_{i}=\delta
t_{i}$ and $Y=1$h and different curves corresponding to different thresholds
$\delta M_{th}$. (c) The quantity $E\left[\Delta P_{+}(Q,Y)\right]$ for the
SC, with $\delta M_{th}=0.75$, $\delta y_{i}=\delta r_{i}$ and differnt values
of $Y$. (c) The quantity $E\left[\Delta P_{+}(Q,Y)\right]$ for the SC, with
$\delta M_{th}=0.75$, $\delta y_{i}=\delta t_{i}$ and different values of $Y$.
In panels (a,c,d) error bars represent $2\Sigma(Q,Y)$.

In summary, our findings unequivocally demonstrate the presence of non-trivial
correlations between subsequent magnitude differences, suggesting potential
modifications to be integrated into the ETAS model. One possible approach, as
proposed by [29, 30], is to consider two distinct $\beta$ values for
aftershocks, $\beta\pm\delta\beta$, depending on the magnitude of the
triggering mainshock. Alternatively, a time-varying $\beta$ value during the
aftershock sequence could also be implemented [50]. This novel research
direction holds promise for extracting valuable insights from the intricate
patterns of seismicity, ultimately enhancing our ability to forecast seismic
events and advancing our understanding of the mechanisms governing earthquake
triggering.

## References

  * Schorlemmer and Woessner [2008] D. Schorlemmer and J. Woessner, Bulletin of the Seismological Society of America 98, 2103 (2008).
  * Mignan  _et al._ [2011] A. Mignan, M. J. Werner, S. Wiemer, C.-C. Chen, and Y.-M. Wu, Bulletin of the Seismological Society of America 101, 1371 (2011).
  * Mignan and Woessner [2012] A. Mignan and J. Woessner, Community Online Resource for Statistical Seismicity Analysis 10.5078/corssa-00180805 (2012).
  * Kagan [2004] Y. Y. Kagan, Bulletin of the Seismological Society of America 94, 1207 (2004).
  * Helmstetter  _et al._ [2006] A. Helmstetter, Y. Y. Kagan, and D. D. Jackson, Bulletin of the Seismological Society of America 96, 90 (2006).
  * Peng  _et al._ [2007] Z. Peng, J. E. Vidale, M. Ishii, and A. Helmstetter, Journal of Geophysical Research: Solid Earth 112, n/a (2007), b03306.
  * Lippiello  _et al._ [2016] E. Lippiello, A. Cirillo, G. Godano, E. Papadimitriou, and V. Karakostas, Geophysical Research Letters 43, 6252 (2016), 2016GL069748.
  * Hainzl [2016a] S. Hainzl, Journal of Geophysical Research: Solid Earth 121, 6499 (2016a), 2016JB013319.
  * Hainzl [2016b] S. Hainzl, Seismological Research Letters 87, 337 (2016b).
  * de Arcangelis  _et al._ [2018] L. de Arcangelis, C. Godano, and E. Lippiello, Journal of Geophysical Research: Solid Earth 123, 5661 (2018).
  * Petrillo  _et al._ [2020] G. Petrillo, F. Landes, E. Lippiello, and A. Rosso, Nature Communications 11, 3010 (2020).
  * Hainzl [2021] S. Hainzl, Bulletin of the Seismological Society of America 112, 494 (2021).
  * Lippiello  _et al._ [2019] E. Lippiello, C. Petrillo, C. Godano, A. Tramelli, E. Papadimitriou, and V. Karakostas, Nature Communications 10, 2953 (2019).
  * Lippiello  _et al._ [2007a] E. Lippiello, C. Godano, and L. de Arcangelis, Phys. Rev. Lett. 98, 098501 (2007a).
  * Lippiello  _et al._ [2007b] E. Lippiello, M. Bottiglieri, C. Godano, and L. de Arcangelis, Geophysical Research Letters 34, L23301 (2007b).
  * Lippiello  _et al._ [2008] E. Lippiello, L. de Arcangelis, and C. Godano, Phys. Rev. Lett. 100, 038501 (2008).
  * Shcherbakov  _et al._ [2005] R. Shcherbakov, G. Yakovlev, D. L. Turcotte, and J. B. Rundle, Phys. Rev. Lett. 95, 218501 (2005).
  * Varotsos  _et al._ [2005] P. A. Varotsos, N. V. Sarlis, H. K. Tanaka, and E. S. Skordas, Phys. Rev. E 72, 041103 (2005).
  * Sarlis  _et al._ [2010a] N. V. Sarlis, E. S. Skordas, and P. A. Varotsos, Phys. Rev. E 82, 021110 (2010a).
  * Sarlis  _et al._ [2010b] N. V. Sarlis, E. S. Skordas, and P. A. Varotsos, Phys. Rev. E 82, 021110 (2010b).
  * Sarlis [2011] N. V. Sarlis, Phys. Rev. E 84, 022101 (2011).
  * Lippiello  _et al._ [2012] E. Lippiello, C. Godano, and L. de Arcangelis, Geophysical Research Letters 39, L05309 (2012).
  * Davidsen and Green [2011] J. Davidsen and A. Green, Phys. Rev. Lett. 106, 108502 (2011).
  * Nichols and Schoenberg [2014] K. Nichols and F. P. Schoenberg, Environmetrics 25, 143 (2014).
  * de Arcangelis  _et al._ [2016] L. de Arcangelis, C. Godano, J. R. Grasso, and E. Lippiello, Physics Reports 628, 1 (2016).
  * Spassiani and Sebastiani [2016a] I. Spassiani and G. Sebastiani, Phys. Rev. E 93, 042134 (2016a).
  * Gulia  _et al._ [2018] L. Gulia, A. P. Rinaldi, T. Tormann, G. Vannucci, B. Enescu, and S. Wiemer, Geophysical Research Letters 45, 13,277 (2018).
  * Ogata  _et al._ [2018] Y. Ogata, K. Katsura, H. Tsuruoka, and N. Hirata, Seismological Research Letters 89, 1298 (2018).
  * Nandan  _et al._ [2019] S. Nandan, G. Ouillon, and D. Sornette, Journal of Geophysical Research: Solid Earth 124, 2762 (2019).
  * Nandan  _et al._ [2022] S. Nandan, G. Ouillon, and D. Sornette, Journal of Geophysical Research: Solid Earth 127, e2022JB024380 (2022), e2022JB024380 2022JB024380.
  * Petrillo and Zhuang [2022] G. Petrillo and J. Zhuang, Sci Rep 12, 20683 (2022).
  * Petrillo and Zhuang [2023] G. Petrillo and J. Zhuang, Phys. Rev. Lett. 131, 154101 (2023).
  * Xiong  _et al._ [2023] Q. Xiong, M. R. Brudzinski, D. Gossett, _et al._ , Nature Communications 14, 2056 (2023).
  * Taroni [2024] M. Taroni, Geophysical Journal International 236, 1596 (2024).
  * Ogata [1985] Y. Ogata, Research Memo. Technical report Inst. Statist. Math., Tokyo. 288 (1985).
  * Ogata [1988a] Y. Ogata, J. Amer. Statist. Assoc. 83, 9 – 27 (1988a).
  * Ogata [1988b] Y. Ogata, Ann. Inst. Math.Statist. 50, 379–402 (1988b).
  * Ogata [1989] Y. Ogata, Numerische Mathematik 55, 137 (1989).
  * Spassiani and Sebastiani [2016b] I. Spassiani and G. Sebastiani, Journal of Geophysical Research: Solid Earth 121, 903 (2016b), 2015JB012398.
  * van der Elst [2021] N. J. van der Elst, Journal of Geophysical Research: Solid Earth 126, e2020JB021027 (2021), e2020JB021027 2020JB021027.
  * Lippiello and Petrillo [2024] E. Lippiello and G. Petrillo, Journal of Geophysical Research: Solid Earth 129, e2023JB027849 (2024), e2023JB027849 2023JB027849.
  * van der Elst and Page [2023] N. J. van der Elst and M. T. Page, Journal of Geophysical Research: Solid Earth 128, e2023JB027089 (2023), e2023JB027089 2023JB027089.
  * van der Elst [2023] N. J. van der Elst, AGU assembly (2023).
  * Hauksson  _et al._ [2012] E. Hauksson, W. Yang, and P. M. Shearer, Bulletin of the Seismological Society of America 102, 2239 (2012).
  * Waldhauser and Schaff [2008] F. Waldhauser and D. P. Schaff, Journal of Geophysical Research: Solid Earth 113 (2008).
  * Ogata and Katsura [1993] Y. Ogata and K. Katsura, Geophysical Journal International 113, 727 (1993).
  * Woessner and Wiemer [2005] J. Woessner and S. Wiemer, Bulletin of the Seismological Society of America 95, 684 (2005).
  * Gulia and Wiemer [2019] L. Gulia and S. Wiemer, Nature 574, 193 (2019).
  * Gulia  _et al._ [2020] L. Gulia, S. Wiemer, and G. Vannucci, Seismological Research Letters 91, 2828 (2020).
  * Gulia  _et al._ [2016] L. Gulia, T. Tormann, S. Wiemer, M. Herrmann, and S. Seif, Geophysical Research Letters 43, 1100 (2016).

###### Acknowledgements.

E.L. and L.de A. acknowledge support from the MIUR PRIN 2022 PNRR P202247YKL,
C.G. acknowledges support from the MIUR PRIN 2022 PNRR P20222B5P9