Low-Frequency Seismic Noise Properties in the Japanese Islands

The records of seismic noise in Japan for the period of 1997–2020, which includes the Tohoku seismic catastrophe on 11 March 2011, are considered. The following properties of noise are analyzed: The wavelet-based Donoho–Johnston index, the singularity spectrum support width, and the entropy of the wavelet coefficients. The question of whether precursors of strong earthquakes can be formulated on their basis is investigated. Attention is paid to the time interval after the Tohoku mega-earthquake to the trends in the mean properties of low-frequency seismic noise, which reflect the constant simplification of the statistical structure of seismic vibrations. Estimates of two-dimensional probability densities of extreme values are presented, which highlight the places in which extreme values of seismic noise properties are most often realized. The estimates of the probability densities of extreme values coincide with each other and have a maximum in the region: 30° N ≤ Lat ≤ 34° N, 136° E ≤ Lon≤ 140° E. The main conclusions of the conducted studies are that the preparation of a strong earthquake is accompanied by a simplification of the structure of seismic noise. It is shown that bursts of coherence between the time series of the day length and the noise properties within annual time window precede bursts of released seismic energy. The value of the lag in the release of seismic energy relative to bursts of coherence is about 1.5 years, which can be used to declare a time interval of high seismic hazard after reaching the peak of coherence.


Introduction
The article is devoted to the processing of seismic noise data recorded at the network F-net of stations on the Japanese Islands for 24 years, 1997-2020. During this time period, on 11 March 2011, a mega-earthquake with a magnitude of 9.1 occurred in Japan. Japan is a region with a dense, open-access seismic network F-net provided by the National Research Institute for Earth Science and Disaster Resilience (NIED). This makes it possible to test various hypotheses about how the preparation of a strong seismic event can affect the statistical properties of seismic noise.
In works [1][2][3][4][5][6], it was shown that the processes of earthquake preparation are preceded by certain changes in the statistical structure of seismic noise. The main changes are in the simplification of noise, namely in the growth of entropy and loss of multifractality. The main sources of seismic noise energy are not earthquakes, but the movements of cyclones in the atmosphere and the impact of ocean waves on the shelf and coast [7][8][9]. Thus, the sources of noise energy are located outside the earth's crust. However, the crust is a medium for the propagation of seismic waves. As a result, processes inside the earth's crust, including the preparation of strong seismic events, are reflected in changes in the properties of seismic noise.
The article investigates three properties of seismic noise: Entropy, determined through the distribution of the squares of orthogonal wavelet coefficients, the multifractal singularity spectrum support width, and the Donoho-Johnstone (DJ) index, defined as the proportion of the total number of wavelet coefficients that can be considered as "informa-tion carriers." These properties are estimated daily from seismic noise records at a network of stations. The choice of these parameters is due to the fact that their changes reflect the complication or simplification of the noise structure. With the simplification of the structure, the entropy increases, while the support width and the DJ index decrease. The methods outlined in the article are designed to detect spatial areas and time intervals, where and when the simplification of the noise structure is observed. These areas and time intervals are interpreted as manifestations of increased seismic hazard.
In addition, a possible trigger of the unevenness of the Earth's rotation in relation to the release of seismic energy is being investigated. For this, the quadratic coherence between the length of the day and the first principal component of the analyzed properties of seismic noise is estimated within the annual sliding time window. The correlation function between bursts of coherence and the release of seismic energy turned out to be significantly shifted by time delays corresponding to the advance of the coherence maxima to strong earthquakes.

