Extracting Correlations in Earthquake Time Series Using Visibility Graph Analysis

Recent observation studies have revealed that earthquakes are classified into several different categories. Each category might be characterized by the unique statistical feature in the time series, but the present understanding is still limited due to their non-linear and non-stationary nature. Here we utilize complex network theory to shed new light on the statistical properties of earthquake time series. We investigate two kinds of time series, which are magnitude and inter-event time (IET), for three different categories of earthquakes: regular earthquakes, earthquake swarms, and tectonic tremors. Following the criterion of visibility graph, earthquake time series are mapped into a complex network by considering each seismic event as a node and determining the links. As opposed to the current common belief, it is found that the magnitude time series are not statistically equivalent to random time series. The IET series exhibit correlations similar to fractional Brownian motion for all the categories of earthquakes. Furthermore, we show that the time series of three different categories of earthquakes can be distinguished by the topology of the associated visibility graph. Analysis on the assortativity coefficient also reveals that the swarms are more intermittent than the tremors.


Network-Theoretical Time Series Analysis
Inspired by the exceptional success of the network theory in recent years [1][2][3][4][5][6][7], the analysis of time series from the perspective of complex network has received considerable attention due to the standing requirement of understanding the dynamical processes behind time series data [8][9][10][11][12]. Often a real-world time series arises from non-linear processes and their precise identification is important for modeling purposes. Recently, a merging trend has been observed coupling ideas both from the field of non-linear time series analysis and complex network theory [13]. If a time series is mapped into a complex network, one may expect that such a network reflects some inherent properties of the original time series. Thus, one can utilize the recent graph-theoretical tools to extract novel properties hidden in the time series.
Among several other methods [11,12], the visibility graph [10] has become popular due to its simplicity and wide range of applicability. This method has demonstrated its potential in extracting several characteristic features of the time series such as the periodicity, fractality, chaoticity, non-linearity, and more [10,14,15]. A merit of the visibility graph method is its ability to capture non-trivial correlations in non-stationary time series without introducing elaborate algorithms such as detrending. For instance, it has been shown that the visibility graph corresponding to the time series generated from a fractional Brownian motion (fBm) is scale-free. Moreover, the exponent γ for the degree distribution corresponds to the Hurst exponent (H) of the fBm as [16]: Since the fBm generates f −β power spectrum with β = 1 + 2H, the exponent γ of the visibility graph should correspond to β as The network-theoretical method enables us to estimate H and β more easily than other standard methods such as calculating power spectrum [16]. Therefore, it has been applied to extract the fBm-like nature of time series in several contexts such as finance [17], health science [18,19], image processing [20], and geophysics [21,22].
In this paper, we study the nature of correlation in earthquake time series by means of visibility graph. In particular, we focus on the two important quantities: the magnitude and the inter-event time (IET) between two consecutive earthquakes.

Characteristics of the Seismic Sequences: Three Categories of Earthquakes
Thanks to the continuous progress in observation technologies, various kinds of earthquakes have been known to date. Aiming at the statistical characterization of earthquakes belonging to different categories, here we choose to analyze three wellestablished categories: regular earthquakes, earthquake swarms, and tectonic tremors. The fundamental difference among these three categories lies in their generation mechanisms and the time scale of energy release.
A time series of regular earthquakes includes mainshockaftershock sequences and the background activity. While the latter is a Poissonian process, the former is generally clustered in space and time. Aftershocks are triggered usually by the static stress change associated with the mainshock, as well as some other post-seismic relaxation processes such as afterslip or fluid flow. Major fraction of the total energy is released almost instantaneously at the time of the mainshock and slowly decreases in time. It is observed that the magnitude-frequency distribution P(M) obeys an exponential distribution, namely, the Gutenberg-Richter (GR) law [23]: P(M) ∝ 10 −bM , with b taking a value around 1 in the active fault zones [24]. On the other hand, the temporal decay of the frequency of aftershocks is described by the Omori-Utsu law [25,26].
The same phenomenology is not observed for the other two categories of earthquakes. In contrast to mainshock-aftershock sequence, a seismic swarm is defined as a cluster of earthquakes with similar magnitudes, which usually occur in a volcanic or geothermal tectonic setting. The intrusion of fluids can reduce the resistance of faults and redistribute the stress in such a manner that the energy is released gradually and almost equally among the largest shocks [27]. The Omori-Utsu law does not generally hold for swarms.
Tectonic tremors represent weak and repetitive seismic signals emitted from a plate boundary in a subduction zone. To the current belief, fluids generated by slab dehydration may be a cause of tremors [28]. Similar to swarm earthquakes, the tectonic tremor activity is characterized by hypocentre migration but on a different spatial and temporal scale: tremors migrate up to several hundreds kilometers, whereas swarms are more local. The statistical laws are largely unknown for remors.

