Correlations and invariance of seismicity under renormalization-group transformations

The effect of transformations analogous to those of the real-space renormalization group are analyzed for the temporal occurrence of earthquakes. The distribution of recurrence times turns out to be invariant under such transformations, for which the role of the correlations between the magnitudes and the recurrence times are fundamental. A general form for the distribution is derived imposing only the self-similarity of the process, which also yields a scaling relation between the Gutenberg-Richter b-value, the exponent characterizing the correlations, and the recurrence-time exponent. This approach puts the study of the structure of seismicity in the context of critical phenomena.

The study of statistical seismology has a long history, exemplified by the Omori law of aftershocks and the Gutenberg-Richter relation for the number of earthquakes above a given magnitude, and more recently, by the fractal properties of earthquake spatial occurrence [1,2,3,4]. Less attention however has been paid to the timing of individual earthquakes, for which a unifying picture was missing until the work of Bak et al. [5,6,7]. The main relevance of that work was the seminal introduction of scaling concepts in earthquake statistics, which provide a powerful tool to unify descriptions and to derive relations between different quantities (in this case, as we will see, by scaling we do not mean just power-law relations, as it is sometimes assumed). Later, Bak el al.'s procedure was modified to study the distribution of times between consecutive events in a single spatial region [8].
Let us consider the seismicity of an arbitrary spatial region. Given a lower bound M c for the magnitude, and intentionally disregarding the spatial degrees of freedom, a marked point process in time of the form (t 0 , M 0 ), (t 1 , M 1 ), . . . is obtained, where t i denotes the time of occurrence of event i, with a magnitude M i ≥ M c . The recurrence times are defined as the time intervals between nearest-neighbor (i.e., consecutive) events, τ i ≡ t i − t i−1 . In the case of stationary seismicity (characterized by a long-term linear relation between the accumulated number of earthquakes and time), for spatial regions of linear size ranging from 20 km to the whole world, and for magnitude bounds from 1.5 to 7.5, the probability densities D(τ) of the recurrence time were found to verify a universal scaling law [8], where f is a universal scaling function and the scaling factor R is the rate of seismic occurrence, defined as the mean number per unit time of events with M ≥ M c in the region, and given by the Gutenberg-Richter law, R(M c ) = R 0 10 −bM c , with the b−value usually close to 1. As no separation of mainshocks and aftershocks is performed, the recurrence-time distribution consists of a mixture of different aftershock sequences and more or less independent events; therefore, it is very surprising that distinct regions and earthquakes of disparate sizes present such a extreme degree of regularity.
But even more surprising, the scaling relation for the recurrence-time distribution reveals that seismicity is in a highly orchestrated state, in which the removal of events (when the lower bound M c is raised) does not affect the properties of seismic occurrence, as the distribution keeps the same shape (with only a different mean) independently of M c . In general, when some events are removed from a point process, the properties of the process do change; therefore, the distribution of recurrence times constitutes a very special case, invariant under a transformation akin to those of the renormalization group (RG) in real space [3,9,10].
The first step of our renormalization-group transformation consists on the raising of the lower bound M c . This implies that only a fraction of events survives the transformation, which defines a different recurrence-time distribution. (If M c is increased in one unit, we are dealing with an authentic decimation, as only about 1 tenth of events are kept, due to the Gutenberg-Richter law.) The second part of the procedure is the scale transformation, which changes the time scale to make the new system comparable with the original one.
The Poisson process, characterized by an exponential recurrence-time distribution, represents a trivial solution to this problem when there are no correlations in the process and therefore events are randomly removed. Indeed, it has been argued that the scaling function f can only be an exponential function [11]. However, the scaling function clearly departs from an exponential [8], and in this way the relevance of correlations in the structure of seismicity becomes apparent.
Indeed, the scaling function f can be described by a decreasing power law for intermediate times, Rτ < 1, with an exponent about 0.3, and a faster decay for long times, Rτ > 1, well approximated in both cases by a gamma distribution. No model of earthquake occurrence is assumed to obtain these results, they are a fundamental characteristic of seismicity. For short times, Rτ < 0.01, the condition of stationarity is usually not fulfilled and the behavior is not universal. Nevertheless, in the non-stationary case the process can be transformed into a stationary one with an appropriate transformation of the time axis, and then the same scaling relation is found to hold again [8].
It should be noted that the term "correlations" can be understood in two forms: If the recurrence-time distribution is not an exponential, this implies the existence of a memory effect in the process, as events do not occur independently at any time, as it would be the case for a Poisson process. But further, the recurrence time may depend on the history of the process, in particular the occurrence time and magnitude of previous events. We shall see that this type of correlations are responsible of the breaking of the memoryless character of seismicity.
In general, the time between two consecutive earthquakes, τ i , may depend on the magnitude of the previous event, M pre = M i−1 , the previous recurrence time, τ i−1 , and also on the occurrence of preceding events, i − 2, i − 3, etc. In their turn, the magnitude of the i−th event M i can depend on τ i , M i−1 , τ i−1 , and so on [12]. We shall only consider here the dependence of τ i with M i−1 , as we have verified it is the most important (together perhaps with the dependence of τ i with τ i−1 ). Note also that although the dependence of the recurrence time and the magnitude with the distance between events can be important, as we have not considered spatial degrees of freedom, we do not need to take this effect into account.
So, in what follows we study the effect in the structure of seismicity of the simplified case of correlations between the recurrence time and the magnitude of the previous earthquake. If we raise the magnitude threshold from M c to M ′ c the distribution of recurrence times for events with magnitudes M ≥ M ′ c can be obtained from the distribution for events with M ≥ M c . Assuming that an event with magnitude M 0 ≥ M ′ c has occurred, we can write for the next event above (or at) M ′ c , where the sum has to be extended up to infinity. P denotes probability and the n−th return time is defined, for events with M ≥ M c , as t i − t i−n , that is, as the elapsed time between any event and the n−th event after it (of course, the 1st return time is the recurrence time). As the recurrence time depends only on the previous magnitude, but the magnitude is independent on any other variable we can write where the last equality comes from the Gutenberg-Richter law . Derivation in Eq. (3) with respect τ yields the probability densities D of the different return times; as the recurrence times are independent on each other, we use that the n−threturn-time distribution is given by n convolutions of the first-return-time distributions (denoted by the symbol * ) to get (5) where ⊤ 1/2 D(τ) denotes the probability density for events with M ≥ M ′ c as a transformation ⊤ 1/2 of the probability density for events with M ≥ M c , D(τ); more precisely, The subscript 1/2 refers to the fact that this is only the first half of the RG transformation. D ↑ and D ↓ denote the recurrence-time probability densities for events above M c conditioned to the fact that the magnitude of the previous event is above or below M ′ c (↑ or ↓), respectively. To be precise, . In Laplace space the things are simpler, Notice that we have used the same symbol D for both the probability densities and for their Laplace transforms (which we may call generating functions), although they are different functions, of course. As q and D ↓ (s) are smaller than one (this is general for generating functions), the infinite sum can be performed, turning out that using that D(τ) is in fact a mixture of the distributions D ↑ and D ↓ , of the form D = pD ↑ + qD ↓ . The previous equation (7) describes the first part of the transformation. The second part is the scale transformation, which puts the distributions corresponding to M c and M ′ c on the same scale. We will obtain this by removing the effect of the decreasing of the rate, which has to be proportional to p, so, and in Laplace space we get Therefore, the combined effect of both transformations leads to In order to get some understanding of this transformation we can consider first the case in which there are no correlations between the magnitude and the subsequent recurrence time. Then, D ↑ = D ↓ = D ≡ D 0 and ⊤D 0 (s) = pD 0 (ps) 1 − qD 0 (ps) .
The fixed points of the transformation are obtained by the solutions of ⊤D 0 (s) = D 0 (s). If we introduce ω ≡ ps and substitute p = ω/s, q = 1 − ω/s, we get, separating variables, where k is a constant, due to the fact that p and s are independent variables and so are s and ω. The solution is then which corresponds to the Laplace transform of an exponential distribution. Indeed, So, in the case in which there are no correlations in the process, the only scale invariant distribution is, as one could have expected, the exponential distribution, characteristic of the Poisson process. Let us see how the existence of correlations between the magnitudes and the recurrence times changes this picture.
Correlations introduce new functions in the process. In our case, in order to iterate the transformation ⊤ we need to know how D ↑ transforms as well. It turns out that D ↑ verifies an equation very similar to Eq. (10), which depends on . So, in order to apply again ⊤ we also need an equation to transform D ⇑ , which in turn depends on higher values of the magnitude threshold. In this way we obtain a hierarchy of equations. An easy way to break this hierarchy is to assume that, at least at the fixed point, D ↑ has the same functional form as D, but in a different scale. So, let us assume where Λ depends on M ′ c − M c with Λ(0) = 1. Figure 2 illustrates this behavior using worldwide earthquakes from the NEIC PDE catalog [8].
In order to find the fixed point of the infinitesimal transformation we impose that the coefficient of δ is zero, obtaining, whose integration yields to where the exponent 1 + ε comes from the definition ε ≡ C/B. We immediately see that in the case of no correlations, C = 0, and therefore ε = 0, recovering Eq. (13) and then the exponential form for D(τ). But there are other values of ε for which the previous equation can be easily solved.
Let us consider first the case ε = 1, corresponding to B = C. The solution of Eq. (19) is whose inverse Laplace transform can be calculated [13], turning out to be, with er f the error function. The asymptotic behavior of D(τ) is very clear: for τ → 0 it diverges as a power law, D(τ) → 1/ √ πkτ, whereas for τ → ∞, D(τ) decays exponentially, using the expansion of the error function [13] and the fact that a power-law factor varies much more slowly than an exponential, for large arguments.
It is interesting also to study the case ε → 0, characteristic of weak correlations, C ≃ 0. We can write the solution D(s) as a perturbation of the Poisson behavior corresponding to ε = 0, i.e., D(s) = D 0 (s) + u(s), with D 0 (s) = (1 + ks) −1 , see Eq. (13). Substituting into Eq. (19) and using (1 + ks) ε ≃ 1 + ε ln(1 + ks) we get and therefore D(s) is just a power of the transform of the exponential density, which means that D(τ) is a gamma distribution, i.e., In this way, for ε ≃ 0, we also get a power law for short times and an exponential decay for long times. This behavior is by no means exclusive of ε = 1 or ε ≃ 0. If in Eq. (19) we consider the limit s → ∞ we get, as D(s) → 0 (which is general for generating functions), and, if 1 + ε > 0, by means of a Tauberian theorem [14], This implies that for all ε > 0, which corresponds to Λ > 1, and then to a shortening of the recurrence times after large earthquakes, we get a decreasing power law, which is a signature of clustering, and in concordance with real data [8]. (The situation would be reversed if −1 < ε < 0.) Further insight can be obtained from Eq. (19), writing it as D = 1 − ksD 1+ε , which can be iterated to get the Taylor expansion of D(s). Up to second order, it turns out that D(s) = 1 − ks + (1 + ε)k 2 s 2 , and as the coefficients of s n yield the moments of the distribution (with a (−1) n /n! factor [14]) we conclude that τ = k, τ 2 = 2(1 + ε)k 2 and the coefficient of variation is cv ≡ τ 2 − τ 2 / τ = √ 1 + 2ε. In fact, the Taylor expansion can be obtained up to any order, and so, all the moments τ n exist and D(τ) decays faster than any power law for τ → ∞, in agreement again with the observations [8].
To conclude, we note that both B and C can be understood as critical exponents. As the energy radiated by an earthquake is about E = E 0 10 3M/2 [3], the Gutenberg-Richter law writes If we define ν = 2C/(3 ln 10), this exponent, the b−value, and the exponent of the recurrence-time distribution for short times given by Eq. (25) fulfill the following scaling relation, ε just using ε = C/B and B = b ln 10. In fact, we must understand this relation as the contribution of the M pre − τ correlations to the recurrence-time distribution, as the value obtained for C (about 0.2) is too small to account for the value of the exponent, about 0.3, using a b−value ≃ 1. We believe the inclusion of other correlations in the calculation will yield to a better quantitative agreement. Our approach mainly consists of a simplification of real seismicity, which allows to understand the complex structure of seismic occurrence in the time and magnitude domains. It is remarkable that simply by imposing the self-similarity of the process and with the only assumption of the scaling of the conditional distributions D ↑ (which is in agreement with the observations), we get such a general characterization of the recurrence-time distribution. With this study we have shown how the structure of seismic occurrence in time, space, and magnitude can be understood as a critical phenomenon and then constitutes a statistical-physics problem [15].