Description of Seismic Noise Properties
Minimum normalized entropy of wavelet coefficients En. Let x(t) be a random signal, and t = 1, . . . , N is the discrete time index. Let c (k) j be the wavelet coefficients of the analyzed signal. The superscript k is the number of the detail level of the orthogonal wavelet decomposition and the subscript j numbers the sequence of centers of time intervals in the vicinity of which the convolution of the signal with finite elements of the basis is calculated. The bases of 17 orthogonal Daubechies wavelets were used: 10 ordinary bases with a minimum support with the number of vanished moments from 1 to 10 and 7 Daubechies symlets [10] with the number of vanishing moments from 4 to 10. For each of the bases, the normalized entropy of the distribution of the squares of the coefficients was calculated. (1) where m is the number of levels of detail accepted for consideration; M k is the number of wavelet coefficients at the detail with the number k. The number of levels m depends on the length N of the analyzed sample. For example, if N = 2 m , then m = n, M k = 2 (n−k) . The condition N = 2 m is necessary to apply the fast wavelet transform. If the length of N is not equal to a power of two, then the signal x(t) is padded with zeros to the minimum length L, which is greater than or equal to N: L = 2 m ≥ N. In this case, among the number M k = 2 (n−k) of all wavelet coefficients at the level k, only N · 2 −k coefficients correspond to the decomposition of the real signal, while the remaining coefficients are equal to zero due to the addition of zeros to the signal x(t). Thus, in formula (1) M k = N · 2 −k and only "real" coefficients c (k) j are used to calculate the entropy. The number N r in formula (1) is equal to the number of "real" coefficients, that is, N r = ∑ m k=1 M k . The optimal wavelet basis is found as those that provide minimum to entropy among all 17 tested bases: En → min . By construction, 0 ≤ En ≤ 1.
The minimum entropy (1) was proposed in [11]. In works [3,5,6,[11][12][13], it was used to analyze the predictive properties of seismic noise. Entropy (1) has common features with the multiscale entropy considered in [14,15]. The multiscale follows from the use of the wavelet transform, which provides decomposition into discrete energetic atoms c (k) j 2 . A similar entropy construction was proposed in [16][17][18] for the natural time method. In [19,20], the non-extensive Tsallis entropy was used to study seismic noise.
Donoho-Johnstone index γ. In wavelet filtering, there is a procedure of "thresholding", which means setting wavelet coefficients that are smallest by the absolute values to zero [10,21]. An assumption that noise is mainly accumulated at the first detail level is accepted. The variance of the wavelet coefficients is equal to the variance of the initial signal following from orthogonality of the wavelet transform. Thus, the noise standard deviation σ could be calculated as standard deviation of the wavelet coefficients at the first detail level. In order to provide a robust estimate, the median of absolute wavelet coefficients at the first detail level is used σ = med c (1) k , k = 1, . . . , N/2 /0.6745 (2) where N/2 is the number of these coefficients. The formula (2) follows from the relation between median and standard deviation for normal random value: Med ≈ 0.6745 · σ. The estimate of standard deviation σ from formula (2) determines the quantity σ √ 2 · ln N as a "natural" threshold for separating the noisy wavelet coefficients [10,21]. The formula σ √ 2 · ln N is following from asymptotic probability of maximal deviations of Gaussian white noise [10]. Let us define the index γ, 0 < γ < 1 as the ratio of number of wavelet coefficients satisfying inequality c k > σ √ 2 · ln N to the total number N of wavelet coefficients. The index γ is calculated for the optimal wavelet basis, which was found from minimum of entropy (1). The lower the index γ, the noisier the signal.
Singularity spectrum support width ∆α. The measure of variability µ X (t, δ) of the random signal x(t) on the time interval [t, t + δ] is defined as its range: x(u). Let us consider the mean value of its power de- The signal is scale-invariant [22] if M(δ, q) ∼ δ ρ(q) when δ → 0 , that is, the following limit exists: The process is monofractal if ρ(q) = Hq, where H = const, 0 < H < 1; otherwise, when ρ(q) is a nonlinear concave function of q, the signal is multifractal. The value of ρ(q) for a finite sample x(t), t = 0, 1, . . . , N − 1 could be calculated using the method of detrended fluctuation analysis (DFA) [23,24], which was modified for estimating multifractals in [25]. The time series is split into adjacent intervals of length s: Let us consider a part of signal x(t), corresponding to interval I (s) k : Let us fit a polynomial of the order m p (s, m) k (t) to the signal y (s) k (t) and consider the deflections: ∆y and the sum: The quantity (7) could be regarded as the estimate of (M(δ s , q)) 1/q . Let us introduce the function h(q) as a linear regression coefficient between ln(Z (m) (q, s)) and ln(s): Z (m) (q, s) ∼ s h(q) within scales ranging s min ≤ s ≤ s max . The minimum value of scale s min within Formulae (4)-(7) was chosen to be 20 samples, and the maximum scale equals s max = N/5. For the monofractal signal h(q) = H = const, but in the general case, ρ(q) = qh(q). The multifractal singularity spectrum F(α) is defined as the fractal dimen-sionality of the set of time moments t for which the Hölder-Lipschitz exponent is equal to α, which means |x(t + δ) − x(t)| ∼ |δ| α , δ → 0 [26]. Let us calculate a Gibbs sum: The mass exponent τ(q) is defined by the condition W(q, s) ∼ s τ(q) . A formula τ(q) = ρ(q) − 1 = qh(q) − 1 follows from (7). The values of exponent q in the Formula (7) were taken from interval [−Q, +Q] where Q is some large number. The value Q = 10 is used. The values F(α) = min +Q] dτ(q)/dq. The derivative dτ(q)/dq is calculated numerically. The accuracy of its calculation is not very important, because this derivative is used for a rough determination of an a priori interval of α values. The value of α min and α max are determined as minimum and maximum values of α for which F(α) ≥ 0. Thus, the spectrum F(α) is defined according to the formula: Let us consider the estimates of singularity spectrum F(α) in a sliding window. For this case, its evolution can provide important information about the structure of the chaotic pulsations of the series. The support width of the singularity spectrum ∆α = α max − α min is an important characteristic of the signal and it is regarded as a measure of variety of stochastic behavior.
Multifractal analysis is often used in geophysics [27,28]. The natural time approach [17] has its own toolboxes using multifractals and multiscale entropy for the analysis of seismicity [18,[29][30][31]. The multifractal property ∆α of seismic noise was used for the purposes of predicting earthquakes and assessing seismic hazard in [1][2][3][4][5]13,32]. The singularity spectrum support width ∆α is used to study the behavior of various nonlinear systems. A decrease in the parameter ∆α is a well-known effect that anticipates changes in the properties of biological and medical systems [33][34][35]. It was shown in [36] that the "loss of multifractality" is also universal in nature in physical systems.