Outline of the Paper
Based on the analysis of the visibility graph, we argue against the current popular belief that earthquake magnitude time series are indistinguishable from random time series. The same method is applied for the IET time series, showing fBm-like correlation clearly. We also show that the time series of three different types of earthquakes can be distinguished in the topology of the associated visibility graph.
The paper is organized as follows. We start by describing the visibility graph algorithm and the characteristics of the three categories of earthquakes including the specifications of the studied seismogenic zones in section 2. The existence of memory in the time series of magnitudes and IETs have been investigated in sections 3 and 4, respectively. We discuss the topology of the visibility graph for both magnitude and inter-event time series in section 5. Finally, we summarize in sections 6 and 7.

Construction of Visibility Graph From Seismic Catalog
Given the time sequence of the occurrence of seismic events, the visibility graph is constructed by considering each event as a node and linking the nodes based on mutual visibility of the corresponding data heights. The data recorded at time t k is represented as the height h k of the k-th node. Specifically, any arbitrary pair of data values (t i , h i ) and (t j , h j ) (t i < t j ) are visible to each other if the straight line joining the two data points does not intersect any intermediate data heights, as illustrated in Figure 1.
If there exists visibility, the slope s ij of the line between the nodes i and j must be the maximum of the slopes s ik for all i < k < j. Therefore, a link is placed between two nodes i and j in the visibility graph if and only if for all t i < t k < t j the following criterion is satisfied: Clearly, every node is visible at least from its left and right nearest neighbors and thus one obtains a completely connected network. The "divide & conquer" algorithm [29] has been used to efficiently transform a time series into its corresponding visibility graph. This algorithm takes advantage of the fact that the node with the maximum height divides the time series into two segments in the sense that the nodes situated at one side of the maximum are not visible from the another side. Therefore, it is not required to check the visibility between the two sides of each separated segments. In each step, the visibility of the node with the maximum height to the other nodes at its right and left sides is determined. Each new segment is then treated independently and the same procedure is repeated until every segment contains one single node. The CPU time taken by the algorithm scales with the size N of a time series as N log N.

Description of the Seismic Catalog
In a seismic catalog, an event is described by the location of the hypocenter, the time of occurrence, and the magnitude (M). We select several representative regions from Japan and California since these two areas are well-known for intense seismic activity and dense monitoring networks. The catalog data we analyze here are provided by the Japanese Meteorological Agency 1 , the Hot Spring Research Institute [30], the Southern California Earthquake Center [31], the World Tremor Database [32], and Slow Earthquake Database [33], respectively. 1 ftp://ftp.eri.u-tokyo.ac.jp/pub/data/jma/mirror/JMA_HYP/ FIGURE 1 | Visibility graph representation of a synthetic time series with 20 data values drawn randomly from an exponential distribution, where t i = i is the time t corresponding to the i-th data. Each vertical bar representing the height variable h is considered as a node and if the top of one bar is visible from the top of the another then a link is placed between the corresponding pair of nodes.
A selected region is described by the minimum and the maximum of the latitude (θ ) and longitude (φ) coordinates, i.e., the values of (θ min , φ min ) and (θ max , φ max ). We consider only the crustal events within the depth of 50 km. For the regular and the swarm earthquakes, we also indicate the magnitude of completeness M c i.e., the lowest magnitude above which the GR law holds. Above this completeness magnitude, missing events in a catalog should be rare and therefore, effects of missing events should be minimized. We determined these values using the Zmap software tool [34]. For tremors, we consider all detected events recorded in the two previously mentioned database [35,36]. The total number of events in a catalog is denoted by N t . The detailed specifications of these catalogs data are given in Table 1.

