Relationship of Multifractal and Entropic Properties of Global Seismic Noise with Major Earthquakes, 1997–2025

1. Introduction

The introduction should briefly place the study in a broad context and highlight why it is important. It should define the purpose of the work and its significance. The current state of the research field should be carefully reviewed and key publications cited. Please highlight controversial and diverging hypotheses when necessary. Finally, briefly mention the main aim of the work and highlight the principal conclusions. As far as possible, please keep the introduction comprehensible to scientists outside your particular field of research. References should be numbered in order of appearance and indicated by a numeral or numerals in square brackets—e.g., [1] or [2,3], or [4,5,6]. See the end of the document for further details on references.

Seismic noise is a reflection of the inner life of the planet and is an important source of information that makes it possible to study processes in the lithosphere, including those preceding strong earthquakes [1,2,3]. The main source of energy for the seismic background of the Earth is the movement of cyclones in the atmosphere and the impact of ocean waves on the shelf and coast [4,5,6,7,8,9,10]. Since the earth's crust is a medium for the propagation of seismic waves from the ocean and the atmosphere, it is logical to assume that the processes inside the earth's crust are reflected in changes in the statistical properties of seismic noise, and the study of these properties makes it possible to determine the structural features of the earth's crust [11].

The global seismic noise, recorded on a network of 229 broadband seismic stations located around the globe for 25 years, from the beginning of 1997 to the end of 2025, has been investigated. To study the noise properties a set of statistics estimated daily have been used: the multifractal singularity spectrum support width, the minimum entropy of the squared wavelet coefficients, and the wavelet-based Donoho-Johnstone (DJ) index. It is shown that the time points of local extrema of the average values of the analyzed noise properties (minima for singularity spectrum support width and DJ-index and maxima for entropy) tend to occur before strong earthquakes.

The article continues the studies carried out in [12,13] to investigate the correlation and coherence properties of low-frequency seismic noise on a global scale covering the entire planet.

2. Data

The article continues the studies carried out in [12,13] to investigate the correlation and coherence properties of low-frequency seismic noise on a global scale covering the entire planet.

The data used are the vertical components of continuous seismic noise records with 1 second sampling intervals, which were downloaded from the Incorporated Research Institutions for Seismology (IRIS) website [14] from 229 broadband seismic stations of 3 networks: http://www.iris.edu/mda/_GSN , http://www.iris.edu/mda/G , http://www.iris.edu/mda/GE

Seismic noise records with a sampling rate of 1 Hz (LHZ records) were considered for 29 years of registration (from January 1, 1997 to December 31, 2025). These data were converted to 1-minute time series by calculating averages for successive 60-second time intervals.

Figure 1. Global seismic noise observation network. Blue dots represent the locations of 229 broadband seismic stations, and numbered red circles represent the locations of 50 reference points.

Figure 1. Global seismic noise observation network. Blue dots represent the locations of 229 broadband seismic stations, and numbered red circles represent the locations of 50 reference points.

3. Properties of Seismic Noise

Seismic noise records with a sampling rate of 1 Hz (LHZ records) were considered for 29 years of registration (from January 1, 1997 to December 31, 2025). These data were converted to 1-minute time series by calculating averages for successive 60-second time intervals.

The minimum entropy of a time series is determined by the formula $En=-{\displaystyle \sum _k{q}_k\cdot \ln (}{q}_k)/\ln (N)$, where $q_k=w_k^{ 2}/{\displaystyle {\sum }_j{w}_j^{ 2}}$, $w_k$ are the wavelet coefficients of the signal decomposition, and $N$ is the total number of coefficients $w_k$. Seventeen orthogonal Daubechies wavelets were used: 10 ordinary minimum support bases with a number of vanishing moments from 1 to 10 and 7 so-called Daubechies symlets [15] with a number of vanishing moments from 4 to 10. In each time window, the wavelet for which the value $En$ is minimal is selected.

The Donoho-Johnston (DJ) wavelet-based index $\gamma$ is the ratio of "large" wavelet coefficients in absolute value to their total number. By definition $0\le \gamma \le 1$. The threshold separating "large" wavelet coefficients is $\sigma \sqrt{2\cdot \ln N}$, where $\sigma =med\left\{\left|w_k^{(1)}\right|\hspace{0.33em},\hspace{0.33em}k=1,\mathrm{...},N/2\right\}/0.6745$ is a robust estimate of the standard deviation of the normal distribution, $w_k^{(1)}$ are the wavelet coefficients at the first detail level of decomposition, and $N/2$ is the number of such coefficients [15,16].

The singularity spectrum support width $\Delta \alpha$ is considered as a measure of the diversity of the stochastic behavior of the signal $u(t)$. It is defined as $\Delta \alpha ={\alpha }_{\max }-{\alpha }_{\min }$, where ${\alpha }_{\max }$ and ${\alpha }_{\max }$ are the minimum and maximum values of the Hölder-Lipschitz exponent [17], which governs the behavior of the signal in the vicinity of the time instant $t$: $|u(t+\frac{\delta }{2})-u(t-\frac{\delta }{2})|\hspace{0.33em}\sim \hspace{0.17em}|\delta {|}^{\hspace{0.17em}\alpha },\hspace{1em}\delta \to 0$. For a monofractal signal, the exponent $\alpha$ is the same for all time instants $t$. If this exponent differs, then the signal is multifractal [18].

After switching to a 1-minute sampling time step, the seismic records from each station were divided into 1-day-long time fragments (1440 samples), and the parameters $(\gamma ,\Delta \alpha ,En)$ of daily seismic noise signals were calculated for each fragment. The DFA method [19] was used to calculate the values of the singularity spectrum support width. The methods for estimating the parameters $\gamma$, $\Delta \alpha$ and $En$ for seismic noise records in a sliding time window are described in detail in [20,21]. Thus, time series of values with a 1-day time step were obtained for each of the seismic stations.

The entropy value $En$ used in this paper, defined through the coefficients of the orthogonal wavelet decomposition, has features in common with multiscale entropy used, for example, in biology [22,23]. In [24,25], entropy is used within the framework of the natural time approach for seismic data analysis. In [26], non-extensive Tsallis entropy is used for processing seismic noise data.

The use of multifractal analysis to study the behavior of various complex systems has a long history. Particular attention is paid to the effect of parameter $\Delta \alpha$ reduction (loss of multifractality) preceding changes in system properties. In medicine, a decrease in the $\Delta \alpha$ value of various parameters accompanies age-related changes [27,28,29]. In [30,31], multifractal analysis is used to analyze geoelectric signals and wind speed. Within the framework of the natural time approach, multifractal analysis is used to study both seismicity and other time sequences [32].

Let us denote by ${\tilde{\gamma }}_j(t)$, $\Delta {\tilde{\alpha }}_j(t)$, $\tilde{E}{n}_j(t)$, $j=1,\mathrm{...50}$, the daily time series of parameters at 50 reference points of the seismic noise observation network, which are calculated as the medians of the values at the 5 nearest operational seismic stations.