Seismic Noise Data
Vertical seismic oscillations data with sampling rate 1 Hz were used for the analysis. These data are accessible from the source [37] for 78 seismic stations of the network F-net in Japan. For the analysis, a time interval from the beginning of 1997 up to 31 March 2021 was selected. Figure 1 presents positions of the network stations and the location of the Nankai Deep Trench, which is the northern boundary of the Philippine tectonic plate.
The seismic data with a sampling frequency of 1 Hz were reduced to a time step of 1 minute by calculating the mean values in adjacent time intervals of 60 values. The seismic records from each station after coming to a 1-minute sampling time step were split into time fragments of the length of 1 day (1440 samples) and for each fragment, parameters (En, γ, ∆α) of daily seismic noise waveforms were calculated. Scale-dependent trends in the Formulae (6) and (7) were removed by polynomials of the eighth order. Removing trends from seismic noise waveforms by the polynomial of the eighth order was used before computing entropy En and index γ in each daily time window. Thus, the time series of (En, γ, ∆α) values with a sampling time step of 1 day was obtained from each of the seismic stations, which are presented in Figure 1. The seismic data with a sampling frequency of 1 Hz were reduced to a time step of 1 minute by calculating the mean values in adjacent time intervals of 60 values. The seismic records from each station after coming to a 1-minute sampling time step were split into time fragments of the length of 1 day (1440 samples) and for each fragment, parameters ( , , ) En   of daily seismic noise waveforms were calculated. Scale-dependent trends in the Formulae (6) and (7) were removed by polynomials of the eighth order. Removing trends from seismic noise waveforms by the polynomial of the eighth order was used before computing entropy En and index  in each daily time window. Thus, the time series of ( , , ) En   values with a sampling time step of 1 day was obtained from each of the seismic stations, which are presented in Figure 1. Figure 2 presents the graph of working stations each day during the considered time interval. The seismic station is considered working within the current day if there are no gaps during this day. One can see that the number of operable stations at the initial time fragment (1997-2001) is rather small. This influences the reliability of further estimates. Nevertheless, the result of data processing should not be removed from the analysis. Particularly, median values of seismic noise properties are rather stable to the number of working stations.    The seismic data with a sampling frequency of 1 Hz were reduced to a time step of 1 minute by calculating the mean values in adjacent time intervals of 60 values. The seismic records from each station after coming to a 1-minute sampling time step were split into time fragments of the length of 1 day (1440 samples) and for each fragment, parameters ( , , ) En   of daily seismic noise waveforms were calculated. Scale-dependent trends in the Formulae (6) and (7) were removed by polynomials of the eighth order. Removing trends from seismic noise waveforms by the polynomial of the eighth order was used before computing entropy En and index  in each daily time window. Thus, the time series of ( , , ) En   values with a sampling time step of 1 day was obtained from each of the seismic stations, which are presented in Figure 1. Figure 2 presents the graph of working stations each day during the considered time interval. The seismic station is considered working within the current day if there are no gaps during this day. One can see that the number of operable stations at the initial time fragment (1997-2001) is rather small. This influences the reliability of further estimates. Nevertheless, the result of data processing should not be removed from the analysis. Particularly, median values of seismic noise properties are rather stable to the number of working stations.  The initial seismic records with a 1 Hz sampling rate were taken "as is", which includes background seismic oscillations and records immediately after earthquakes. It is important to underline that these 1 Hz records were transformed into a sampling time step of 1 min by averaging within adjacent time fragments of the length of 60 values. This smoothing operation gets rid of the influence of high-frequency immediate reactions of earthquakes.
The issue of predicting of strong earthquakes in Japan using entropy and multifractal properties of seismic noise was investigated in [1][2][3][4][5]. The natural time approach was applied to estimate seismic danger in Japan in [38,39].