Remarks on Regional Specifics
For time series of regular earthquakes, we analyzed three active seismic regions located in different tectonic settings: subduction, compression, and active faulting. The region named Tohoku corresponds to an offshore area of the Japan Trench subduction zone where the 2011 earthquake of moment magnitude M w 9.0 and its aftershocks were recorded. Time series before and after the M w 9.0 event are referred here as Tohoku 1 and Tohoku 2 , respectively. The Southern California region is located in a complex compressional tectonic setting dominated by the southern part of the San Andreas Fault system, but also includes earthquakes generated by the slow uplifting of the Sierra Nevada Mountain range, as well as volcanic and geothermal related activity. The Kumamoto region mostly includes the recent seismic activity generated by the 2016 M w 7.0 Kumamoto earthquake around the active Futagawa-Hinagu fault and the surrounding active volcanic region of Aso-Yufuin-Beppu. Thus, most earthquakes in the Kumamoto catalog are aftershocks. In the Hakone volcanic region, significant swarm activity was detected since 2001 [37]. Although many different swarm episodes were recorded, they don't exhibit any specific temporal pattern. An increase in the seismicity level was observed in 2015 due to a volcanic eruption [38]. The Izu volcanic region is characterized by magma-intrusion episodes which generate frequent swarm activity [39]. Concerning the tremor activity, we selected two areas where the largest number of detected events is available, such as Cascadia in North America and Shikoku around the Nankai Trough in Japan.

Stretched Exponential Nature of Degree Distribution
To investigate whether the magnitude of earthquakes has any correlations, we study the degree distribution of the visibility graph constructed from the magnitude time series. First we check if the degree distribution is power law. Typically, a power law distribution is characterized by a long tail that develops with the network size N in such a manner that the average maximum nodal degree k max (N) grows as k max (N) ∝ N α . This signifies the existence of power-law degree distribution for the infinitely large network, N → ∞. In order to do this analysis, the original time series is divided into several segments such that each segment contains exactly N number of events.
We start with our results for regular earthquakes in the Tohoku region. Since the period of Tohoku 2 is exceptionally active after the occurrence of the magnitude 9.0 earthquake, we have analyzed the data for Tohoku 1 and Tohoku 2 separately. In Figures 2A,B, the degree distribution of the visibility graph is shown on a double logarithmic scale for four values of N starting from 2 10 to 2 13 , at each step N being increased by a factor of 2. For all the four values of N in both the cases (Tohoku 1 and Tohoku 2 ), the curves have certain amount of curvature and the tails of the degree distributions do not elongate significantly as N increases. To see this dependence more clearly, we have plotted the average maximum nodal degree k max (N) against N on a semilog scale in Figure 2C. Clearly, this implies that k max (N) ∼ ln N, demonstrating that the degree distribution is not a power law: namely, the absence of fBm-like structure in the magnitude time series.
Specifically, the degree distribution appears to follow a stretched exponential function: In Figure 3A, we have plotted the degree distribution p(k) of the visibility graph on a log-log scale for the whole time series of Tohoku 1 and Tohoku 2 containing 55824 and 91197 events, respectively. The logarithmically binned data for both the series fits quite well with the above functional form in the range of k between 6 to approximately 100. This is shown more explicitly in Figure 3B, where p(k) is replotted against k τ on a semilog scale. The curves are straight in the intermediate region, indicating that the distribution follows an exponentially decaying function of k τ . This behavior is also evident from the cumulative degree distribution shown in Figure 3B (inset).
To confirm the ubiquity of the stretched exponential nature of degree distribution, we analyze the other six earthquake catalogs. Figure 4 shows degree distributions presented similarly to those in Figure 3B for Southern California (regular), Hakone (swarms), and Shikoku (tremors). Apparently, these degree distributions are fitted with the stretched-exponential function irrespective of the region or the earthquake type. The cumulative degree distributions are also shown in Figure 4 (inset).
Furthermore, two important points should be remarked regarding the robustness of the above result. First, the stretched exponential nature does not significantly change even when the cutoff magnitude M c is set to be lower or slightly higher than the completeness magnitude: Namely, the result is rather insensitive to some undetected smaller events. This may be because the tail of the degree distribution is controlled by events of larger magnitude, which generally have higher visibility. Second, we confirm that the shape of degree distribution is unaltered even if the time series is with respect to the event index instead of the real occurrence time. Namely, the degree distribution remains stretched exponential even if the event time t i is replaced by an integer i in the visibility criterion, Equation (3).