4. First Principal Component and Weighted Mean of Multiple Time Series

In the future, we will need to calculate the first principal components and weighted averaging of multivariate time series. The values of the parameters $(\gamma ,\Delta \alpha ,En)$ introduced above are calculated in successive time windows of a certain length, resulting in a 3-dimensional time series, the properties of which are further studied jointly. The used properties of seismic noise reflect the change in its structure, in particular, we are interested in the phenomenon of noise structure simplification as a sign that precedes strong earthquakes. Attempts to determine the "best" property of noise led to the idea of using the principal component approach [33,34,35] to aggregate time series $(\gamma ,\Delta \alpha ,En)$ into one scalar time series. Since the purpose of the analysis is to study the variability of noise properties both in time and space, the principal component method was applied in a sliding time window.

Let’s consider a multiple time series $p(t)={(p_1(t),\mathrm{...},p_m(t))}^T$, $t=0,1,\mathrm{...}$ of the dimensionality $m$. It is necessary to estimate two types of its first principal component in the moving time window of the length $L$ samples. For this purpose let’s consider samples with time indices $t$ under the condition $s-L+1\le t\le s$ where $s$ is the right-most end of moving time window. Let’s perform normalization of multiple time series components within current time window:

$$q_{\hspace{0.33em}k}^{(s)}(t)=(p_k(t)-{\overline{p}}_{\hspace{0.33em}k}^{(s)}),\hspace{1em}{\overline{p}}_{\hspace{0.33em}k}^{(s)}={\displaystyle \sum _{t=s-L+1}^s{p}_k(t)/L},\hspace{1em}k=1,\mathrm{...},m$$

(1)

$$r_{\hspace{0.33em}k}^{(s)}(t)=q_{\hspace{0.33em}k}^{(s)}(t)/{\sigma }_{\hspace{0.33em}k}^{(s)},\hspace{1em}{\left({\sigma }_{\hspace{0.33em}k}^{(s)}\right)}^2={\displaystyle \sum _{t=s-L+1}^s(q_{\hspace{0.33em}k}^{(s)}(t)}{)}^2/(L-1),\hspace{1em}k=1,\mathrm{...},m$$

(2)

Let’s consider two correlation matrices $Q\hspace{0.17em}(s)$ and $R (s)$ of the size $m\times m$ which are calculated by the formulae:

$$Q(s)=\left(Q_{\hspace{0.33em}kj}^{(s)}\right),\hspace{0.33em}\hspace{0.33em}{Q}_{\hspace{0.33em}kj}^{(s)}={\displaystyle \sum _{t=s-L+1}^s{q}_{\hspace{0.33em}k}^{(s)}(t)q_{\hspace{0.33em}j}^{(s)}(t)/L,\hspace{1em}k,j=1,\mathrm{...},m}$$

(3)

$$R(s)=\left(R_{\hspace{0.33em}kj}^{(s)}\right),\hspace{0.33em}\hspace{0.33em}{R}_{\hspace{0.33em}kj}^{(s)}={\displaystyle \sum _{t=s-L+1}^s{r}_{\hspace{0.33em}k}^{(s)}(t)r_{\hspace{0.33em}j}^{(s)}(t)/L,\hspace{1em}k,j=1,\mathrm{...},m}$$

(4)

Let ${\varphi }^{(s)}={({\varphi }_{\hspace{0.33em}1}^{(s)},\mathrm{...},{\varphi }_{\hspace{0.33em}m}^{(s)})}^T$ and ${\phi }^{(s)}={({\phi }_{\hspace{0.33em}1}^{(s)},\mathrm{...},{\phi }_{\hspace{0.33em}m}^{(s)})}^T$ be eigenvectors of the matrices $Q(s)$ and $R(s)$ corresponding to their maximum eigenvalues. Let’s define first principal component $P^{(s)}(t)$ and weighted mean $W^{(s)}(t)$ of multiple time series $p(t)$ within current time window by formulae:

$$P^{(s)}(t)={\displaystyle \sum _{k=1}^m{\varphi }_{\hspace{0.33em}k}^{(s)}\cdot r_{\hspace{0.33em}k}^{(s)}(t)},\hspace{1em}{W}^{(s)}(t)={\displaystyle \sum _{k=1}^m{({\phi }_{\hspace{0.33em}k}^{(s)})}^2\cdot q_{\hspace{0.33em}k}^{(s)}(t)}$$

(5)

and define scalar time series $P(t)$ and $W(t)$ of the adaptive first principal component and weighted mean in the moving time window by the formulae:

$$P(t)=\left\{\begin{array}{l}{P}^{(L-1)}(t), 0\le t\le (L-1)\\ P^{(t)}(t), t\ge L\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\end{array}\right.,\hspace{1em}W(t)=\left\{\begin{array}{l}{W}^{(L-1)}(t), 0\le t\le (L-1)\\ W^{(t)}(t), t\ge L\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\end{array}\right.$$

(6)

Formulae (1-6) are applied independently within each time window. According to these formulae within first time window time series $P(t)$ and $W(t)$ consist of $L$ values corresponding to (3-5). In all subsequent windows $P(t)$ and $W(t)$ correspond the only sample in the right-most time index. Thus, outside the first time window $P(t)$ and $W(t)$ are dependent on the past values of $p (t)$.

The standard first component $P(t)$ is used when it is necessary to aggregate heterogeneous components of a multidimensional time series, for example, a daily time series of seismic noise properties $(\gamma ,\Delta \alpha ,En)$ at a single reference point. In this case the time series dimension is $m=3$. Weighted averaging $W(t)$ of time series components is used when the time series components are homogeneous, for example, to average the time series of principal components at all reference points. In this case the time series dimension is $m=50$. In the following, we will agree to call these two variants of using the principal component method P-averaging and W-averaging.

Hereinafter, a window with a length of 365 days will be used. The choice of this length is quite natural because of presence of annual periodicity in almost all background processes in Earth’s crust. At each reference point with a number $j$ we calculate the first principal component of three time series ${\tilde{\gamma }}_j(t)$, $\Delta {\tilde{\alpha }}_j(t)$, $\tilde{E}{n}_j(t)$ using the P-averaging method (1-6) in a sliding time window of 365 days in length and denote it as $P_{\gamma ,\Delta \alpha ,En}^{\hspace{0.33em}\hspace{0.33em}(j)}(t)$. It should be noticed that values of seismic noise parameters and their first principal component are dimensionless. Figure 2 shows graphs of time series of seismic noise properties and their first principal components, calculated in a sliding adaptation time window of 365 days for 6 reference points.

Let's estimate the temporal variability of strong correlations between the first principal components of seismic noise properties at different reference points. To do this, we calculate the average absolute correlations of the principal components $P_{\gamma ,\Delta \alpha ,En}^{\hspace{0.33em}\hspace{0.33em}(j)}(t)$ between all reference points over a 365-day sliding time window. Figure 3(a) shows a graph of these average correlations, which shows a trend of increasing spatial correlations since 2003, superimposed by oscillations with a period of approximately 2.5 years. We set a threshold of 0.8 for correlations and calculate the change in the number of reference point pairs for which the correlation exceeds this threshold. This dependence is shown in Figure 3(b), which shows that the number of strong spatial correlations increased essentially after 2013.

The time point 2013, corresponding to the right end of the time window, was found formally by the method of variance analysis [34] by comparing the mean values in two time intervals between the trial time point separating these intervals for the fragment 1998-2025. The final fragment of the positions of the right end of the time window 2025-2026 stands out especially due to a remarkable anomaly in the behavior of the number of pairs of reference points with strong correlations: in the second half of 2025 the number of strong correlations tripled in six months. Figure 3(c) shows the dependence of the maximum distances between reference points for which a strong correlation greater than 0.8 emerged on the position of the right end of a 365-day sliding time window. This dependence also reveals a step in the distances of strong spatial correlations after 2012.

In Figure 4, strong correlations between the data points are visualized as two connectivity diagrams for time windows spanning two 2-years’ time intervals. Figure 4(b) corresponds to the time window with the largest number of strong correlations.

5. Probability Densities of Extreme Values

In what follows, we will distinguish two cases of considering extreme values of the seismic noise statistics used. We will agree to denote the first case by $({\gamma }_{\min },\Delta {\alpha }_{\min },E{n}_{\max })$, in which we analyze the joint spatial distribution or local extrema in the time realizations of the minimum values of the Donoho-Johnston index $\gamma$ and the singularity spectrum support width $\Delta \alpha$ and the maximum values of entropy $En$. We will denote the second case by $({\gamma }_{\max },\Delta {\alpha }_{\max },E{n}_{\min })$, for which we consider the joint spatial distribution or local extrema in the time realizations of the maximum values of the Donoho-Johnston index $\gamma$ and the singularity spectrum support width $\Delta \alpha$ and the minimum values of entropy $En$. These two cases are antagonistic. Case 1 corresponds to a simplification of the seismic noise structure and its properties approaching white noise. The opposite case 2 corresponds to enrichment of the noise structure and, as a rule, corresponds to the appearance of chaotic high-amplitude spikes in noise realizations.

Let us study the variability of the spatial distribution of extreme values of seismic noise properties. For this purpose, we consider a regular grid of 100 nodes by longitudes and 50 nodes by latitudes, covering all area of the Earth. Let $U$ be any value of $\gamma$, $\Delta \alpha$ or $En$. For each grid node $(i,j)$ and for each day with number $t$, we find the 5 nearest working seismic stations, which provides 5 values of $U$. Let us denote by $U_{ij}^{(t)}$ the median value of these 5 properties at a node $(i,j)$ on a day with number $t$. For each daily set of $U_{ij}^{(t)}$ values with a discrete time index $t$, we find the coordinates of the nodes ${\zeta }_{\hspace{0.17em}mn}^{(t)}=(Lo{n}_m^{(t)},La{t}_n^{(t)})$ at which extreme values of $U$ is reached relative to all other nodes of the regular grid for 2 cases which we introduced above, $({\gamma }_{\min },\Delta {\alpha }_{\min },E{n}_{\max })$ and $({\gamma }_{\max },\Delta {\alpha }_{\max },E{n}_{\min })$. We estimate the two-dimensional probability distribution function of the vectors ${\zeta }_{\hspace{0.17em}mn}^{(t)}$ within the time interval $t\in [t_0,t_1]$ for each node ${\zeta }_{ij}$ of the regular grid using the Parzen-Rosenblatt estimate with the Gaussian kernel function [35]:

$$p({\zeta }_{ij}|t_{ 0},t_{ 1})={\displaystyle \frac{1}{2\pi D^2(h)M_{t_0,t_1}}}{\displaystyle \sum _{t=t_{ 0}}^{t_{ 1}}{\displaystyle \sum _{mn}\exp \left(-\frac{{\rho }^2({\zeta }_{ij},{\zeta }_{\hspace{0.17em}mn}^{(t)})}{2D^2(h)}\right)}}$$

(7)

Here $h$ is the kernel averaging radius in spherical degrees, $t_{ 0},t_{ 1}$ are integer indices that number the daily elementary maps, and $M_{t_0,t_1}=(t_{ 1}-t_{ 0}+1)$ is the number of daily maps in the time interval under consideration, ${\rho }^2(\xi ,\zeta )$ is a squared spherical distance between points $\xi$ and $\zeta$ at the Earth surface, $D^2(h)$ is the square of the spherical distance between points on the Earth's surface, corresponding to the value $h$. A smoothing bandwidth of $h={15}^{\circ }$ was used.

Let us calculate the 2-dimensional distribution densities (7) of extreme values in successive time windows of 10 days in length ($M_{t_0,t_1}=10$). Since we have three such distributions, we will calculate their weighted average by taking as weights the squares of the components of the eigenvector of the correlation matrix of probability densities corresponding to the maximum eigenvalue of the matrix [33,34,35]. By construction, the sum of the squares of the components is equal to 1.

Figure 5 presents averaged distribution maps of extreme seismic noise properties for two cases $({\gamma }_{\min },\Delta {\alpha }_{\min },E{n}_{\max })$ and $({\gamma }_{\max },\Delta {\alpha }_{\max },E{n}_{\min })$ for three consecutive time intervals. The first interval covers the time interval from the beginning of observations to the time of the Tohoku mega-earthquake in Japan on March 11, 2011. The next time interval extends from the Tohoku event to the end of 2024. The final time interval begins in early 2024 and extends until the end of 2025. The final time interval, the shortest, was chosen because it contains a sharp increase in the number of pairs of control points with correlations exceeding the 0.8 threshold (Figure 3(b)).

As can be seen from Figures 5 (a1, a2, a3), the concentration regions of the most "quiet" seismic noise behavior (variant $({\gamma }_{\min },\Delta {\alpha }_{\min },E{n}_{\max })$) are stable throughout the entire observation interval, with the exception with the exception of the middle of the Atlantic Ocean in the vicinity of the island of Saint Edena and the south of the Indian Ocean in the vicinity of the Kerguelen archipelago, where a gradual increase in probability density is observed. As for the variant $({\gamma }_{\max },\Delta {\alpha }_{\max },E{n}_{\min })$ of the most "violent" seismic noise behavior, the concentration region of this behavior remains stable until 2024 and occupies northeastern Eurasia, centered on the Putarana Plateau, which formed after the largest supervolcanic eruption 250 million years ago. However, after 2024, this region shifted to the Middle East.

The graphs in Figure 6 provide a more detailed representation of the spatiotemporal features of the changes in the distribution densities of the extreme values of seismic noise properties. The graphs in Figure 6(a1) and Figure 6(b1) represent the changes in the latitude of those nodes of a regular 100×50 grid for which the maximum distribution of extreme properties is realized for the cases $({\gamma }_{\min },\Delta {\alpha }_{\min },E{n}_{\max })$ and $({\gamma }_{\max },\Delta {\alpha }_{\max },E{n}_{\min })$.

Figure 6(a1) shows a strong annual periodicity, with the point of maximum probability density $({\gamma }_{\min },\Delta {\alpha }_{\min },E{n}_{\max })$ making a sudden transition to the Arctic region in December–January. Furthermore, Figure 6(a1) shows a transition from the point of maximum probability density from 15 degrees north latitude to 30 degrees south latitude, occurring around 2011, which is noticeable when comparing Figure 5(a1) and Figure 5(a2).

As for Figure 6(b1), it shows a gradual evolution of the region of maximum probability densities for the case $({\gamma }_{\max },\Delta {\alpha }_{\max },E{n}_{\min })$ from Northeast Eurasia to the Middle East, which began around 2019.

Figure 6(a2) and Figure 6(b2) present graphs of probability histograms of the corresponding latitude values with maximum values of the distribution densities of extreme values of seismic noise properties.

6. Sequence of Major Earthquakes and Its Periodic Components

The following sections of the article will analyze the relationships between the properties of global seismic noise and the sequence of strong earthquakes. Figure 7 shows the time sequence of 433 seismic events with a magnitude of at least 7 for the time interval 1997–2025.

As Figure 7 shows, the regime of the strongest earthquakes is non-stationary. In particular, we will be interested in the periodic components of the intensity change in the earthquake sequence. Below, the method proposed in [37] is used to estimate the periodic components of the intensity of the sequence of events. In [38,39], this method was used to calculate the periodic component of the stepwise variations in the time series of the displacement of the earth’s surface measured by GPS and periodic components of earthquake sequence.

Let $t_i,\hspace{0.33em}i=1,\mathrm{...},N$ be the times of the sequence of events observed on the interval $(0,T]$. Consider the following intensity model containing a periodic component:

$$\lambda (t)=\mu \cdot (1+a\cos (\omega t+\phi ))$$

(8)

where frequency $\omega$, amplitude $a,\hspace{0.33em}0\le \hspace{0.17em}a\hspace{0.17em}\le \hspace{0.17em}1$, phase angle $\phi \in [0,2\pi ]$ and multiplier $\mu >0$ (describing the Poisson part of the intensity) are parameters of the model. Thus, the Poisson part of the intensity is modulated by a harmonic oscillation. Let us fix some value of frequency $\omega$. The logarithmic likelihood function [40] in this case for a series of observed events is equal to

$$\begin{array}{l}\ln L(\mu ,a,\phi |\omega )={\displaystyle \sum _i\ln (\lambda (t_i))}-{\displaystyle {\int }_0^T\lambda (s)ds}=\\ =N\ln (\mu )+{\displaystyle \sum _i\ln (1+a\cos (\omega t_i+\phi ))}-\mu T-{\displaystyle \frac{\mu a}{\omega }}[\sin (\omega T+\phi )-\sin (\phi )]\end{array}$$

(9)

Taking the maximum of expression (9) with respect to the parameter $\mu$, it is easy to find that

$$\ln (L(\widehat{\mu },a,\phi |\omega ))={\displaystyle \sum _i\ln (1+a\cos (\omega t_i}+\phi ))+N\cdot \ln (\widehat{\mu }(a,\phi |\omega ))-N$$

(10)

where $\widehat{\mu }(a,\phi |\omega )=N/(T+a\cdot ((\sin (\omega T+\phi )-\sin \phi )/\omega )$. It should be noted that the expression $\widehat{\mu }(a=0,\phi |\omega )\equiv {\mu }_0=N/T$ is an estimate of the intensity of the process under the condition that it is purely random. Thus, the increment of the logarithmic likelihood function due to the consideration of a richer intensity model with a harmonic component with a given frequency $\omega$ than for a purely random flow of events is equal to

$$\Delta \ln L(a,\phi |\omega )={\displaystyle \sum _i\ln (1+a\cos (\omega t_i}+\phi ))+N\cdot \ln (\widehat{\mu }(a,\phi |\omega )/{\mu }_0)$$

(11)

Let

$$R(\omega )=\underset{a,\phi }{\max }\Delta \ln L(a,\phi |\omega ),\hspace{0.33em}\hspace{0.33em}0\le \hspace{0.17em}a\hspace{0.17em}\le \hspace{0.17em}1,\hspace{0.33em}\phi \in [0,2\pi ]$$

(12)

As shown in the works [37,38,39], based on the application of Wilks’ theorem [34], if the hypothesis about the presence of a periodic component of the intensity of a point process is valid, then an asymptotic distribution takes place:

$$\mathrm{P}(R(\omega )<x)=1-e^{-x},\hspace{1em}N\to \infty$$

(13)

Expression (13) allows us to set thresholds for statistics that allow us to assert that only when they are exceeded does the sequence of time moments differ from the Poisson sequence with a given probability.

For a time sequence of seismic events with a magnitude of at least 7, we calculated the increments of the logarithmic likelihood function (12) in a sliding time window of 5 years with a shift of 0.01 years for 100 frequency values $\omega$ corresponding to the values of periods varying from 1 to 5 years with a uniform step in a logarithmic scale. The number of events within moving time window of the length 5 years is presented at Figure 8(a). The average number of events within the sliding window is 75 and it varies from 59 to 97. Therefore, we believe that the conditions for applicability of the asymptotic inequality (13) are applicable. It is interesting to note in Figure 8(a) that the number of events $M\ge 7$ in the 5-year time window started to increase after the Maule mega-earthquake in Chile on February 27, 2010, M=8.8, and reached its maximum values after the Tohoku mega-earthquake in Japan on March 11, 2011, M=9.1.

The resulting time–frequency dependence is shown in Figure 8(b). Figure 8(c) shows a plot of the mean values of the log-likelihood increments for all 5-year time windows. It shows a spectral peak at a period of 2.6 years. This plot was obtained by summing over all time windows, and the maximum increment value $\Delta \ln L$ is 1.55, which gives a probability of 0.787 for accepting the hypothesis of the presence of a periodic intensity component for a period of 2.6 years, according to formula (15). However, if we calculate the value $\Delta \ln L$ using the entire sample of seismic event times, we obtain the plot shown in Figure 8(d), in which the maximum value $\Delta \ln L$ at a period of 2.6 years is 5.55, which gives a probability of accepting the periodicity hypothesis of 0.996.

To determine the dominant period in the earthquake sequence, one could rely solely on the graph in Figure 8(d). However, the time-frequency diagram in Figure 8(b) is of independent interest, as it provides information on the evolution of the periodic intensity components. It shows that the period gradually increases until approximately 2013, after which the periodic component virtually disappears until 2018. After 2018, two branches of the periodic component's evolution emerge—low-frequency and high-frequency—which merge at the end of the analyzed time interval.

An analysis of the periodic components of the seismic process yielded an interesting result: it turned out that the dominant period of seismicity, 2.6 years, is close to the 2.5-year period of oscillations of the average value of pairwise absolute correlations between the properties of seismic noise at reference points – Figure 3(a).

7. Weighted Means of Seismic Noise Properties and Their Local Extrema

We apply the method for calculating weighted averages (“W-averaging” in formulae (5-6)) to the analyzed seismic noise statistics. The weights for averaging are the squared absolute values of the first eigenvectors of the 50×50 correlation matrices calculated for daily values of seismic noise properties ${\tilde{\gamma }}_j(t)$, $\Delta {\tilde{\alpha }}_j(t)$, $\tilde{E}{n}_j(t)$ at 50 reference points, corresponding to the maximum eigenvalues of these matrices. The weighted average values of these properties of seismic noise, calculated in a sliding time window of 365 days, are denoted by $W_{\gamma }(t)$, $W_{\Delta \alpha }(t)$ and $W_{En}(t)$. The results of such a “left-oriented” weighted “W-average” are presented in Figure 9 (a1, a2, a3).

Three time intervals can be distinguished in the behavior of the weighted average properties of seismic noise: 1997-2004, 2004-2016, and 2016-2025, when their linear trend changes. For the last two intervals, 2004-2016 and 2016-2025, the linear trend is virtually absent and coincides with the mean value. For the 2004-2016 time interval, a linear trend of decreasing weighted average values $W_{\gamma }(t)$, $W_{\Delta \alpha }(t)$ is observed, while for $W_{En}(t)$ an increasing trend is present. In Figure 9 (a1, a2, a3), these trends are shown by the red lines.

In the future, we will be interested in the positions of the time points of a certain number of the largest local maxima or the smallest local minima of daily seismic noise properties in comparison with the times of earthquakes with a magnitude of at least 7 (Figure 7). To eliminate the influence of low-frequency components of changes in statistical values on the determination of the times of local extrema, the time series of noise properties were subjected to the operation of removing low frequencies using Gaussian kernel smoothing. Let $z(t)$ be a time series with discrete time $t$. Gaussian kernel averaging of a time series $z(t)$ with radius (scale parameter) $a>0$ at time $t$, is calculated using the formula [41]:

$$\overline{z}(t|a)={\displaystyle \sum _{\xi }z(\xi )\cdot e^{-{\left(\frac{\xi -t}{a}\right)}^2}}/{\displaystyle \sum _{\xi }{e}^{-{\left(\frac{\xi -t}{a}\right)}^2}}$$

(14)

Let us denote by $W_{\hspace{0.33em}\gamma }^{(a)}(t)$, $W_{\hspace{0.33em}\Delta \alpha }^{(a)}(t)$ and $W_{\hspace{0.33em}En}^{(a)}(t)$ the result of removing local trends from signals $W_{\gamma }(t)$, $W_{\Delta \alpha }(t)$ and $W_{En}(t)$ by Gaussian window of radius $a$, equal to 2 days. We will be interested in the local extremum points of signals $W_{\hspace{0.33em}\gamma }^{(a)}(t)$, $W_{\hspace{0.33em}\Delta \alpha }^{(a)}(t)$ and $W_{\hspace{0.33em}En}^{(a)}(t)$. The results of removing local trends and determining local extremum points are presented in Figure 9 (b1, b2, b3). The number 433 of most prominent local extrema was chosen equal to the number of earthquakes with magnitude not less than 7.

Let's consider how the time points of the most significant local extrema of the statistics $W_{\hspace{0.33em}\gamma }^{(a)}(t)$, $W_{\hspace{0.33em}\Delta \alpha }^{(a)}(t)$ and $W_{\hspace{0.33em}En}^{(a)}(t)$ are distributed. To do this, let’s calculate the probability histograms separately for 433 smallest local minima $W_{\hspace{0.33em}\gamma }^{(a)}(t)$, $W_{\hspace{0.33em}\Delta \alpha }^{(a)}(t)$ and 433 largest local maxima of entropy $W_{\hspace{0.33em}En}^{(a)}(t)$ (it means case $({\gamma }_{\min },\Delta {\alpha }_{\min },E{n}_{\max })$), i.e. empirical probabilities of time values falling within successive six-month intervals. Since the total length of observations is 29 years, the number of such intervals is 58. Figures 10(a1, b1, c1) show graphs of these histograms. Figures 10(a2, b2, c2) correspond to probability histograms for 433 largest local maxima $W_{\hspace{0.33em}\gamma }^{(a)}(t)$, $W_{\hspace{0.33em}\Delta \alpha }^{(a)}(t)$ and 433 smallest local minima of entropy $W_{\hspace{0.33em}En}^{(a)}(t)$ (case $({\gamma }_{\max },\Delta {\alpha }_{\max },E{n}_{\min })$). It shows that the last six-month interval has the highest empirical probability of containing local extremum points except Figure 10(a2).

Thus, from the graphs in Figure 10 it is evident that the times of prominent local extremes of almost all seismic noise statistics are significantly concentrated in the second half of 2025.

8. Influence Matrices

Further analysis involves calculating a measure describing the lead or lag of two time sequences, one of which represents the times of strong earthquakes, and the other represents the times of certain characteristic features of seismic noise statistics. Such features could be, for example, the times of the most prominent local extrema.

To solve this problem, we will apply the influence matrices method, which was used in [42,43,44,45] to analyze the prognostic properties of the earth's surface tremor measured by GPS, to analyze the relationship between magnetic field fluctuations and strong earthquakes, and to analyze the relationship between anomalies in meteorological time series and seismicity.

Let $t_{\hspace{0.17em}j}^{(p)},\hspace{0.33em}j=1,\mathrm{...},N_p;\hspace{0.33em}p=1,2$ represent the moments of time of 2 sequences of events. Let us represent their intensities as follows:

$$\hspace{0.33em}{\mu }^{(p)}(t)=b_{\hspace{0.17em}0}^{(p)}\hspace{0.17em}+{\displaystyle \sum _{q=1}^2{b}_{\hspace{0.17em}q}^{(p)}}\cdot g^{(q)}(t)$$

(15)

where $b_{\hspace{0.17em}0}^{(p)}\ge 0,b_q^{(p)}\ge 0$ are parameters, $g^{(q)}(t)\hspace{0.17em}$ - function of influence of time moment $t_j^{(q)}$ of the sequence with number $q$:

$$g^{(q)}(t)={\displaystyle {\sum }_{t_j^{(q)}<t}{e}^{-(t-t_{\hspace{0.33em}j}^{(q)})/\tau }}$$

(16)

According to formula (16), the weight of the event with number $j$ becomes non-zero for times $t\hspace{0.17em}>t_j^{(q)}$ and decays with characteristic time $\tau$. The parameter $b_q^{(p)}$ determines the degree of influence of the flow $q$ on the flow $p$. The parameter $b_p^{(p)}$ determines the degree of influence of the flow $p$ on itself (self-excitation), and the parameter $b_{\hspace{0.17em}0}^{(p)}$ reflects a purely random (Poisson) component of intensity. Let us fix the parameter $\tau$ and consider the problem of determining the parameters $b_{\hspace{0.17em}0}^{(p)},b_q^{(p)}$.

The log-likelihood function for a non-stationary Poisson process is equal to over the time interval $[0,T]$ [40]:

$$\ln (L_p)={\displaystyle \sum _{j=1}^{N_p}\ln ({\mu }^{(p)}(t_j^{(p)}))-{\displaystyle \underset{0}{\overset{T}{\int }}{\mu }^{(p)}(s)ds}},\hspace{1em}p=1,2$$

(17)

It is necessary to find the maximum of functions (17) with respect to the parameters $b_{\hspace{0.17em}0}^{(p)},b_q^{(p)}$. As shown in detail in [13,42,43,44], the problem of finding parameter values $b_{\hspace{0.17em}0}^{(p)},b_q^{(p)}$ from the maximum likelihood principle is reduced to solving problems for the maximum:

$${\displaystyle \sum _{j=1}^{N_p}\ln ({\mu }_{\hspace{0.17em}0}^{(p)}+{\displaystyle \sum _{q=1}^2{b}_q^{(p)}\cdot \Delta g^{(q)}(t_j^{(p)}))\to }}\underset{b_1^{(p)},b_2^{(p)}}{\max },\hspace{1em}p=1,2$$

(18)

where ${\mu }_{\hspace{0.17em}0}^{(p)}=N_p/T$, $\Delta g^{(\beta )}(t)=g^{(\beta )}(t)-{\overline{g}}^{(\beta )}$, ${\overline{g}}^{(q)}={\displaystyle {\int }_0^T{g}^{(q)}(s)ds}/T$, under restrictions:

$$b_1^{(p)}\ge 0,\hspace{1em}{b}_2^{(p)}\ge 0,\hspace{1em}{b}_1^{(p)}{\overline{g}}^{(1)}+b_2^{(p)}{\overline{g}}^{(2)}\le {\mu }_{\hspace{0.17em}0}^{(p)}$$

(19)

Having solved problem (18-19) numerically for a given $\tau$, we can introduce the elements of the influence matrix ${\kappa }_q^{(p)},p=1,2;\hspace{0.33em}q=0,1,2$ according to the formulas:

$${\kappa }_0^{(p)}=b_{\hspace{0.17em}0}^{(p)}/{\mu }_{\hspace{0.17em}0}^{(p)}\ge 0,\hspace{1em}{\kappa }_q^{(p)}=b_q^{(p)}\cdot {\overline{g}}^{(q)}/{\mu }_{\hspace{0.17em}0}^{(p)}\ge 0$$

(20)

The quantity ${\kappa }_0^{(p)}$ is a part of the average intensity ${\mu }_{\hspace{0.17em}0}^{(p)}$ of the process with number $p$, which is purely random, the part ${\kappa }_p^{(p)}$ is caused by the influence of self-excitation and ${\kappa }_q^{(p)},\hspace{0.33em}p\ne q$ is determined by the external influence $q\to p$. From formula (23) follows the normalization condition:

$${\kappa }_0^{(p)}+{\kappa }_1^{(p)}+{\kappa }_2^{(p)}=1,\hspace{1em}p=1,2$$

(21)

As a result, we can determine the influence matrix:

$$\left(\begin{array}{l}{\kappa }_{\hspace{0.33em}\hspace{0.17em}0}^{(1)}\hspace{0.33em}{\kappa }_{\hspace{0.17em}\hspace{0.17em}1}^{(1)}\hspace{0.33em}{\kappa }_{\hspace{0.17em}\hspace{0.17em}2}^{(1)}\\ {\kappa }_{\hspace{0.17em}\hspace{0.17em}0}^{(2)}\hspace{0.33em}{\kappa }_{\hspace{0.33em}1}^{(2)}\hspace{0.33em}{\kappa }_{\hspace{0.33em}2}^{(2)}\end{array}\right)$$

(22)

The first column of matrix (22) is composed of Poisson shares of mean intensities. The diagonal elements of the right submatrix of size 2×2 consist of self-excited elements of mean intensity, while the off-diagonal elements correspond to mutual excitation. The sums of the component rows of the influence matrix (22) are equal to 1. The influence matrices are estimated in a certain sliding time window of length with offset and with a given value of the attenuation parameter $\tau$.

We will investigate the combined effect of the advance or lag of local extrema of the functions $W_{\hspace{0.33em}\gamma }^{(a)}(t)$, $W_{\hspace{0.33em}\Delta \alpha }^{(a)}(t)$ and $W_{\hspace{0.33em}En}^{(a)}(t)$ relative to the earthquake moments. To do this, we will calculate the first principal component of these functions in a sliding time window of 365 days according to formulas (1-6) for the P-averaging and denote it by $P_{\gamma ,\Delta \alpha ,En}^{\hspace{0.33em}\hspace{0.33em}(Wa)}(t)$.

The components of the influence matrix (22) will be calculated in a sliding time window of 2.6 years with an offset of 0.01 years. This choice of window length is due to the fact that it is equal to the dominant period of the sequence of strong earthquakes (Figure 8).

The parameter $\tau$ is searched from maximizing the average value of the difference $\Delta \kappa ={\kappa }_{\hspace{0.17em}\hspace{0.17em}2}^{(1)}-{\kappa }_{\hspace{0.33em}1}^{(2)}\to \underset{\tau }{\max }$, where ${\kappa }_{\hspace{0.17em}\hspace{0.17em}2}^{(1)}$ is the portion of the intensity of the earthquake sequence caused by the "influence" of the sequence of local extrema, and ${\kappa }_{\hspace{0.33em}1}^{(2)}$ is the portion of the intensity of the sequence of local extrema caused by the leading influence of seismic events.

The graphs in Figure 11 show that local minima of the function $P_{\gamma ,\Delta \alpha ,En}^{\hspace{0.33em}\hspace{0.33em}(Wa)}(t)$ lead earthquake times significantly more than vice versa (average lead measures are 0.714 versus 0.228). Moreover, the lead measure of local minima $P_{\gamma ,\Delta \alpha ,En}^{\hspace{0.33em}\hspace{0.33em}(Wa)}(t)$ relative to seismic event times is subject to cyclic variations with a period of about 9 years, as can be seen from Figure 11(b2). The cyclic trend in Figure 11(b2) (purple line) is determined by minimizing the variance of the residual using a trial cyclic trend with fitted parameters of the linear trend and oscillation period. As for the local maxima of the function $P_{\gamma ,\Delta \alpha ,En}^{\hspace{0.33em}\hspace{0.33em}(Wa)}(t)$, they almost always occur after earthquakes, and the reverse effect is negligible (mean lag measure is 0.802 versus 0.023).

9. The Step Component of the Entropy Probability Density and Its Relationship with the Strongest Earthquakes

To quantitatively describe the sequence of values of the weighted probability density function (PDF) for the distribution of extreme seismic noise properties in 10-day time windows, we consider the entropy of these densities. Let $M=100\times 50=50000$ be the number of grid nodes. We calculate the histogram of the weighted probability density functions for each 10-day interval $i$ as a sequence of values:

$$g_{\hspace{0.33em}k}^{(i)}=l_{\hspace{0.33em}k}^{(i)}/L_b,\hspace{1em}k=1,\mathrm{...},L_b$$

(23)

where $L_b$ is the number of intervals dividing the range of density variation within a given 10-day window into intervals of equal length, $l_{\hspace{0.33em}k}^{(i)}$ is the number of PDF values falling within the interval with number $k$. According to recommendations [34,41], we choose the value $L_b=\sqrt{M}$. Obviously ${\displaystyle {\sum }_k{g}_{\hspace{0.33em}k}^{(i)}=1}$, the values of the histogram (23) also have the properties of a discrete probability distribution. We define the normalized entropy of this distribution by the formula:

$${{\rm E}}_i=-{\displaystyle \sum _{k=1}^{L_b}{g}_{\hspace{0.33em}k}^{(i)}\ln (g_{\hspace{0.33em}k}^{(i)})/\ln (L_b),\hspace{1em}0\le {{\rm E}}_i\le 1}$$

(24)

In Figure 12(a), the gray line represents the graph of entropy change (24) depending on the position of the right end of 10-day time windows. The blue line represents a stepwise approximation using so-called long WTMM (Wavelet Transform Modulus Maxima) chains of the signal skeleton. Scale-dependent stepwise approximation allows one to formally determine the moments in time of significant changes in the mean entropy values. This approximation is constructed using continuous wavelet transforms with a kernel in the form of derivatives of a Gaussian function, which allows one to formalize the selection of moments in time at which a significant, scale-dependent change in the mean value of a noisy signal occurs [46,47].

Let us define smoothed functions of entropy (24):

$$\overline{{\rm E}}(t,a)={\displaystyle \sum _{\xi }{{\rm E}}_{\xi }\cdot e^{-{\left(\frac{\xi -t}{a}\right)}^2}}/{\displaystyle \sum _{\xi }{e}^{-{\left(\frac{\xi -t}{a}\right)}^2}}$$

(25)

In formula (25) a Gaussian is used as the smoothing kernel, the parameter $a>0$ being the smoothing radius. Due to the rapid decay of the Gaussian and its symmetry, the first derivative of the smoothed entropy can be written in the form:

$${\overline{{\rm E}}}^{\prime }(t,a)=2\cdot {\displaystyle \sum _{\xi }{{\rm E}}_{\xi }\cdot \left(\frac{\xi -t}{a}\right)\cdot e^{-{\left(\frac{\xi -t}{a}\right)}^2}}/{\displaystyle \sum _{\xi }{e}^{-{\left(\frac{\xi -t}{a}\right)}^2}}$$

(26)

We define a WTMM point $t_a$ as the argument of the modulus $|{\overline{{\rm E}}}^{\prime }(t,a)|$ at which a local maximum is achieved with respect to $t$. As the smoothing parameter $a$ changes, the points $t_a$ on the plane $(t,a)$ form so-called WTMM chains. The set of all WTMM chains forms the WTMM skeleton of the signal. Since ${\overline{{\rm E}}}^{\prime }(t,a)=d\overline{{\rm E}}(t,a)/dt$, the points $t_a$ determine the time instants for which the maximum trends (decrease or increase) of the smoothed function $\overline{{\rm E}}(t,a)$ are observed. If we set a parameter $a=a^{*}$ , a sequence of WTMM points $t_{a^{*}}^{(k)}$ will arise, using which we can formalize a stepwise approximation of the original signal ${\rm E}(t)$ as a sequence of its average values between adjacent time instants $t_{a^{*}}^{(k)}$ and $t_{a^{*}}^{(k+1)}$.

Let us consider the sequence of time instants of stepwise changes in the WTMM approximation of entropy (24) of the values of the weighted average density distribution of extreme values $({\gamma }_{\min },\Delta {\alpha }_{\min },E{n}_{\max })$. Figure 12(a) shows a graph of the change in this entropy (gray line) and a graph of the stepwise WTMM approximation. Since entropy (24) is defined as a sequence of values with a time step of 10 days, then when selecting a parameter $a^{*}=10$, the dimensional length of this parameter is equal to 100 days.

Figure 12(b) shows the time sequence of absolute values of steps in the WTMM approximation. We examine the relationship between steps in the WTMM approximation of entropy and the time instants of the strongest earthquakes with a magnitude of at least 7.8, which is shown in Figure 12(c).

Let's calculate the influence matrix between the two time sequences shown in Figure 12(b) and Figure 12(c) over the entire length of the available sample. The result depends on the choice of parameter $\tau$, which we select based on the condition of maximizing the difference between the influence matrix components responsible for mutual influence. This value is $\tau =0.45$ years. The calculation results are presented in the Table 1.

The results presented in the Table indicate that the abrupt changes in the entropy values of the weighted average probability densities of the extreme value distributions for the case $({\gamma }_{\min },\Delta {\alpha }_{\min },E{n}_{\max })$ significantly precede the times of the strongest earthquakes (the lead measure is 0.382 versus 0.027). It is interesting to note that, when similarly considering the "antagonistic" case $({\gamma }_{\max },\Delta {\alpha }_{\max },E{n}_{\min })$, the mutual influence is negligible.

10. Estimates of Coherence

For the further analysis it is necessary to estimate coherence spectra between 2 time series within sliding time window. A parametric model of vector autoregression will be used which provides a better frequency resolution with respect to methods of spectra and cross-spectra estimates based on Fourier expansion [46]. For the multiple time series $X(t)$ of dimensionality $q$ AR-model is defined by a formula:

$$X(t)+{\displaystyle \sum _{k=1}^p\hspace{0.17em}{B}_{\hspace{0.33em}k}\cdot }X(t-k)=\epsilon (t)\hspace{0.33em}$$

(27)

Here $t$ is a discrete time index, $p$ is autoregression order, $B_{\hspace{0.33em}k}$ are matrices of autoregression coefficients of the size $q\times q$, $P=M(\epsilon (t){\epsilon }^T(t))$ is covariance matrix of the size $q\times q$ of residual signal $\epsilon (t)$. Matrices $B_{\hspace{0.33em}k}$ and $P$ are computed by Durbin-Levinson procedure [46]. Parametric estimate of spectral matrix is defined by the formula:

$$S_{XX}(\omega )={\Phi }^{-1}(\omega )\cdot P\cdot {\Phi }^{-H}(\omega ),\hspace{1em}\Phi (\omega )=E+{\displaystyle \sum _{k=1}^p{B}_{\hspace{0.33em}k}}{e}^{-i\omega k}$$

(28)

where $E$ is a unit matrix of the size $q\times q$, $\omega$ is a frequency. For dimensionality $q=2$ squared coherence spectrum is computed by the formula:

$$c^2(\omega )=|S_{12}(\omega ){|}^2/S_{11}(\omega )\cdot S_{22}(\omega )$$

(29)

Here $S_{11}(\omega )$ and $S_{22}(\omega )$ are diagonal elements of the matrix (28) whereas $S_{12}(\omega )$ is cross-spectrum.

11. Connection of Seismic Noise Response to Irregularity of Earth Rotation with Major Earthquakes

Some researchers pointed to the connection between the uneven rotation of the Earth and seismicity [49]. The trigger mechanism of the effect of variations in the Earth's rotation on the seismic process was studied in [50]. Estimates of the effect of a strong earthquake on the length of the day are given in [51]. The graph 13(a) presents length of day time series for interval 1997-2025 [52].

Relations between seismic noise properties and irregularity of Earth rotations was investigated in [13,53]. We calculate maximum quadratics coherence spectrum between LOD and first principal component of properties $P_{\gamma ,\Delta \alpha ,En}^{\hspace{0.33em}\hspace{0.33em}(j)}(t)$ for each of 50 reference points in a moving time window of the length 365 days with mutual shift 3 days. We will call these maximum quadratic coherences as response $R_{ \gamma ,\Delta \alpha ,En}^{\hspace{0.33em}\hspace{0.33em}(j)}(t)$ of seismic noise to LOD. We used 2-dimensional AR-model (34) for $p=5$ within moving time window of the length 365 days with mutual shift 3 days.

Examples of such response functions are presented at Figure 13 (b1)-(b6) for 6 reference points. Having response functions $R_{ \gamma ,\Delta \alpha ,En}^{\hspace{0.33em}\hspace{0.33em}(j)}(t)$ from all reference points we calculated their weighted mean value where we used squared components of the first eigenvector (corresponding to maximum eigenvalue) of the 50×50 size covariance matrix. This weighted mean of all response functions we will call integral response of 3 seismic noise properties all over the world to irregularity of Earth rotation.

Graph 13(c) of the integrated response of seismic noise properties to the Earth's rotation irregularity shows that the strongest earthquakes $M\ge 8.5$ occur primarily after the response reaches its maximum. This response behavior highlights a distinctive pattern before and after the Kamchatka mega-earthquake of July 29, 2025, $M=8.8$: a prolonged increase prior to the event and a sharp decline immediately afterward.

To analyze the relationship between local extrema of the integrated seismic noise response to LOD and seismic events $M\ge 7$, we remove local trends from the response function using a Gaussian window. The goal of this operation, as before, is to isolate the effect of the leading time of the most prominent local extrema relative to the earthquake time points. However, it turned out that the optimal window width, for which the leading effect in this case is maximal, is 24 days (8 time steps with a step length of 3 days). Figure 14(a) shows the response function graphs after removing local trends and the positions of the 433 largest local maxima. Graphs 14(b1) and 14(b2) show graphs of the change in the influence matrix components responsible for the mutual leading of the analyzed time sequences. The assessments were performed in a sliding time window of 2.6 years with an offset of 0.01 years. The damping parameter value $\tau$ was chosen to maximize the difference between the "forward" and "reverse" lead measures. This value turned out to be equal to $\tau =0.148$ year.

Comparing the graphs in Figure 14(b1) and Figure 14(b2) shows how much the lead measure of the local maxima of the LOD response function relative to earthquake times outperforms the inverse measure (average value of 0.253 versus 0.049). Figure 14(b2) shows that the lead measure is subject to a cyclical trend with a period of 7.5 years. This period was found by minimizing the residual variance for cyclical trends with trial periods.

12. Discussion

A method for analyzing long-term continuous seismic noise records from a global network of broadband seismic stations is proposed. This method utilizes their multifractal properties and entropy estimates based on wavelet decompositions. The goal of this research is to establish patterns in the evolution of noise properties and their relationships to the sequence of strong earthquakes and the irregularity of the Earth's rotation.

As a result of the analysis, the following facts were established:

1. From 2003 to the end of 2025, spatial correlations in global seismic noise properties increased. This increase follows a linear trend, superimposed by oscillations with a period of approximately 2.5 years – Figure 3(a).

2. Since 2011, following the Tohoku mega-earthquake of March 11, 2011 in Japan, there has been a sharp increase in the largest distances at which strong correlations between seismic noise properties occur in a set of 50 control points covering the surface of the globe – Figure 3(c).

3. Along with the growth of the largest distances of strong correlations, the number of pairs of control points between which strong correlations arise also increases. In the second half of 2025, following the Kamchatka mega-earthquake of July 29, 2025 (M=8.8), the number of pairs of control points with strong correlations exploded, tripling over six months – Figure 3(b), Figure 4.

4. On the Earth’s surface there are stable areas of concentration of the highest values of probability densities of extreme values of seismic noise properties – Figure 5, Figure 6.

5. The sequence of strong earthquakes with a magnitude of at least 7 demonstrates a complex non-stationary nature of the evolution of its periodic components (an increase in frequency before 2014 and a decrease after 2018) with the identification of a dominant period of 2.6 years, close to the period of 2.5 years of oscillations of the average value of pairwise correlations between the properties of seismic noise at reference points – Figure 8.

6. Probability histograms of time points at which extreme values of seismic noise properties are observed have significant maxima in the second half of 2025 – Figure 10.

7. The average lead measure for the smallest local minima of the Donoho-Johnstone index and the singularity spectrum support width and the largest local maxima of the entropy values of the wavelet coefficients relative to the time moments of earthquakes with a magnitude of at least 7 significantly leads the inverse lead measure (0.714 versus 0.228) – Figure 11.

8. The full measure of the lead time of the step changes in the entropy of the values of the average weighted distribution densities of the minimum values of the Donoho-Johnston index and the the singularity spectrum support width and the maximum values of the entropy of the wavelet coefficients relative to the moments of time of earthquakes with a magnitude of at least 7.8 significantly leads the inverse measure of the lead (0.382 versus 0.000) – Figure 12, Table 1.

9. The average lead measure for the largest local maxima of the value of the generalized response of seismic noise properties to the unevenness of the Earth's rotation relative to the moments of time of earthquakes with magnitudes of at least 7 significantly leads the inverse lead measure (0.253 versus 0.049) – Figure 14.

13. Conclusions

Taken together, the established facts allow us to assert that three important features, which can be termed "events," occurred in the Earth's dynamics, as manifested by low-frequency seismic noise, between 1997 and 2025. The first "event" occurred around 2003, when the process of increasing global spatial correlation of seismic noise began. The second "event" occurred around 2013, when the increase in spatial correlation of noise increased essentially. Finally, the third "event" occurred in the second half of 2025, when spatial correlation of noise increased dramatically. The final "event" can be associated with the Kamchatka mega-earthquake of July 29, 2025.

The results also indicate that the minimum values of the Donoho-Johnston index and the multifractal singularity spectrum support width, as well as the maximum values of the entropy of the wavelet coefficient distribution, possess significant predictive power for the sequence of strong earthquakes. Furthermore, the maximum values of the response of seismic noise properties to the irregularity of the Earth's rotation also possess predictive power.