Averaged Maps of Seismic Noise Properties
Let us consider the regular grid of the size 30 × 30 nodes covering the region with latitudes between 28 • N and 46 • N and longitudes between 128 • E and 146 • E (see Figure 1). Let U be any value of ∆α, En or γ. For each node (i, j) of the grid and for each day t, the five nearest operable seismic stations are found, which provides five values of U. Let us take a median U ij form an "elementary" daily map. These daily maps could be averaged: for daily time indexes t between two given dates t 0 and t 1 . The averaged map of ∆α before the Tohoku earthquake is presented in Figure 3a2. One can notice that the region of future seismic events is extracted by relatively low values of ∆α. If Figure 3a1,a2 are compared, one can see that after the earthquake on 25 September 2003, the domain with low values of singularity spectrum support width was split into two parts. The northern part turns out to be the region of the mega-earthquake on 11 March 2011 whereas the southern part preserves low values of ∆α (Figure 3a2,a3).
Based on the assumption that low values of ∆α correspond to high seismic danger, a hypothesis that the Tohoku earthquake drops only part of the accumulated tectonic energy from the southern region could be considered, and the region corresponding to the Nankai Trough could be regarded as the region of a future strong earthquake. Comparing the averaged maps of entropy En in Figure 3b1-b3 with similar maps of ∆α in Figure 3a1-a3, one can notice that they are antipodes of each other. It means that high values of entropy En correspond to regions with high seismic danger. It should be noted that "high" and "low" values of seismic noise statistics are not considered absolute, but rather relative to the values within different time intervals. Each time interval is characterized by its own mean value, and differences with respect to mean value within the time interval define whether the value is "high" or "low" in the considered time interval.
A possible mechanism as to why low values of ∆α and high values of En extract seismically dangerous regions was given in [2][3][4][5][6]. It is considered as the consequence of consolidation of small blocks of the Earth's crust into the large one before the strong earthquake. Mutual movements of small blocks follow the existence of high-amplitude, irregular spikes in low-frequency seismic noise waveforms. After consolidation, these irregular spikes disappear, which cause the decrease of ∆α and increase of entropy.