Degree Distribution for Shuffled Data
Hereafter we refine the analysis and argue if there are any other correlations in the magnitude time series. To this end, we first analyze the visibility graph produced from the shuffled time series. Namely, by randomly choosing a pair of events, their respective magnitudes are swapped. This process is repeated by N t times (the number of events in the catalog), leading to one shuffled sequence. This procedure preserves the probability density function of magnitude but destroys any potential correlations between them. Then, for a shuffled sequence, the visibility graph is constructed and the degree distribution is calculated. This process is repeated for many times and the degree distributions are averaged over these shuffled sequences. The averaged degree distribution is shown in Figure 5 (main panel). This is again fitted with the stretched exponential distribution. The same is true for the degree distribution of each shuffled sequences, and the curves are not distinguishable from the original time series (inset of Figure 5). This again validates the absence of fBm-like correlations in the original sequence.  with k being the degree of the nodes, for the shuffled series corresponding to Tohoku 1 (black) and Tohoku 2 (red). The plot is based on 10 6 independent shuffled series. Inset: Log-log plot of the cumulative degree distribution P(k) for six individual shuffled series of Tohoku 1 (different colors are used to represent different shuffled series) along with the original one (black).

Kolmogorov-Smirnov Test
More importantly, however, the above analysis does not mean that there are no correlations in earthquake magnitude, since the averaging process may mask some subtle short-range irregular correlations. To scrutinize the statistical difference in the visibility graph structure of the original and the shuffled time series, we perform the Kolmogorov-Smirnov (KS) test. Here the null hypothesis is that two empirical degree distributions originate from the same function for the original time series and its shuffled sequence. We adopt the 0.05 significance level and reject this null hypothesis if the p-value is smaller than 0.05. In this formulation, rejecting the null hypothesis means that the degree distributions are different for two visibility graphs produced from the original time series and its shuffled one.  Specifically, the KS test statistic is computed as a distance between two empirical degree distributions produced from the original time series and its shuffled one. Then the p-value is calculated from the distance. This procedure is repeated for many shuffled sequences to yield the distribution of the p-value. They are shown in Figure 6A for Tohoku 1 (regular earthquakes) and Figure 6C for Cascadia (tremors). Apparently, the null hypothesis is rejected for both the cases. Namely, the degree distributions are not the same for the original time series and the shuffled surrogates. We also find that the null hypothesis is rejected for all the other catalogs shown in Table 1. This implies that the visibility structure in the original time series is somewhat altered if shuffled. In other words, the original time series can be discriminated among many other shuffled data.
To support the above statement from another aspect, we again perform the KS test by comparing a specific shuffled time series with many other shuffled ones. The distribution functions for the p-value are shown in Figures 6B,D. In this case, the null hypothesis is not rejected at the 0.05 significance level. Namely, shuffled time series are indistinguishable in terms of the degree distribution of their visibility graphs. This makes a quite contrast to the original time series, which is distinguishable from shuffled ones.