Trends of Seismic Noise Properties
From the point of view of studying the processes of preparation of large earthquakes, the behavior of the integral characteristics of the field of seismic vibrations is of particular interest. As such characteristics, let us consider the median values of the properties of seismic noise, calculated daily for all operational stations of the seismic network. The three graphs on the left in Figure 4 show the median values of ( , , ) En   , and green lines present running average values in a moving time window at the length of 57 days. The use of a moving average over a 57-day window is intended to facilitate visual perception of the daily median values of seismic noise properties. The window length of 57 days was chosen experimentally as a value that, on the one hand, smooths high-frequency pulsations of daily median seismic noise properties, and on the other hand, allows one to see the annual frequency of changes in these properties. The number The peculiarities of spatial distribution of the wavelet-based index γ, as can be seen from a comparison of Figure 3a1-a3,c1-c3, basically coincides with the main features of spatial distribution of the ∆α parameter. However, it should be noted that their calculation is based on completely different principles. Therefore, the γ parameter seems to provide additional information and their mutual consideration is important because estimations of these parameters are based on different approaches. It should be noticed that index γ in Figure 3c3  In addition, the Tohoku aftershocks region has lower values of ∆α in Figure 3a3 and higher values of entropy in Figure 3b3 for the same reason.

Trends of Seismic Noise Properties
From the point of view of studying the processes of preparation of large earthquakes, the behavior of the integral characteristics of the field of seismic vibrations is of particular interest. As such characteristics, let us consider the median values of the properties of seismic noise, calculated daily for all operational stations of the seismic network. The three graphs on the left in Figure 4 show the median values of (En, γ, ∆α), and green lines present running average values in a moving time window at the length of 57 days. The use of a moving average over a 57-day window is intended to facilitate visual perception of the daily median values of seismic noise properties. The window length of 57 days was chosen experimentally as a value that, on the one hand, smooths high-frequency pulsations of daily median seismic noise properties, and on the other hand, allows one to see the annual frequency of changes in these properties. The number 57 equals to double the length of the Moon month (28 days) plus 1. An additional unit is needed to make the window length odd.   To highlight possible predictive signs in the behavior of the median values of seismic noise properties, a deeper smoothing of high-amplitude random fluctuations should be performed. For this purpose, a kernel Gaussian smoothing [40,41] with a bandwidth of 182 days (half a year) will be applied. The result of this operation is shown in Figure 4 in the right column of the graphs. The red lines represent linear trends of smoothed values for the final observation interval, starting from 2012. From the behavior of the linear trends, it can be seen that, since 2012, there has been a systematic decrease in the singularity spectrum support width ∆α and the DJ index γ, as well as an increase in the entropy En of seismic noise. It should be noted that the final ∆α and γ values are less than any previous smoothed values, and the final entropy En value is superior to any previous value. Such behavior of the smoothed values is interpreted as an indicator of the growth of seismic hazard for the entire region of the Japanese Islands. The smoothed values of ∆α and En at right panel of Figure 4 have a rather explicit reaction on the Tohoku mega-earthquake on 11 March 2011: The singularity spectrum width ∆α increases whereas the entropy En decreases. As for the DJ index γ, its increasing is observed but its amplitude is rather weak. Thus, the statistic γ serves as some kind of additional characteristic of seismic noise to the "main" properties ∆α and En. Nevertheless, the Pearson correlation coefficient between smoothed γ and ∆α curves in the right panel of Figure 4 equals 0.77. Such a high value of correlation points out the fact that both these statistics reflect similar variations in seismic noise structure. Besides smoothed curves, it is interesting to compare simple mean values that are calculated for the sequence of three time intervals, which are discussed in Figure 3. These mean values are presented in the right panel of Figure 4 by blue horizontal lines. One can notice that mean values of γ and ∆α progressively decrease, whereas the mean value of entropy En increases. This fact is interpreted as the general seismic danger in Japan permanently increasing.

Maps of Probability Densities for Extreme Values
Let us consider values of the parameter as a function of 2D vectors of longitudes and ij ≡ U (t) (z ij ). For each daily "elementary map" with a discrete time index t let us find coordinates z n ) of the nodes where U attains a given number n m of extreme values with respect to all other nodes of the regular grid. If U = ∆α or U = γ, then the minimum values will be sought. Otherwise, for U = En, the maximum values will be sought. Further on, n m = 10 extreme values are used. The cloud of 2D vectors z (t) mn , which are regarded within some time interval t ∈ [t 0 , t 1 ], forms some random set. Let us estimate their 2D probability distribution function for each node z ij of the regular grid. For this purpose, the Parzen-Rosenblatt estimate with Gaussian kernel function [41] will be applied: Here, h is the radius of kernel averaging (smoothing bandwidth), t 0 , t 1 are intege indices that numerate daily "elementary" maps. Thus, (t 1 − t 0 + 1) is the number of daily maps within the considered time interval. The smoothing bandwidth h = 1 • was used. Figure 5 presents maps of probability density estimate (11) for time indices t corresponding to three time fragments similar to the maps presented at Figure 2. Kernel estimates (11) of the probability densities of extreme values of statistics of random fluctuations of geophysical fields in a moving time window were used in [6,42]. The distribution of probability densities of extreme values of seismic noise properties is shown only in the vicinity of the Japanese Islands in the union of circles with a radius of 250 km, built around each seismic station. The main purpose of constructing two-dimensional maps of the probability distribution of extreme values of the studied statistics is an attempt to more accurately localize those areas where their minima or maxima are most often realized, compared to using simple maps of property values, similar to those shown in Figure 3. The probability density distribution maps of extreme values are presented in Figure 5 for the same three time intervals as in Figure 2. From a comparison of Figures 3 and 5, it is noticeable that the maxima of the probability densities distinguish compact regions. Note that for the maps in Figure 5a2  The main purpose of constructing two-dimensional maps of the probability distribution of extreme values of the studied statistics is an attempt to more accurately localize those areas where their minima or maxima are most often realized, compared to using simple maps of property values, similar to those shown in Figure 3. The probability density distribution maps of extreme values are presented in Figure 5 for the same three time intervals as in Figure 2. From a comparison of Figures 3 and 5, it is noticeable that the maxima of the probability densities distinguish compact regions. Note that for the maps in  Figure 3. This area is interpreted as a possible epicenter for the next mega-earthquake. This hypothesis will be discussed in more details in the Discussion section.
It is noteworthy that the maps of the probability densities of extreme values of the properties of seismic noise in the left column of the maps in Figure 5 do not distinguish the epicenter of the 2003 earthquake. As a "justification", one can cite the consideration that the presented method for assessing "seismic hazard spots" is applicable only for mega-earthquakes such as the Tohoku event on 11 March 2011. The 2003 event off the coast of Hokkaido with a magnitude of 8.3 has an energy of almost 1.5 orders of magnitude less than the energy of the Tohoku event. Besides that, this could be the consequence of the small number of operable seismic stations during the initial time fragment of 1997-2001 (see Figure 2).

First Principal Component in Moving Time Window and Estimates of Coherence Spectrum
Furthermore, the method of principal components [43] for aggregating median daily time series of (En, ∆α, γ) into one scaled time series will be applied. This will be necessary for investigating the connection of seismic noise properties with the irregularity of Earth's rotation. A modification of the principal components method, which was proposed in [5], will be used.
Let us consider several time series of general dimensionality m P(t) = (P 1 (t), . . . , P m (t)) T , t = 0, 1, . . . In our case m = 3. It is necessary to estimate first principal component of P(t) in a moving time window of the length L samples. For this purpose, we will consider samples with time indices t under the condition s − L + 1 ≤ t ≤ s where s is the most right-hand end of time window. Correlation matrix Φ(s) of the size m × m is calculated according to formulae and define the scalar time series ψ(t) of the adaptive first principal component in a moving window according to the formula: Formulae (12)-(15) are applied independently within each time window. According to them, in the first time window, the time series ψ(t) consists of L values calculated according to (14)- (15). In all subsequent time windows, ψ(t) corresponds to the single most right-hand end. Thus, outside the first time window ψ(t) depends on past values of P(t) only.
For further analysis, it is necessary to estimate coherence spectra between two time series in a moving time window. A parametric model of vector autoregression that has a better frequency resolution than Fourier-based methods for estimating the spectra and cross-spectra [44] will be used. For a time series X(t) of dimensionality q, the AR-model is given by the formula: Here, t is the discrete time index, p is the order of autoregression, B k is the matrix of autoregression coefficients of the size q × q, and ε(t) is the residual signal covariance matrix P = M ε(t)ε T (t) of size q × q. The matrices B k and P are calculated by the Durbin-Levinson procedure [44]. The parametric estimate of the spectral matrix is defined by the formula: where E is the unit size q × q matrix. For dimension q = 2, the quadratic coherence spectrum is calculated according to the formula: Here, S 11 (ω) and S 22 (ω) are the diagonal elements of matrix (17) whereas S 12 (ω) is cross-spectrum. The coherence estimation was performed using a fifth-order vector autoregressive model with preliminary removal of linear trends and transition to increments.