Analyses on Three Other Surrogates
In addition to shuffled time series investigated above, we inspect three other surrogate data. The first and the second ones are the random time series, in which the height values {h i } are drawn randomly and independently from an exponential distribution p(h) ∼ e −λh between [2,9]. For the first surrogate data, the time is set to be the event index: i.e., t i = i for the i-th event. Note that λ is proportional to the b-value in the GR law as λ = 2.303b. In Figure 7, the degree distribution p(k) are shown for several values of λ. Each curve is seen to follow the stretched exponential form. Similar to the original earthquake data, k max (N) grows logarithmically with N (not shown). This makes a contrast to the exponential degree distribution observed for the uniformly distributed heights [10]. The second surrogate data is the Poisson model, where events occur according to the Poisson process, and the height values are again drawn randomly from the GR law. We confirm that this surrogate data also produces the stretched exponential degree distribution. However, in the KS test that compares the surrogate data and the original magnitude series, the null hypothesis is rejected. Namely, they don't yield the same degree distribution.
The third surrogate data we wish to inspect is a time series with a short memory. Here the time series is generated by simulating a Brownian particle in one dimension subjected to a linear potential: U(x) = c|x|. Starting from x = 0 at time t = 0, the position of the particle is updated in steps of dt = 10 −6 according to the following Langevin equation: where ξ is a Gaussian white noise with zero mean and unit variance. The height distribution for x(t) follows the Boltzmann-Gibbs distribution at equilibrium. Since the potential is linear, the distribution function is exponential, as confirmed in Figure 8A.
We construct the visibility graph using the time series of x(t) and compute the degree distribution. As shown in Figure 8B with four different system sizes N, the degree distribution is observed to follow a power law. Additionally, we confirm that k max (N) grows as a power-law with N: i.e., k max (N) ∼ N 0.486(5) (not shown). This signifies that a systematic single step memory in the time series leads to a scale-free network. All the findings above lead us to conclude that the time series of earthquake magnitude are not statistically identical to uncorrelated time series, although no apparent systematic memories exist, either long-ranged (fBm-like) or short-ranged.

Power-Law Nature of Degree Distribution for Tohoku Data
To characterize the temporal correlations between seismic events and to understand whether they are dependent on specific details of the seismic activity, we focus on studying the inter-event time (IET) series of earthquakes. Here the IET series is obtained from an earthquake catalog by calculating time intervals between two consecutive events and labeling them with the event index i. Namely, the IET series is represented as (i, h i ), where h i = t i+1 − t i , and t i is the real occurrence time of the i-th event in a catalog. Here the threshold is set as the completeness magnitude M c (listed in Table 1). Figure 9A shows the cumulative degree distribution P(k) for the IET series of Tohoku 1 and Tohoku 2 . This is the probability of finding a node with degree at least k in the visibility graph. For both the cases, the degree distribution is found to be heavytailed distribution and the tail extending more than one decade can be described by an approximate power law. We estimate the exponent: γ = 2.34(5) for the Tohoku 1 and γ = 2.60(8) for the Tohoku 2 . The average maximum nodal degree also varies as a power law: k max (N) ∼ k α , where α = 0.77(3) and 0.53(4) for the Tohoku 1 and Tohoku 2 , respectively (not shown here). This behavior supports the power law nature of the degree distribution. Thus, the visibility graphs constructed from the IET series exhibit typical signatures of a scale-free network, indicating the existence of fBm-like correlations in the time series.
To validate the presence of correlation in a contrasting manner, we analyze the shuffled sequences of the IET data and find that the degree distribution p(k) is fitted with the stretched exponential function given in Equation (4). In Figure 9B, the degree distributions p(k) are plotted with k τ for the shuffled IET series of Tohoku 1 and Tohoku 2 . The straight line here confirms the stretched exponential form of the degree distribution. In addition, we find that k max (N) ∼ ln N (not shown). Evidently, the shuffled data produces the properties of a random time series and therefore provides evidence on the existence of correlation in the original time series.

Power-Law Nature of Degree Distribution: Other Regions
The same analyses are carried out for regular earthquakes in different regions, as well as for swarms and tremors. The results are shown in Figure 10. In Figures 10A-C, the cumulative degree distribution is plotted for regular earthquakes, swarms, and tremors. For every case, a heavy-tailed distribution has been observed. While for regular earthquakes and tremors a power law regime extending more than one decade is quite apparent, the data for swarms shows more complex behavior. However, an approximate power law variation can fit the data in the intermediate region. For each case, the data points in the most linear regime (estimated by eyes) starting from a moderate value of k to a value at the tail part upto which they do not fall-off due to the limitations by finite size are fitted to the best straight line. From the slope of the straight line, we estimate the powerlaw exponent γ as 1.73(8), 2.64(5), 1.81(9), 1.79(9), 2.51 (5), and 2.13 (5) for Kumamoto (regular), Southern California (regular), Hakone (swarm), Izu (swarm), Cascadia (tremor), and Shikoku (tremor), respectively. In addition, the power law dependence of the average largest degree k max (N) with N has been observed for every set of data (not shown), supporting the power law nature of the degree distribution.
The tail part of the degree distribution is characterized by the exponent γ , which seems to depend on the seismic activity of the specific region: (i) Earthquake swarms (Izu and Hakone) have a common value, γ ≃ 1.8. (ii) Regular earthquakes may also have a common value, γ ≃ 2.6 (Tohoku 2 and Southern California), while it is somewhat smaller (2.3) before the Tohoku M w 9.0 earthquake (Tohoku 1 ). (iii) Kumamoto is exceptional with γ ≃ 1.7. This value is rather close to swarms, although the data mainly consist of aftershocks of 2016 Kumamoto earthquake. There may be two reasons for this discrepancy. First, the data is not a usual mainshock-aftershocks sequence, but rather a foreshocksmainshock-aftershocks sequence. Alternatively, we may interpret it as two mainshocks (M w 6.2 and 7.0) that occurred within only 30 h. In any case, it is rather anomalous seismic activity. The second potential reason is an active volcano (Mt. Aso) located in the proximity of the main fault. The M w 7.0 mainshock triggered many earthquakes in the volcanic area, including an M w 5.9 event and its own aftershocks. Thus, the overall seismic activity is influenced by the nearby volcanic field and this may explain the resemblance to swarms.
If we suppose the relation between the fractional Brownian motion and the power-law degree distribution, i.e., Equation (2), the exponent for the power spectrum β can be determined. For example, swarms have β ≃ 2.2 and H ≃ 0.6. They are close to those for standard Brownian motion (β = 2 and H = 0.5) but yet slightly larger, corresponding to superdiffusion. Regular earthquakes (γ ≃ 2.6) have β = 1.4 and H = 0.2, corresponding to subdiffusion. Extraction of these exponents from actual seismic data is difficult using other standard methods such as autocorrelation functions due to the strong non-stationary nature of the seismic record. In this sense, these exponents might not be considered as that for fBM itself, but should represent some counterpart in seismic activities.
Shuffled time series again yield visibility graphs with their degree distributions of stretched-exponential form, resembling the properties of a random time series. Therefore, we confirm that the original IET series possess fBm-like correlations irrespective of the earthquake types: regular, swarms, and tremors. This result does not contradict the previous studies on regular earthquakes obtained using some different methods [40,41]. Here we have confirmed the correlation using complex network based approach, and more importantly, found correlations in tremors and swarms.
Lastly, we wish to add a remark on catalogs on tremors. Since a single event is not as distinct as regular earthquakes, there may be some errors in the IET of tremors. To check the effect of such errors in IET, we add a certain amount of noise to the IET data of tremors and construct the visibility graph from these noisy data. We find that the degree distribution is indeed robust to the noise: it retains the power-law nature against the small noise in IET.