Connection with Irregularity of Earth's Rotation
The irregular rotation of the Earth traditionally is explained by the influence of processes in the atmosphere [45]. At the same time, some researchers pointed out the connection between the irregular rotation of the Earth and seismicity [46,47]. The possible trigger mechanism of the influence of variations in the Earth's rotation on the seismic process was investigated in [48]. According to such an interpretation, a question arises about the impact of the processes in the atmosphere through the irregular rotation of the Earth on the seismic process. Estimates of the influence of the strong earthquake on the length of days are given in [49]. Figure 6a presents a graph of the length of day (LOD), which was taken from the source [50]. Figure 6b Figure 6d is a graph presenting maximum values with respect to the frequencies of the squared coherence. Figure 6e is a graph of the decimal logarithm of the seismic energy release (joules) in the vicinity of the Japan Islands: Information about seismicity was taken from the source [51].
In previous papers [11,52,53], maximum values of the coherence spectrum between LOD and daily median values of seismic noise properties were investigated in seeking the reasons for the existence of a break point in 2003 for trends and correlations of global seismic noise. It was shown that after 2003, trends of global seismic noise became typical for areas with increasing seismic hazard. Note that after the Sumatran mega-earthquake of 26 December 2004, M = 9, there was a sharp increase in the number of the strongest earthquakes around the world.

Discussion
The problem of a possible strong earthquake in the area of the Nankai Trough is traditional in Japanese seismology [54,55]. In [55], an estimate of the probability of 0.35-0.45 of the occurrence of an earthquake with a magnitude of more than 8.5 in this region during the time interval 2000-2010 is given. In [56], after the mega-earthquake of 11 March 2011, based on the analysis of GPS data, it was concluded that the probability of a strong earthquake to the south of the Tohoku aftershock region increased. Based on a retrospective analysis of seismic catalogs in [57], it was concluded that a magnitude 9 earthquake off the coast of Japan should not have come as a surprise. However, this event came as a surprise to traditional earthquake prediction methods. It led to an overestimation of the value of the maximum possible magnitude for seismic events in the Japan Trench, and in [58], an estimate is given up to a magnitude of 10.
A serious obstacle to the successful prediction of earthquakes is the presence of such a mechanism for the discharge of accumulated tectonic energy as "slow earthquakes" Thus, the structure of Figure 6 is the following. Figure 6c in the right-hand part of Figure 6 is a time-frequency map of the squared coherence function between two curves at the left-hand part presented by Figure 6a (LOD curve) and Figure 6b (first principal component of three daily median seismic noise properties). The map in Figure 6c illustrates a series of bursts of coherence within narrow frequency bands. This series of coherence bursts is presented as a one-dimensional graph in Figure 6d. The question of how much the bursts of coherence between LOD and the first principal component of the daily median properties of seismic noise can, on average, outpace the release of seismic energy is of interest. The estimate of this lag is presented in Figure 6e. To clarify this issue, the crosscorrelation function for time shifts of ±1200 days was calculated. The graph of correlation function is presented in Figure 6f.
According to the correlation function estimate, one can notice that its values for negative time shifts significantly exceed the values for positive shifts. This confirms the fact that the release of seismic energy is delayed relative to bursts of the coherence measure.
As for the lag time, it is estimated from the correlation maximum to be 432 days. This result also relates to the possibility of predicting strong earthquakes. Namely, if there is a burst of coherence between the LOD and the seismic noise properties, then this may be a harbinger of the strong release of seismic energy in about 1.5 years.
Although the value of correlation near 0.4 is rather small and it could not be an argument for a statistically significant linear connection between two random variables, I suppose that the main purpose of the correlation function estimate is establishing the fact that one variable shifted with respect to the other. The strong asymmetry of the correlation function in Figure 6f confirms the preceding of a coherence burst to a burst of seismic energy, which could be noticed visually by comparing the graphs in Figure 5d,e.