DETAILED STRUCTURE OF VISIBILITY GRAPH
The detailed characterization of the topology of the network has served to identify several non-trivial features exhibited by diverse types of real-world systems including the basic principles that played role in the network formation [1,3,4]. In order to extract more properties hidden in the seismic records, the following graph-theoretical quantities have been analyzed after obtaining the visibility graph using "divide & conquer" algorithm.
Since our visibility graph is connected and undirected, there always exists at least one path between any arbitrary pair of nodes i and j through the links of intermediate nodes. The path with the minimal links traversed is called the shortest path length d ij , and the average shortest path length is defined as, In Figures 11A,B, we show the variation of l(N) with N on a semilog scale for both the magnitude and IET series, respectively. The best fit of the data by a straight line indicates its logarithmic scaling and hence, the network is small-world. Although the data for IET series of Shikoku has some curvature, the linear behavior is quite apparent for large values of N. For IET series of swarms, l(N) grows more slower than ln N. Another important quantity associated with the network is the clustering coefficient which measures the three point correlation among the neighbors. Specifically, the clustering coefficient C i of node i measures the probability that the two neighbors of i are connected. If there exists E i links among the k i neighbors of node i then, C i = 2E i /k i (k i − 1). In the case of k i < 2, C i = 0. The global clustering coefficient is expressed as, By varying N from 2 9 to 2 16 we have observed that C is almost independent of N (values differ only at 4-th decimal place) for both magnitude and IET time series of different types of earthquakes. Further, the clustering coefficient C(k) for the nodes with degree k has been found to decay as C(k) ∼ k −ν with ν ≈ 1, as shown in Figure 12. This is the universal feature of a hierarchical network observed in many real-world networks [42]. The clustering coefficient C assumes its highest value for the IET series of tremors. We have also calculated the Pearson correlation coefficient r to investigate whether a high degree node tends to be linked with a high degree node (assortative mixing, r > 0) or a low degree node (disassortative mixing, r < 0). We have calculated r using the following formula [43], where, k 1i and k 2i are the degrees of nodes at the ends of link i with i = 1, 2, · · · , L. We found that for all earthquake types, the magnitude series shows assortative nature (last column of Table 2). In contrast, in case of IET series we obtain a value of r ≈ 0 for the regular earthquakes and for swarms and tremors r < 0 (last column of Table 3). Moreover, the graph associated with the IET series of swarms has been found to be more disassortative than that of tremors. This means that for swarms the high degree nodes show more preference toward linking with the low degree nodes. This indicates that the smaller heights are abundant in both the time series, however, there are a few very large heights (i.e., long quiescence periods) in the swarms series which are even larger than the largest height in the tremor series. Therefore, swarms are more intermittent than tremors. For a detailed comparison of the characteristic differences among the three different types of earthquakes, the above quantities have been calculated for a fixed value of N = 2 12 and the obtained values are listed in Tables 2, 3 for the magnitude and the IET series, respectively. Clearly, they can be distinguished by the values of different graph-theoretical quantities obtained from their individual IET series.

DISCUSSION
Finding and characterizing any correlations in the time series of earthquake magnitude is a subject of great importance as it may be useful in forecasting major earthquakes. However, to date, existence of correlations is somewhat controversial and has not been settled [40,[44][45][46]. For instance, it was reported that regular earthquakes occurring close in space and time are correlated in their magnitudes [44]. A counterargument was given in [45] that these were pseudo correlations due to the magnitude incompleteness and the modified Omori law. To shed new light to this long-standing problem, we have made use of the complex network theory and analyzed the visibility graph to extract correlations in magnitude time series.
The previous studies [47,48] in this context involve regular earthquakes only. Here we extend the analysis to two other types of earthquakes [49] to consider this problem in a more general perspective. By using the method of visibility graph, we have analyzed seismic time series in seven seismogenic zones including the regular earthquakes in Southern California in common with [48] but for more extended time period. The degree distribution appears to be fitted with a stretched exponential function for all the types of earthquakes analyzed here.
Visibility graphs are also constructed from shuffled catalogs ( Figure 5) or synthetic data drawn randomly from the GR law (Figure 7). On average, the degree distribution appears to be fitted with the stretched exponential function. However, the Kolmogorov-Smirnov test rejects the null hypothesis that these degree distributions are identical. Namely, the degree distributions for these surrogate data are indeed distinguishable from that of the original data. This means that the original series have some special characters that are lost in their surrogates: shuffled or synthetic catalogs.
Since the criterion for the visibility graph involves both magnitude and IET, one might argue that the difference in the degree distribution detected by the KS test is a mere by-product resulting from the correlation in IET. To exclude this possibility, we also perform the KS test by constructing the visibility graph using the event index i instead of the occurrence time t i . We find that the null hypothesis is again rejected. Namely, the magnitude series (i, M i ) leads to slightly different degree distributions if they are shuffled, although the difference is detectable only by the KS test. Thus, the memory should exist in magnitudes alone.
The memoryless nature of earthquake magnitudes is a basic assumption in the epidemic-type aftershock sequences (ETAS) model, which is the most successful statistical model for earthquake time series [50]. The results given here implies that the memoryless assumption in earthquake magnitude is rather approximate. Thus, if one wishes to improve statistical models for earthquake occurrence, the correlation in magnitude should be taken more seriously. To this end, the correlation found here should be defined and quantified more clearly.
The degree distributions of stretched exponential form appear to contradict some previous studies [47,48], in which the power law tails are concluded for the magnitudes of regular earthquakes. In view of Equation (1), this may imply a fBm-like correlation in the magnitude time series. Interestingly, however, they also analyzed randomly shuffled sequences of magnitudes and did not find any significant difference in the degree distributions. This rather contradicts the existence of a fBm-like correlation. Additionally, the degree distribution obtained in [47] spans approximately one decade only, and the tails are noisy. Thus,  one needs to be careful to draw a conclusion based on these data alone. In [48], the tails of the degree distributions are less noisy, but they appear to fall off from a power law at their tails. Thus, their degree distributions might be fitted with a stretched exponential function. However, the degree distribution produced from Mexican catalog appears to develop a tail that is still different from stretched exponential. We noticed that the magnitude data in the Mexican catalog do not always obey the GR law, and this may be the reason for the deviation from the stretched exponential function. However, the Mexican data require more careful and dedicated analyses to draw any decisive conclusions on specific type of magnitude correlation. We apply the visibility-graph analysis for the characterization of the inter-event times (IET) between consecutive earthquakes. Contrary to the magnitude time series, we find an evidence of fBm-like correlations between the inter-event times. The network associated with the IET series has a scale-free nature with the exponents γ , which depends on the essential characteristics of seismic activity. In the context of the f −β noise, the exponent γ is directly related to β. These exponents may work as a generalized and unified quantification of the intermittent nature of seismic time series. For instance, we find that the IET series for swarms are similar to superdiffusive Brownian motion, whereas those for regular earthquakes correspond subdiffusion. However, the interpretation of superdiffusive or subdiffusive nature in the IET series is yet unclear from the mechanical point of view, and should be pursued in the subsequent studies.
We have also analyzed the whole set of data using the horizontal visibility graph algorithm. For both magnitude and IET series, however, the degree distribution results in an exponential distribution and no appreciable change has been observed between these different data sets, making it harder to draw any conclusive remarks on the distinction of different time series.

CONCLUSION
In conclusion, we have investigated the correlations in the time series of magnitudes and of inter-event times (IETs) for three different categories of earthquakes in seven seismogenic zones in the world. By applying the methods of visibility graph, we show that the IET series possess correlations similar to fractional Brownian motion, and that the three categories of earthquakes have different exponents. While such an apparent correlation is absent in the magnitude series, the Kolmogorov-Smirnov test on the degree distribution reveals that the earthquake magnitudes are not statistically equivalent to an uncorrelated (random or shuffled) time series. This challenges a current popular belief that magnitude time series are random. Since current major statistical models for earthquake rate are based on this belief, these results provide us with useful constraints in developing better statistical models.
Different temporal behaviors of three categories of earthquakes are also reflected in various graph-theoretical quantities. As found from the analysis of the assortativity coefficient, the swarms are more intermittent than tremors. More graph-theoretical techniques including horizontal visibility graph [51], multiplex visibility graph [52], and recurrence networks [11] would give new criteria for categorizing or unifying different seismic activities. A novel approach for forecasting time series based on visibility graph [53] might find potential application for earthquakes. Our study therefore shows with affirmation that the visibility graph algorithm has the potentiality to capture the non-trivial complexity inherent in a time series which is non-linear and non-stationary in nature.
One can also consider more elaborated methods for the graph construction [54]. For instance, the visibility graph constructed here is undirected and unweighted. Time directionality and weighted links based on the inter-event distances would be interesting subjects. Additionally, since the spatial information of the seismic events has been disregarded here, the extension of the visibility graph method to space-time may be a promising attempt.
Together with the present results, such graph-theoretical approaches would bring benefits to statistical modeling of various types of seismic activities that cannot be reproduced by the well-established ETAS model for regular earthquakes.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.