Discussion
The problem of a possible strong earthquake in the area of the Nankai Trough is traditional in Japanese seismology [54,55]. In [55], an estimate of the probability of 0.35-0.45 of the occurrence of an earthquake with a magnitude of more than 8.5 in this region during the time interval 2000-2010 is given. In [56], after the mega-earthquake of 11 March 2011, based on the analysis of GPS data, it was concluded that the probability of a strong earthquake to the south of the Tohoku aftershock region increased. Based on a retrospective analysis of seismic catalogs in [57], it was concluded that a magnitude 9 earthquake off the coast of Japan should not have come as a surprise. However, this event came as a surprise to traditional earthquake prediction methods. It led to an overestimation of the value of the maximum possible magnitude for seismic events in the Japan Trench, and in [58], an estimate is given up to a magnitude of 10.
A serious obstacle to the successful prediction of earthquakes is the presence of such a mechanism for the discharge of accumulated tectonic energy as "slow earthquakes" [59]. In fact, forecasting methods are only capable of identifying areas and time intervals where and when there is a noticeable accumulation of energy. The methods for studying low-frequency seismic noise discussed in the article can capture the accumulation of energy, which is associated with the consolidation of small blocks of the earth's crust into large blocks and, as a consequence of this consolidation, the simplification of the statistical structure of the noise due to the disappearance of high-amplitude bursts arising from the mutual motion of small blocks. The stored energy can be discharged both as a result of an ordinary earthquake, and during a "slow earthquake", the duration of which can be from several hours to several days. For the most part, slow earthquakes occur completely imperceptibly and are not recorded. Nevertheless, in terms of the efficiency of dissipation of accumulated tectonic energy, they are not inferior to ordinary "fast" earthquakes. Modern methods of analyzing geophysical data are still unable to give an estimate according to which of the two scenarios the accumulated energy will be released in. Due to this uncertainty in the practice of forecasting, there are many cases when the behavior of various characteristics of observations behaves similarly to previous cases before a seismic event; however, a new expected event does not occur, which means a possible release of energy according to a "quiet scenario".
Identification of areas in which minima or maxima of seismic noise statistics are most often realized, like the maps in Figure 5, provides a dynamic assessment of seismic hazard. The regions of maxima of two-dimensional distribution densities of extreme values of seismic noise statistics can be called "seismic hazard spots". Seismic hazard spots can arise and disappear without the manifestation in the form of an ordinary earthquake. At the same time, the detection of a stable "seismic hazard patch" like an area 30 • N ≤ Lat ≤ 34 • N, 136 • E ≤ Lon ≤ 140 • E means that there are persistent tectonic causes, and such areas should be given close attention.
As for the estimation of the event time, this part of the earthquake forecast is the most difficult. At present, one can only talk about the assessment of the trend, that is, whether the seismic hazard is increasing or decreasing. In this regard, the linear trends presented in the right column of the graphs in Figure 3 indicate an increase in the seismic hazard of the next mega-earthquake. In addition, the graph of the cross-correlation function presented in Figure 6f can give an estimate of the time of the next large earthquake from the occurrence of a burst of the coherence function between the median properties of seismic noise and the time series of the length of day reflecting the irregular rotation of the Earth.

Conclusions
A method for analyzing continuous records of low-frequency seismic noise on a network of broadband stations in a seismically active region is presented. The method is based on the calculation of multifractal singularity spectrum support width and waveletbased entropy and Donoho-Johnstone index in successive daily time intervals. Methods have been developed for assessing the spatial distribution of the values of these seismic noise properties and the identification of "seismic hazard spots" as areas of increased values of two-dimensional probability densities of extreme values of seismic noise statistics. The method is applied to the analysis of data from the F-net seismic network on the Japanese Islands for the time interval from the beginning of 1997 to the end of March 2021. The question of a possible trigger effect of the irregularity of Earth's rotation on the seismic energy release is considered. An estimate of the lag between strong bursts of seismic energy release and bursts of coherence between the time series of the length of the day and the median values of seismic noise properties in Japan is obtained. The value of this lag is close to 1.5 years, which could be used for announcing a high seismic danger time interval after the peak of coherence is achieved.
Funding: This research received no external funding.