text:
# Diffusion of epicenters of earthquake aftershock, Omori law and generalized continuous-time random walk models

## I Introduction

The spatio-temporal complexity of earthquakes is often invoked as an illustration of the phenomenon of critical self-organization with scale-invariant properties Geilo ; Rundleklein ; Main ; Sorreview ; Turcottereview . This concept points to the importance of developing a system approach in which large scale properties can emerge from the repeating interactions occurring at smaller scales. Such ideas are implemented in models proposing links between the physics of earthquakes and concepts of statistical physics, such as critical points, self-organized criticality, spinodal decomposition, critical depinning, etc., in order to explain the most solidly established facts in the phenomenology of earthquakes, of which we cite the three most important.

* LAW 1: The Gutenberg-Richer law GR  states that the cumulative distribution of earthquake magnitudes $`m`$ sampled over broad regions and large time intervals is proportional to $`10^{bm}`$, with a $`b`$-value $`b1`$. Translating into energies $`E`$ with the correspondence $`m=(2/3)\mathrm{log}_{10}E+`$ constant leads to a power law $`1/E^B`$ with $`B2/3`$.
* LAW 2: Omori law for aftershocks Omori  states that the rate of earthquakes triggered by a mainshock decays with time according to an inverse power $`1/t^p`$ of time with an exponent $`p1`$.
* LAW 3: The earthquakes are clustered in space along hierarchical fault structures Ouillonetal  and their spatial distribution over long times can be approximately described by a fractal dimension close to $`2.2`$ (in three dimensions) KK80 .

There are many other empirical “laws” but these three characterize the very fundamentals of seismicity in size, time and space.

We should immediately point out that these three laws come with significant caveats.

1. There have been on-going controversies on the universality of the exponent $`B`$ or $`b`$-value of the Gutenberg-Richter law Pissor ; kaganuniv .
2. The exponent $`p`$ of Omori law exhibits a large variability from one aftershock sequence to another aftershock sequence and is found typically in the range from $`0.3`$ to $`2`$. We note however that not all these values, especially the extreme ones, automatically reflect a bona-fide power law decay and one should exert caution in attributing too much confidence to them.
3. The view that geological faults and earthquake hypocenters are fractal objects is now recognized to be a naive description of a much more complex reality in which a hierarchy of scales occur with possibly different organizations at different scales Ouillonetal .

In addition, a major difficulty for making progress in modeling and predicting earthquakes is that these three and other laws may be “explained” by a large variety of models, with many different mechanisms. For instance, with respect to the first two laws, we observe the following.

* There are many mechanisms that create a power law distribution of earthquake sizes (see for instance the list of mechanisms described in chapter 14 of Ref. Sorbook ).
* Omori law is essentially a slowly decaying “propagator” describing a long-time memory of past events impacting on the future seismic activity. Such slow power law time decay of the Omori propagator may result from several and not necessarily exclusive mechanisms (see Harrisreview  and references therein): pore-pressure changes due to pore-fluid flows coupled with stress variations, slow redistribution of stress by aseismic creep, rate-and-state dependent friction within faults, coupling between the viscoelastic lower crust and the brittle upper crust, stress-assisted micro-crack corrosion Yamakno ; Leesor , slow tectonic driving of a hierarchical geometry with avalanche relaxation dynamics Huangetal , etc.

The zeroth order description of earthquakes is to consider a single isolated homogeneous fault on which earthquakes are recurrent to accommodate the long-term slow tectonic loading. But faults are not isolated and the most conspicuous observation is that earthquakes interact and influence each other on complex fault structures. Understanding these interactions is essential for understanding earthquakes and fault self-organization. However, the full impact of interactions between earthquakes is still far from being well understood. The simplest and clearest observation of earthquake interaction is provided by aftershocks whose phenomenology is captured by Omori law (LAW 2). Indeed, aftershocks are the most obvious and striking signature of the clustering of the seismicity in time and space, and are observed after all large shallow earthquakes. Most aftershocks are triggered a few hours or days after the mainshock. However, due to the very slow power law decay of the rate of aftershock, known as the Omori law Omori , aftershocks can be triggered up to a hundred years after the mainshock Utsu . Aftershocks often occur near the rupture zone of the mainshock with a variety of focal mechanisms suggesting that they are actually on separate structures BathRichter ; Beroza . They are also sometimes triggered at very large distances from the mainshock Hill ; Steeples ; KJ1 ; Meltzner1 ; Dreger . As an example, Hill et al. Hill  observed aftershocks of the Landers earthquake as far as 1250 km from the epicenter. Similarly to the temporal distribution of aftershocks, a power-law distribution seems to describe well the distribution of distances between pairs of events KJ1  . Since a power-law decays slowly, it describes a slow decay of the probability of observing aftershocks at large distances to the mainshock.

Thus, Omori law can be considered as the simplest and best established description of earthquake interactions of a certain kind. The question we investigate is whether it can be used fruitfully to explain a larger variety of earthquake interactions beyond the class of observations that were used to establish it. In a series of papers SS ; HS1 ; SH2 , we find that Omori law for aftershocks plus the constrain that aftershocks are distributed according to the Gutenberg-Richter power law for earthquake size distribution independently of the magnitude of their progenitor is enough to derive many of the other empirical “laws,” as well the variability of the $`p`$-exponent. Here, we test the potential of this approach to account for and to quantify observations on aftershock diffusion.

Aftershock diffusion refers to the phenomenon of expansion or migration of aftershock zone with time Mogi ; Imoto ; Chatelain ; Tajima1 ; Tajima2 ; Wesson ; Ouchi ; Noir ; Jacques . Immediately after the mainshock occurrence, most aftershocks are located close to the rupture plane of the mainshock, then aftershocks seem to migrate away from the mainshock, at velocity ranging from 1 km/hour to 1 km/year Jacques ; Rydelek . Note that this expansion is not universally observed, but is more important in some areas than in others Tajima1 .

The diffusion of aftershocks is usually interpreted as a diffusion of the stress induced by the mainshock, either by a viscous relaxation process Rydelek , or due to fluid transfer in the crust Nur ; Hudnut ; Noir . Another interpretation of the expansion of aftershocks is given by Dieterich Diete , who reproduces the Omori law decay of aftershocks and the expansion of the aftershock zone with time, using a rate and state friction law and assuming that the rate of aftershocks is proportional to the stress rate. In his model, the expansion of aftershock zone arises from the non-uniform stress induced by the mainshock. Another alternative explanation is that the diffusion of aftershocks is mainly due to the occurrence of large aftershocks, and to the localization of secondary aftershock close to the largest aftershocks, as observed by Ouchi Ouchi . The apparent diffusion of the seismicity may thus result from a cascade process; the mainshock triggers aftershocks that in turn trigger their own aftershocks, and thus lead to an expansion of the aftershock zone.

In the present paper, we investigate the epidemic time aftershock sequence model (ETAS), and show that the cascade of secondary aftershocks can indeed explain the reported diffusion of aftershocks. The ETAS model was introduced by Kagan and Knopoff KK  (in a slightly different form than used here) and Ogata Ogata1  to describe the temporal and spatial clustering of seismicity. This model provides a tool for understanding the clustering of the seismic activity, without arbitrary distinction between aftershocks, foreshocks and mainshocks. In this model, all earthquakes are assumed to be simultaneously mainshocks, aftershocks and possibly foreshocks. Each earthquake generates aftershocks that decay with time according to Omori law, which will in turn generate their own aftershocks. The seismicity rate at any given time and location is given by the superposition of aftershock sequences of all events impacting that region at that time according to space-time “propagators.” The additional ingredient in the version of the ETAS model that we study is that the number of aftershocks per earthquake increases exponentially $`10^{\alpha m}`$ with the magnitude $`m`$ of the mainshock (i.e., as a power law $`E^{2\alpha /3}`$ of the energy released by the mainshock), in agreement with the observations YaShi ; Drakatos . Since the energy of an earthquake is a power law of its rupture length, this law expresses the very reasonable idea that the number of events related to a given earthquake is proportional to a power of its volume of influence. The value of the exponent $`\alpha `$ controls the nature of the seismic activity, that is, the relative role of small compared to large earthquakes. Few studies have measured $`\alpha `$ in seismicity data YaShi ; Guo2 ; H02 . This parameter $`\alpha `$ is often found close to $`b`$ YaShi  or fixed arbitrarily equal to $`b`$ KK ; Felzer . In the case where $`\alpha `$ is close to the Gutenberg-Richter $`b`$-value, this law also reproduces Felzer  the self-similar empirical Bath’ s law Bath , which states that the average difference $`m_Mm_A`$ in size between a mainshock and its largest aftershock is $`1.2`$ magnitude units, regardless of the mainshock magnitude: $`m_A=m_M1.2`$. If $`\alpha <b`$, small earthquakes, taken together, trigger more aftershocks than larger earthquakes. In contrast, large earthquakes dominate earthquake triggering if $`\alpha b`$. This case $`\alpha b`$ has been studied analytically in the framework of the ETAS model by Ref. SH2  and has been shown to eventually lead to a finite time singularity of the seismicity rate. This explosive regime cannot however describe a stationary seismic activity.

A natural way to tame this singular behavior is to introduce an upper cut-off for the magnitude distribution at large magnitudes, mirroring the cut-off $`m_0`$ used for the low-magnitude range. The physical argument for introducing this cut-off is based on the finiteness of the maximum earthquake that the earth is capable of carrying. The specific way of introducing such a cut-off (abrupt or smooth with a transition to a power law with larger exponent or to an exponential taper) is not very important qualitatively because all these laws will regularize the singular behavior and make the average branching ratio finite. Such regularization with a maximum upper magnitude then allows $`\alpha b`$. The special case $`\alpha =b`$ required for Bath’s law to hold exactly can not therefore be excluded.

However, based on a recent re-analysis of seismic catalogs using the powerful collapse technique, one of us H02  has presented strong evidence that $`\alpha `$ is strictly smaller than $`b`$. In this paper, we will therefore consider only the case $`\alpha <b`$ and take $`\alpha =0.5`$ specifically in our numerical simulations. In this regime $`\alpha <b`$, Bath’ s law cannot be reproduced because the average difference in size between a mainshock and its largest aftershock increases with the mainshock magnitude. For $`\alpha <b`$, it is easy to show that Bath’s law is replaced by $`m_A=(\alpha /b)m_M`$ constant, where $`m_M`$ and $`m_A`$ are the magnitudes of the mainshock and of the largest aftershock. Tests of this prediction will be reported in a future publication but we expect that distinguishing this modified Bath’s law from Bath’s law will be a difficult task due to the limited range of the studied magnitudes as well as the dependence of the distribution of $`m_Mm_A`$ on the magnitude thresholds chosen for the mainshocks and for the aftershocks consolebath .

We assume that the distribution of all earthquakes follow the Gutenberg-Richter distribution and take this distribution of aftershock sizes to be independent of the magnitude of the mainshock. Therefore, an earthquake can trigger a larger earthquake, albeit with a small probability. This model can thus describe a priori both aftershock and foreshock sequences. The ETAS model has been calibrated to real seismicity catalogs to retrieve its parameters Ogata1 ; Ogata2 ; Ogata3 ; Ogata5 ; Ogata6 ; Kagan ; Guo2 ; Felzer  and to give short-term probabilistic forecast of seismic activity by extrapolating past seismicity into the future via the use of its space-time propagator KK ; KJ2 ; Console .

The ETAS model is a branching model which exhibits different regimes HS1  depending upon the value of the branching ratio $`n`$, defined by the average number of primary aftershocks per earthquake. The critical case $`n=1`$ corresponds to exactly one primary aftershock per earthquake, when averaging over all mainshock magnitudes larger than a threshold $`m_0`$. Let us stress that $`n`$ is an average quantity which does not reflect adequately the large variability of the number of aftershocks per main shock, as a function of its magnitude. Indeed, the number of aftershocks per mainshock increases exponentially fast as a function of the mainshock magnitude, so that large mainshocks will have significantly more than $`n`$ aftershocks. For $`\alpha =0.5`$, a magnitude $`7`$-earthquake gives typically $`10`$ times more direct aftershocks than a magnitude $`5`$, and $`100`$ times more direct aftershocks than a magnitude $`3`$-earthquake. The increase in triggered seismic activity with the magnitude of the mainshock is obviously stronger for a larger value of $`\alpha `$. Note that these numbers refer to aftershocks of the first generation; the total number of triggered events is larger by the factor $`1/(1n)10`$ (for $`n0.9`$ which is typical), due to the cascades of secondary aftershocks. Notwithstanding this large variability, the average number $`n`$ of primary aftershocks per earthquake controls the global regime. For $`n`$ exactly equal to $`1`$, seismicity is at the border between death and growth. In the sub-critical regime $`n<1`$, since each earthquake triggers on average less that one aftershock, starting from a large event, the seismicity will decrease with time and finally die out. The super-critical $`n>1`$ corresponds to more that one primary aftershock per earthquake on average. Starting from a large earthquake, after a transient regime, the average seismicity will finally increase exponentially with time HS1 , but there is still a finite probability for aftershock sequences to die out.

The numerical simulations reported below have been performed with $`\alpha =0.5`$. It is probable that a good fit to seismic data is obtained by using a value of $`\alpha 0.8`$ larger that the value $`0.5`$, as reviewed and documented recently by one of us H02 . We have checked that results similar to those presented below hold true qualitatively for larger values $`0.5<\alpha <1`$. Such larger values of $`\alpha `$ lead however to stronger fluctuations which are more difficult to handle numerically because the variance of the number $`\rho (m)`$ of direct triggered aftershocks defined below in (3) becomes undefined for $`\alpha >0.5`$. A full understanding of this regime requires a special treatment that will be reported elsewhere.

Sornette and Sornette SS  studied analytically a particular case of this model, without magnitude and spatial dependence, and they considered only the subcritical regime $`n<1`$. Starting with one event at time $`t=0`$ and considering that each earthquake generates an aftershock sequence with a “local” Omori exponent $`p=1+\theta `$, where $`\theta >0`$, they studied the decay law of the “global” aftershock sequence, composed of all secondary aftershock sequences, i.e., by taking into account that the primary aftershocks can create secondary aftershocks which themselves may trigger tertiary aftershocks and so on. They found that the global aftershock rate decays according to an Omori law with an exponent $`p=1\theta <1`$, up to a characteristic time SS ; HS1

$$t^{}=c\left(\frac{n\mathrm{\Gamma }(1\theta )}{|1n|}\right)^{1/\theta },$$
(1)
and then recovers the local Omori exponent $`p=1+\theta `$ for time larger than $`t^{}`$. Helmstetter and Sornette HS1  extended their analysis to the general ETAS model with magnitude dependence, and considered both the sub- and the super-critical regime, but still restricted the analysis to the temporal distribution of the seismicity, without spatial dependence. In the sub-critical regime, they recovered the crossover found by Sornette and Sornette SS . In addition, Helmstetter and Sornette HS1  give the explicit mathematical formula for the gradual transition between the Omori law with exponent $`p=1\theta `$ for $`tt^{}`$ to the Omori law with exponent $`p=1+\theta `$ for $`tt^{}`$. This smooth transition can be observed in figure 2 on the line calculated for $`t^{}=10^9`$ days with $`n<1`$. $`t^{}`$ can thus be viewed as the time where the apparent exponent $`p`$ of the Omori law is approximately in between the two asymptotic values $`1\theta `$ and $`1+\theta `$. A more rigorous mathematical definition HS1  is that $`t^{}`$ is the characteristic time scale such that $`\beta t^{}`$ is the dimensionless variable of the Laplace transform (with variable $`\beta `$) of the seismicity rate.

In the super-critical regime, Helmstetter and Sornette HS1  found a novel transition between a power-law decay with exponent $`p=1\theta `$ at early times, similar to the sub-critical regime, to an exponential increase of the seismicity at large times. The regime where $`\alpha >b`$ or equivalently $`2\alpha /3>B`$ has been found to lead to a new kind of critical stochastic finite-time-singularity SH2 , relying on the interplay between long-memory and extreme fluctuations. Recall that the number of aftershocks per earthquake increases as a power law $`E^{2\alpha /3}`$ of the energy released by the mainshock whereas the number of earthquakes of energy $`E`$ decreases as the Gutenberg-Richter law $`1/E^{1+B}`$. Intuitively, when $`2\alpha /3>B`$, the increase in the rate of creation of aftershocks with the mainshock energy more than compensate the decrease of the probability to get a large mainshock when the mainshock energy increases. This theory based solely on the ETAS model has been found to account for the main observations (power law acceleration and discrete scale invariant structure) of critical rupture of heterogeneous materials, of the largest sequence of starquakes ever attributed to a neutron star as well as of some earthquake sequences SH2 .

In the sequel, we extend the analytical study of the temporal ETAS model SS ; HS1 ; SH2  to the spatio-temporal domain. To model the spatial distribution of aftershocks, we assume that the distance between a mainshock and each of its direct aftershock is drawn from a given distribution, independently of the magnitude of the mainshock and of the delay between the mainshock and its aftershocks. For illustration but without loss of generality for the mapping to the continuous time random walks model (CTRW) discussed later, we shall take a power law distribution of distances between earthquakes. We take the simplest and most parsimonious hypothesis that space, time and magnitude are decoupled in the earthquake propagator. Our first result is to establish a correspondence between the ETAS model and the CTRW, first introduced by Montroll and Weiss Montroll1  and used to model many physical processes. We then build on this analogy to derive the joint probability distribution of the times and locations of aftershocks. We show analytically that, for sufficiently short times $`t<t^{}`$, the average distance between a mainshock and its aftershock increases subdiffusively as $`Rt^H`$, where the exponent $`H`$ depends on the local Omori exponent $`1+\theta `$ and on the distribution of the distances between an earthquake and its aftershocks. We also demonstrate that the local Omori law is not universal, but varies as a function of the distance from the mainshock. Due to the diffusion of aftershocks with time, the decay of aftershock is faster close to the mainshock than at large distances. These non-trivial space-time couplings occur notwithstanding the decoupling between space, time and magnitude in the “bare” propagator, and is due to the existence of cascades of aftershocks.

A recent work of Krishnamurthy et al. Tanguy  substantiates the general modeling strategy used here of representing the space-time dynamics of earthquakes by an effective stochastic process (the ETAS model) entirely defined by two exponents (corresponding to our $`\mu `$ and $`H(\theta ,\mu )`$ defined below), where $`\mu `$ is the exponent of the power law distribution of jumps between successive active sites and $`H`$ is the (sub-)diffusion exponent. Indeed, Krishnamurthy et al. Tanguy  show that the Bak and Sneppen model and the Sneppen model of extremal dynamics (corresponding to a certain class of self-organized critical behavior Sorbook ) can be completely characterized by a suitable stochastic process called “Linear fractional stable motion.” Beyond recovering the scaling exponents of this model, the stochastic process strategy predicts the conditional probabilities of successive activations at different sites and thus offers novel insights. We note that this approach with the Linear fractional stable motion is extremely close in spirit as well as in form to our approach mapping the ETAS model to the CTRW model. The ETAS model can thus be taken to represent an effective stochastic process of the complex self-organization of seismicity.

## II The ETAS model

### II.1 Definitions and specific parameterization of the ETAS model

We assume that a given event (the “mother”) of magnitude $`m_i`$ occurring at time $`t_i`$ and position $`\stackrel{}{r}_i`$ gives birth to other events (“daughters”) of any possible magnitude chosen with some independent Gutenberg-Richter distribution at a later time between $`t`$ and $`t+dt`$ and at point $`\stackrel{}{r}\pm \stackrel{}{d}r`$ to within $`d\stackrel{}{r}`$ at the rate

$$\varphi _{m_i}(tt_i,\stackrel{}{r}\stackrel{}{r}_i)=\rho (m_i)\mathrm{\Psi }(tt_i)\mathrm{\Phi }(\stackrel{}{r}\stackrel{}{r}_i).$$
(2)
We will refer to $`\varphi _{m_i}(tt_i,\stackrel{}{r}\stackrel{}{r}_i)`$ both as the seismic rate induced by a single mother or as the “bare propagator”. It is the product of three independent contributions:

1. $`\rho (m_i)`$ gives the number of daughters born from a mother with magnitude $`m_i`$. This term will in general be chosen to account for the fact that large earthquakes have many more triggered events than small earthquakes. Specifically, we take

$$\rho (m_i)=K10^{\alpha (m_im_0)},$$
(3)
which, as we said earlier, is justified by the power law dependence of the volume of stress perturbation as a function of the earthquake size. $`\alpha `$ quantifies how fast the average number of daughters per mother increases with the magnitude of the mother.
2. $`\mathrm{\Psi }(tt_i)`$ is a normalized waiting time distribution giving the rate of daughters born at time $`tt_i`$ after the mother. The normalization condition reads $`_0^+\mathrm{}𝑑t\mathrm{\Psi }(t)=1`$. $`\mathrm{\Psi }(tt_i)dt`$ can thus be interpreted as the probability for a daughter to be born between $`t`$ and $`t+dt`$ from the mother that was born at time $`t_i`$. $`\mathrm{\Psi }(tt_i)`$ embodies Omori law: it is the “bare” or “direct” Omori law

$$\mathrm{\Psi }(t)=\frac{\theta c^\theta }{(t+c)^{1+\theta }}H(t),$$
(4)
where $`\theta >0`$ and $`H(t)`$ is the Heaviside function.
3. $`\mathrm{\Phi }(\stackrel{}{r}\stackrel{}{r}_i)`$ is a normalized spatial “jump” distribution from the mother to each of her daughter, quantifying the probability for a daughter to be triggered at a distance $`|\stackrel{}{r}\stackrel{}{r}_i|`$ from the mother. Specifically, we take

$$\mathrm{\Phi }(\stackrel{}{r})=\frac{\mu }{d(\frac{|\stackrel{}{r}|}{d}+1)^{1+\mu }},$$
(5)
which has the form of an (isotropic) elastic Green function dependence describing the stress transfer in an elastic upper crust. The exponent $`\mu `$ is left adjustable to account for heterogeneity and the possible complex modes of stress transfers. The normalization condition reads $`𝑑\stackrel{}{r}\mathrm{\Phi }(\stackrel{}{r})=1`$ where the integral is carried out over the whole space.

The physical justification for this decoupled model (2) in which $`\varphi _{m_i}(tt_i,\stackrel{}{r}\stackrel{}{r}_i)`$ is the product of three independent distributions is that elastic waves propagate at kilometers per second and thus almost instantaneously reset the stress field after a large main shock. In other words, there is a well-defined separation of time scales between the time of propagation of seismic waves (seconds to minutes) which control the convergence to a new mechanical equilibrium after the main shock and the time scales involved in aftershock sequences (hours, days, months to many years). The spatial dependence in (2) reflects the stress redistribution. This new stress field then relaxes slowly and more or less independently from point to point leading to the local Omori law $`\mathrm{\Psi }(tt_i)`$. Notwithstanding this argument, the decoupling in (2) between the local responses in magnitudes, space and time is mostly performed because of its simplicity. It constitutes an approximation that should be checked and relaxed in future studies.

We assume a distribution $`P(m)`$ of earthquake sizes expressed in magnitudes $`m`$ which follows the Gutenberg-Richter distribution

$$P(m)=b\mathrm{ln}(10)10^{b(mm_0)},$$
(6)
with a $`b`$-value usually close to $`1`$. $`m_0`$ is a lower bound magnitude below which no daughter is triggered.

### II.2 The branching ratio $`n`$

A key parameter of the ETAS is the average number $`n`$ of daughter-earthquakes created per mother-event, summed over all possible magnitudes. As we shall see, it is also natural to call it the “branching ratio”. To see this, consider the integral of the seismic rate $`\varphi _{m_i}(tt_i,\stackrel{}{r}\stackrel{}{r}_i)`$ induced by one earthquake over all times after $`t_i`$, over all spatial positions and over all magnitudes $`m_im_0`$, which must give by definition the average number $`n`$ of direct (or primary) daughter-earthquakes created per mother-event independently of its magnitude. For $`\alpha <b`$ and using (2), (3) and (6), it is exactly given by

$$n𝑑\stackrel{}{r}_{t_i}^+\mathrm{}𝑑t_{m_0}^+\mathrm{}𝑑m_iP(m_i)\varphi _{m_i}(tt_i,\stackrel{}{r}\stackrel{}{r}_i)=_{m_0}^+\mathrm{}𝑑m_iP(m_i)\rho (m_i)=\frac{Kb}{b\alpha },$$
(7)
since the two integrals over time and space contribute each a factor $`1`$ by the normalization of $`\mathrm{\Psi }`$ and $`\mathrm{\Phi }`$. This result (7) is identical to that found in absence of spatial dependence of $`\varphi _{m_i}(tt_i)`$ with respect to $`\stackrel{}{r}\stackrel{}{r}_i`$ due to the factorization of the rate $`\rho `$, time $`\mathrm{\Psi }`$ and space $`\mathrm{\Phi }`$ dependences HS1 . The branching ratio has also been evaluated in the case where the magnitude distribution follows a gamma distribution Kagan .

We stress again that $`n`$ is an average quantity which does not reflect the large fluctuations in the number of aftershocks from events to events. Indeed, large events with magnitudes $`M`$ produce in general many more aftershocks than small events with magnitude $`m<M`$, simply because $`\rho (M)\rho (m)`$ if $`M>m`$ (see the exponential dependence (3) of $`\rho (m)`$ on the magnitude $`m`$).

### II.3 Numerical simulation of the spatial ETAS model

The ETAS model has been simulated numerically using the algorithm described in Refs. Ogata4 ; Ogata5 . Starting with a large event of magnitude $`M`$ at time $`t=0`$, events are then simulated sequentially. At any given time $`t`$, we calculate the conditional seismic rate $`\lambda (t)`$ defined by

$$\lambda (t)=\underset{t_it}{}K10^{\alpha (m_im_0)}\frac{\theta c^\theta }{(tt_i+c)^{1+\theta }}$$
(8)
where $`K=n(b\alpha )/b`$, and $`t_i`$ and $`m_i`$ are the times and magnitudes of all preceding events that occurred at time $`t_it`$. Note that we use the bare propagator because the sum in (8) is performed exhaustively on the complete catalog of past events. The time of the following event is then determined according to the non-stationary Poisson process of conditional intensity $`\lambda (t)`$, and its magnitude is chosen according to the Gutenberg-Richter distribution with parameter $`b`$. To determine the position in space of this new event, we first choose its mother randomly among all preceding events with a probability proportional to their rate of aftershocks $`\varphi _{m_i}(tt_i)`$ evaluated at the time of the new event. Once the mother has been chosen, we generate the distance $`r`$ between the new earthquake and its mother according to the power-law distribution $`\mathrm{\Phi }(\stackrel{}{r})`$ given by (5). The location of the new event is determined by assuming an isotropic distribution of aftershocks. By this rule, it is clear that new events tend to be close in general to the last large earthquakes, leading to space clustering.

Note that this two-steps procedure is equivalent to but more convenient for a numerical implementation than the one-step method, consisting of calculating at each point on a fine space-covering grid the seismic rate, equal to the sum over all preceding mothers weighted by the bare space $`\mathrm{\Phi }(\stackrel{}{r})`$ and time $`\mathrm{\Psi }(t)`$ propagators given by (5) and (4); after normalizing, these rates then provide to each grid point a probability for the event to occur on that point. The equivalence between our two-step procedure and the direct calculation of the seismic rates is based on the law of conditional probabilities: probability of next event ($`A`$) = probability of next event conditioned on its mother (event $`B`$) $`\times `$ probability of choosing the mother, i.e., $`P(A,B)=P(A|B)\times P(B)`$.

Figure 1 shows the result of a numerical simulation of the ETAS model which exhibits a diffusion of the seismic activity. We simulate a sequence of aftershocks and secondary aftershocks starting from a mainshock of magnitude $`M=7`$, with the following parameters: $`\theta =0.2`$, $`b=1`$, $`\alpha =0.5`$, $`n=1`$ and $`\mu =1`$. At early times, aftershocks are localized close to the mainshock, and then diffuse and cluster close to the largest aftershocks. This (sub-)diffusion is extremely slow, as we shall quantify in the sequel. Our purpose is to provide a theory for this process based on the ETAS model. This theory will be tested by numerical simulations.

The different regimes are illustrated by Figure 2 which shows the seismicity rate $`N(t)`$ for the temporal ETAS model studied by SS ; HS1  obtained by summing the seismic activity over all space for the 3 cases $`n<1`$ (sub-critical), $`n=1`$ (critical) and $`n>1`$ (super-critical). The sub-critical regime is characterized by the existence of the time scale $`t^{}`$ given by (1). There is no difference between the critical case $`n=1`$ and the sub-critical case for $`t<t^{}`$ (see figure 2). Indeed, the difference between the sub-critical regime and the critical regime can be observed only for $`t>t^{}`$. A simple way to see this is to realize that the critical regime $`n=1`$ gives $`t^{}=+\mathrm{}`$, meaning that, in the critical regime, one is always in the situation $`t<t^{}`$.

It is interesting to note that the spatial distribution of epicenters shown in the right panel of figure 1 has the visual appearance of a fractal set of points. This is confirmed by the calculation of the correlation dimension of this set of $`N=3000`$ points generated in the time interval $`[30,70]`$ yrs, which is found approximately equal to $`D_2=1.5\pm 0.05`$ over more than two decades in spatial scales, as shown in figure 3. If we use instead all $`30,000`$ events of the simulation performed up to time $`t=70`$ yrs, we find $`D_2=1.85\pm 0.05`$ while the correlation dimension of the geometrical set made of the epicenters of the $`10,000`$ last events (time interval $`[7,70]`$ yrs) is $`D=1.7\pm 0.05`$, also over more than two decades in scale. These values are similar to those reported for 2D maps of active fault systems ScholzMandel ; Anneexp ; BartonLapointe1 ; Ouillonetal , and are in good agreement with $`D_2`$ values in the range $`[1.65,1.95]`$ measured for aftershocks epicenters Nanjo . The fractal clustering of the earthquake epicenters, according to the ETAS model, occurs because of a self-similar process taking place on many different scales. However, the description of this multi-scale process solely in terms of a single fractal dimension fails to fully embody the complex spatial superposition of local “singularities” associated with each aftershock on the one hand and finite-size effects (stemming from the finite lifetime of each aftershock sequence) on the other hand. Each event indeed creates its cloud of direct aftershocks which can be characterized by its singular exponent $`1\mu `$ for $`\mu 1`$ and $`0`$ for $`\mu >1`$, defined by the scaling $`_0^Rr𝑑r/r^{1+\mu }R^{1\mu }`$ of the “mass” of the cloud with its radius $`R`$. Finite-size effects and randomness have been documented to generate realistic but sometimes spurious fractal signatures Ouillonsor ; Hamburger ; Eneva ; Malcai . This problem requires a special study which is left for another work.

### II.4 Relationship with the space-independent ETAS model

The spatial ETAS model reduces to the space-independent ETAS model solved in HS1  by integrating the dressed propagator obtained below over all space. In the Fourier representation (see expression (21)), this corresponds to putting the wavenumber $`k`$ to zero. Indeed, for $`k=0`$, the Fourier transform amounts to perform a simple integration over all space. Since $`\widehat{\mathrm{\Phi }}(\stackrel{}{k}=\stackrel{}{0})=1`$, expression (21) derived below reduces to the form studied at length in HS1 . Therefore, all results reported previously hold also for the version of the space-dependent ETAS model studied here, when averaging over the whole space. This is an important property that all the solutions discussed below must obey.

## III Mapping of the ETAS model on the CTRW model

In order to study the space-time properties of the ETAS model, it is very useful to use an exact correspondence between the ETAS model and the continuous time random walk (CTRW) that we establish here. In this way, we can adapt and use the wealth of results previously derived for the CTRW. But first, let us demonstrate the correspondence between the ETAS and CTRW models. For this, our strategy is to derive the Master equation for both models and show that they are identical.

### III.1 The Master equation of the ETAS model

The ETAS model can be rephrased by defining the rate $`\varphi _{m_im}(tt_i,\stackrel{}{r}\stackrel{}{r}_i)`$ at which a given event (the “mother”) of magnitude $`m_im_0`$ occurring at time $`t_i`$ and position $`\stackrel{}{r}_i`$ gives birth to other events (“daughters”) of specified magnitude $`m`$ at a later time between $`t`$ and $`t+dt`$ and at point $`\stackrel{}{r}`$ to within an infinitesimal volume $`|d\stackrel{}{r}|`$. Note that the only difference with respect to the previous definition (2) is that we now specify also the magnitude $`m`$ of the daughter. $`\varphi _{m_im}(tt_i,\stackrel{}{r}\stackrel{}{r}_i)`$ is given by

$$\varphi _{m_im}(tt_i,\stackrel{}{r}\stackrel{}{r}_i)=\rho (m_im)\mathrm{\Psi }(tt_i)\mathrm{\Phi }(\stackrel{}{r}\stackrel{}{r}_i),$$
(9)
where $`\mathrm{\Psi }(tt_i)`$ and $`\mathrm{\Phi }(\stackrel{}{r}\stackrel{}{r}_i)`$ are the same as previously while

$$\rho (m_im)=P(m)\rho (m_i).$$
(10)
With the parameterization (3) and (6), this reads

$$\rho (m_im)=n\mathrm{ln}(10)(b\alpha )10^{\alpha (m_im_0)}10^{b(mm_0)}.$$
(11)
Let us consider the case where there is an origin of time $`t=0`$ at which we start recording the rate of earthquakes, assuming that a large earthquake has just occurred at $`t=0`$ and somehow reset the clock. In the following calculation, we will forget about the effect of events at times prior to $`t=0`$ and count all aftershocks that are created only by this main shock.

Let us call $`N_m(t,\stackrel{}{r})dtdmd\stackrel{}{r}`$ the number of earthquakes occurring between $`t`$ and $`t+dt`$ of magnitude between $`m`$ and $`m+dm`$ inside of box of volume $`|d\stackrel{}{r}|`$ centered at point $`\stackrel{}{r}`$. $`N_m(t,\stackrel{}{r})`$ is the solution of a self-consistency equation that formalizes mathematically the following process : an earthquake may trigger aftershocks; these aftershocks may trigger their own aftershocks, and so on. The rate of seismicity at a given time $`t`$ and position $`\stackrel{}{r}`$ is the result of this cascade process. The self-consistency equation that sums up this cascade reads

$$N_m(t,\stackrel{}{r})=S(t,\stackrel{}{r},m)+\stackrel{}{d}r^{}_{m_0}^{\mathrm{}}𝑑m^{}_0^t𝑑\tau \varphi _{m^{}m}(t\tau ,\stackrel{}{r}\stackrel{}{r}^{})N_m^{}(\tau ,\stackrel{}{r}^{}).$$
(12)
The rate $`N_m(t,\stackrel{}{r})`$ at time $`t`$ and position $`\stackrel{}{r}`$ is the sum over all induced rates from all earthquakes of all possible magnitudes that occurred at all previous times and locations propagated to the present time $`t`$ and to the position $`\stackrel{}{r}`$ of observation by the corresponding bare propagator. The induced rate of events per earthquake that occurred at an earlier time $`\tau `$ and position $`\stackrel{}{r}^{}`$ is equal to $`\varphi _{m^{}m}(t\tau ,\stackrel{}{r}\stackrel{}{r}^{})`$. The source term $`S(t,\stackrel{}{r})`$ is the main shock plus the background seismicity, if any. In absence of background seismicity, a main earthquake which occurs at the origin of time $`t=0`$ at position $`\stackrel{}{r}=\stackrel{}{0}`$ with magnitude $`M`$ gives

$$S(t,\stackrel{}{r},m)=\delta (t)\delta (mM)\delta (\stackrel{}{r})$$
(13)
where $`\delta `$ is the Dirac distribution. Other arbitrary source functions can be chosen.

The source term corresponding to a single mainshock is indeed the delta function (13) rather than the direct Omori law created by this mainshock in direct lineage. To see this, notice that the direct Omori law is recovered from (12) by replacing $`N_m^{}(\tau ,\stackrel{}{r}^{})`$ in the integral by $`S(t,\stackrel{}{r},m)`$ given by (13). This shows that the difference between the renormalized and the direct Omori laws comes from taking into account the secondary, tertiary, etc., cascade of aftershocks.

As we have seen, a key assumption of the ETAS model is that the daughters born from a given mother have their magnitude drawn independently of the magnitude of the mother and of the process that give them birth, with a probability given by the Gutenberg-Richter distribution (6). The consequences resulting from relaxing this hypothesis will be reported elsewhere. Keeping this assumption, it can be shown HSG  that for $`\alpha b/2`$ an ensemble of realizations will obey

$$N_m(t,\stackrel{}{r})=P(m)N(t,\stackrel{}{r}),\mathrm{for}t>0,$$
(14)
which makes explicit the separation of the magnitude from the time and space variables. $`N(t,\stackrel{}{r})`$ is the number of events at position $`\stackrel{}{r}`$ at time $`t`$ of any possible magnitude. Expression (14) means that the Gutenberg-Richter distribution is preserved at all times. That (14) holds for the ETAS model stems from the fact that the waiting time $`\mathrm{\Psi }(t)`$ distribution (4) and jump size $`\mathrm{\Phi }(\stackrel{}{r})`$ distribution (5) are independent of the magnitudes and that fluctuations in the seismicity rate are not too wild for $`\alpha b/2`$. Note that, in a more complex model in which time, space and magnitudes are interdependent, expression (14) would become a mean-field approximation, in which the fluctuations of the rates induced by the fluctuations of the realized magnitudes of the daughters factorize from the process.

Putting (14) in (12), for $`t>0`$ when the source term $`S(t,\stackrel{}{r},m)`$ is identically zero, one can simplify by $`P(m)`$ and obtain

$$N(t,\stackrel{}{r})=𝑑\stackrel{}{r}^{}_0^t𝑑\tau \varphi (t\tau ,\stackrel{}{r}\stackrel{}{r}^{})N(\tau ,\stackrel{}{r}^{}),t>0,$$
(15)
where

$$\varphi (t\tau ,\stackrel{}{r}\stackrel{}{r}^{})=_{m_0}^{\mathrm{}}𝑑m^{}P(m^{})\varphi _m^{}(t\tau ,\stackrel{}{r}\stackrel{}{r}^{}).$$
(16)
Equation (15) is nothing but the expectation (or statistical average, i.e., average over an ensemble of realizations) of expression (8), with the definition $`N(t,\stackrel{}{r})\mathrm{E}[\lambda (t)\mathrm{\Phi }(\stackrel{}{r})]`$. Therefore, the Master equation obtained here gives us only the first moment of the space-time dynamics of seismicity. It is not difficult to derive the equations for the variance and covariance of the seismic rate as well as higher moments.

The value of the source term at $`t=0`$ that should be incorporated in (15) requires more care. Indeed, a naive treatment would give a source term $`\delta (t)\delta (mM)\delta (\stackrel{}{r})/P(M)`$ obtained by simply dividing by $`P(m)`$, expressed at $`m=M`$ due to the Dirac distribution $`\delta (mM)`$. However, this source term still depends on $`m`$ via the Dirac distribution $`\delta (mM)`$ and is thus unsuitable as a source term of the equation (15) which is independent of $`m`$. In order to circumvent this difficulty, one has to get rid of the Dirac distribution $`\delta (mM)`$. The corresponding procedure has been described in details in Ref. HS1  and consists in applying the integral operator $`_{m_0}^{\mathrm{}}𝑑m\widehat{\varphi }(\beta ,\stackrel{}{r})`$ to (12), where $`\widehat{\varphi }(\beta ,\stackrel{}{r})`$ is the Laplace transform with respect to the time variable of $`\varphi (t,\stackrel{}{r})`$. In this way, the Dirac distribution $`\delta (mM)`$ is regularized. Identifying with the results of Ref. HS1 , we obtain that $`N(t,\stackrel{}{r})`$ is solution of (15) with a source term

$$S_M(t,\stackrel{}{r})=\delta (r)\delta (t)\rho (M)/n,$$
(17)
where $`\rho (M)`$ is defined in (3) and $`n`$ is given by (7). Thus, the complete Master equation for the number $`N(t,\stackrel{}{r})`$ of events at position $`\stackrel{}{r}`$ at time $`t`$ of any possible magnitude is solution of

$$N(t,\stackrel{}{r})=S_M(t,\stackrel{}{r})+𝑑\stackrel{}{r}^{}_0^t𝑑\tau \varphi (t\tau ,\stackrel{}{r}\stackrel{}{r}^{})N(\tau ,\stackrel{}{r}^{}),t>0,$$
(18)
$`N(t,\stackrel{}{r})`$ is the “dressed” or “renormalized” propagator, obtained by summing the bare Omori propagator over all possible aftershock cascades. $`N(t,\stackrel{}{r})`$ can also be called the renormalized Omori law SS .

The essential assumption used to derive (12) is that the fluctuations of the earthquake magnitudes in a given sequence can be considered to be decoupled from those of the seismic rate. This approximation can be shown to be valid for $`\alpha b/2`$ HSG , for which the random variable $`\rho (m_i)`$ has a finite variance. In this case, any coupling between the fluctuations of the earthquake energies and the instantaneous seismic rate provides only sub-dominant corrections to the equation (12). For $`\alpha >b/2`$, the variance of $`\rho (m_i)`$ is mathematically infinite or undefined as $`\rho (m_i)`$ is distributed according to a power law with exponent $`b/\alpha <2`$. In this case, the Master equation (12) is not completely correct as an additional term must be included to account for the effect of the dependence between the fluctuations of earthquake magnitudes and the instantaneous seismic rate. Our results are presented below for $`\alpha =0.5`$ which belongs to the first regime $`\alpha b/2`$. For $`\alpha >b/2`$, Ref. HSG  has shown that the renormalization of the bare propagator into the dressed propagator is weaker than for $`\alpha b/2`$, all the more so as $`\alpha b`$. Preliminary numerical simulations for $`\alpha >b/2`$ shows that our results presented below hold qualitatively but with a reduction of the observed spatial diffusion exponent compared to the value predicted from the Master equation approach developed here. This regime $`\alpha >b/2`$ is probably relevant to the real seismicity YaShi ; Guo2 ; H02 , even if a precise estimation of $`\alpha `$ is very difficult.

### III.2 A Master equation of the CTRW model

We now demonstrate that the self-consistent mean field equation (18) is identical to the Master equation of a continuous-time random walk (CTRW). Random walks underlie many physical processes and are often the basis of first-order description of natural processes. The CTRW model, which is a generalization of the naive model of a random walker which jumps by $`\pm 1`$ spatial step on a discrete lattice at each time step, was introduced by Montroll1  and investigated by many other workers Montroll12 ; Scher ; Kenkre ; Shlesinger1 ; Weiss1 . The CTRW considers a continuous distribution of spatial steps as well as time steps (which can be seen either as waiting times between steps or as durations of the steps). The CTRW model is thus based on the idea that the length of a given jump, as well as the waiting time $`\tau _i=t_it_{i1}`$ elapsing between two successive jumps are drawn from a joint probability density function (pdf) $`\varphi (\stackrel{}{r},t)`$, which is usually referred to as the jump pdf. From a mathematical point of view, a CTRW is a process subordinated to random walks under the operational time defined by the process $`\{t_i\}`$.

From $`\varphi (\stackrel{}{r},t)`$, the jump length pdf $`\mathrm{\Phi }(\stackrel{}{r})=_0^+\mathrm{}𝑑t\varphi (\stackrel{}{r},t)`$ and the waiting time pdf $`\mathrm{\Psi }(t)=𝑑\stackrel{}{r}\varphi (\stackrel{}{r},t)`$ can be deduced. Thus, $`\mathrm{\Phi }(\stackrel{}{r})d\stackrel{}{r}`$ produces the probability for a jump length in the interval $`(\stackrel{}{r},\stackrel{}{r}+d\stackrel{}{r})`$ and $`\mathrm{\Psi }(t)dt`$ the probability for a waiting time in the interval $`(t,t+dt)`$. When the jump length and waiting time are independent random variables, this corresponds to the decoupled form $`\varphi (\stackrel{}{r},t)=\mathrm{\Psi }(t)\mathrm{\Phi }(\stackrel{}{r})`$. If both are coupled, a jump of a certain length involves a time cost or, vice versa in a given time span the walker can only travel a maximum distance. With these definitions, a CTRW process can be described through a Master equation (see Weiss1 ; Hughes ; Meltzner2  for a review and references therein) which turns out to be given by an equation which is identical to (18).

This connection between the ETAS model of earthquakes and a model of random walks provides an important advance for the understanding of spatio-temporal earthquake processes, as it allows one to borrow for the deep knowledge accumulated in past decades on random walks. In the same spirit, polymer physics acquired its status as a fundamental physical problem from its previous status of an applied field of research in chemistry when Flory, Edwards, de Gennes, des Cloizeaux and others showed how to formulate problems in polymer physics in the language of random walks and how to extract novel results. In the sequel of this article, we use this analogy to provide a wealth of new predictions as well as new questions for earthquake aftershocks.

In the context of the CTRW, we have the following correspondence.

* $`N(t,\stackrel{}{r})`$ is the pdf for the random walker to just arrive at position $`\stackrel{}{r}`$ at time $`t`$.
* The source term $`S_M(t,\stackrel{}{r})`$ given by (17) denotes the initial condition of the random walk, here chosen to be at the origin of space at time $`t=0`$. The constant $`\rho (M)/n`$ adds the possibility via the parameter $`M`$ to have more than one initial walker at the origin.
* In the CTRW context, the Master equation (18) states that the pdf $`N(t,\stackrel{}{r})`$ of just having arrived at position $`\stackrel{}{r}`$ at time $`t`$ comes from all possible paths in number $`N(\tau ,\stackrel{}{r}^{})`$ having crossed a position $`\stackrel{}{r}^{}`$ at an earlier time $`\tau `$, weighted by a transfer or propagator function $`\varphi (t\tau ,\stackrel{}{r}\stackrel{}{r}^{})`$ describing all the possible steps of the random walker from $`(\tau ,\stackrel{}{r}^{})`$ to $`(t,\stackrel{}{r})`$.

It is important to stress that $`N(t,\stackrel{}{r})`$ defined above is different from the standard quantity $`W(t,\stackrel{}{r})`$ usually studied in random walk problems, defined as the probability to find the random walk at position $`\stackrel{}{r}`$ at time $`t`$. The relationship between $`N(t,\stackrel{}{r})`$ and $`W(t,\stackrel{}{r})`$ is

$$W(t,\stackrel{}{r})=_0^t𝑑t^{}\left[1_0^{tt^{}}𝑑t^{\prime \prime }\mathrm{\Psi }(t^{\prime \prime })\right]N(t^{},\stackrel{}{r}).$$
(19)
The term $`1_0^{tt^{}}𝑑t^{\prime \prime }\mathrm{\Psi }(t^{\prime \prime })`$ in bracket is the probability for the walker not to jump in the time interval $`[t^{},t]`$ and the integral in the right-hand-side of (19) means that the probability $`W(t,\stackrel{}{r})`$ for the random walker to be at position $`\stackrel{}{r}`$ at time $`t`$ is the sum over all possible scenarios in which the walker just arrives at $`\stackrel{}{r}`$ at an earlier time $`t^{}`$ and then does not jump until time $`t`$. In the context of earthquake aftershocks, $`W(t,\stackrel{}{r})`$ is the probability that an event at $`\stackrel{}{r}`$ has occurred at a time $`t^{}t`$ and that the whole system has remained quiescent from $`t^{}`$ to $`t`$.

In the Fourier-Laplace domain (see below), expression (19) reads

$$\widehat{W}(\beta ,\stackrel{}{k})=\frac{1\widehat{\mathrm{\Psi }}(\beta )}{\beta }\widehat{N}(\beta ,\stackrel{}{k}).$$
(20)
In general, the CTRW models transport phenomena in any heterogeneous media. It has for instance been used successfully for describing the behavior of chemical species as they migrate through porous media Margolin ; Berko1 . In insight, it is rather natural that it can be applied to the “transport of stress” through the heterogeneous crust and thus to the description of the anomalous diffusion of seismic activity.

Table 1 synthesizes the correspondence between the ETAS and CTRW models and then draws its consequences.

### III.3 Experimental verifications of the cross-over between the two power law Omori decays in photoconductivity in amorphous semi-conductors and in fractal stream chemistry using the correspondence between the ETAS and CTRW model

The crossover from an Omori law $`1/t^{1\theta }`$ for $`t<t^{}`$ to $`1/t^{1+\theta }`$ for $`t>t^{}`$ found in SS ; HS1  with $`t^{}`$ given by (1) has actually a counterpart in the CTRW. This behavior was first studied by Scher and Montroll Scher  in a CTRW with absorbing boundary condition to model photoconductivity in amorphous semi-conductors As<sub>2</sub>Se<sub>3</sub> and an organic compound TNF-PVK finding $`\theta 0.5`$ and $`\theta =0.8`$ respectively. In a semiconductor experiment, electric holes are injected near a positive electrode and then transported to a negative electrode where they are absorbed. The transient current follows exactly the transition $`1/t^{1\theta }`$ for $`t<t^{}`$ to $`1/t^{1+\theta }`$ for $`t>t^{}`$ found for Omori law for earthquake aftershocks in the ETAS model. In the semiconductor context, the finiteness of $`t^{}`$ results from the existence of a force applied to the holes while in the ETAS model it results from a finite distance $`1n`$ to the critical point $`n=1`$ in the subcritical regime. When the force goes to zero or $`n1`$, $`t^{}+\mathrm{}`$.

A similar transition has been recently proposed to model long-term time series measurements of chloride, a natural passive tracer, in rainfall and runoff in catchments Schermarklaber . The quantity analogous to the dressed Omori propagator is the effective travel time distribution $`h(t)`$ which governs the global lag time between injection of the tracer through rainfall and outflow to the stream. $`h(t)`$ has been shown to have a power-law form $`h(t)1/t^{1m}`$ with $`m`$ between -0.3 and 0.2 for different time series Kirchner . This variability may be due to the transition between an exponent $`1\theta `$ at short times to $`1+\theta `$ at long times Schermarklaber , where $`\theta `$ is the exponent of the bare distribution of individual transition times.

### III.4 General and formal solution of the spatial ETAS model

Let us solve (18) for the number $`N(t,\stackrel{}{r})`$ of events at position $`\stackrel{}{r}`$ at time $`t`$ of any possible magnitude. Recall that $`N(t,\stackrel{}{r})`$ can also be interpreted as the dressed Omori propagator. Extending HS1  to the spatial domain and also in analogy with the standard approach to solve the CTRW, the Laplace-in-time Fourier-in-space transform $`\widehat{N}(\beta ,\stackrel{}{k})`$ of $`N(t,\stackrel{}{r})`$ is given by

$$\widehat{N}(\beta ,\stackrel{}{k})=\frac{\widehat{S}_M(\beta ,\stackrel{}{k})}{1n\widehat{\mathrm{\Psi }}(\beta )\widehat{\mathrm{\Phi }}(\stackrel{}{k})},$$
(21)
where $`\widehat{S}_M(\beta ,\stackrel{}{k})`$ is the Laplace Fourier transform of the source $`S_M(t,\stackrel{}{r})`$ given by (17) and $`\widehat{\mathrm{\Psi }}(\beta )`$ (respectively $`\widehat{\mathrm{\Phi }}(\stackrel{}{k})`$) is the Laplace (respectively Fourier) transforms of $`\mathrm{\Psi }(t)`$ (respectively $`\mathrm{\Phi }(\stackrel{}{r})`$). For a mainshock of magnitude $`M`$ occurring at time $`t=0`$ and position $`\stackrel{}{r}=0`$, the source term is thus $`\widehat{S}_M(\beta ,\stackrel{}{k})=\rho (M)/n`$. The only difference between expression (21) and the Laplace-Fourier transform of the pdf of the CTRW of just having arrived at $`\stackrel{}{r}`$ at time $`t`$ occurs when the branching ratio $`n`$ is different from $`1`$. In general, solutions of CTRW models are expressed for $`n=1`$ and for the variable $`W(t,\stackrel{}{r})`$ which is simply related to $`N(t,\stackrel{}{r})`$ according to (19). Using (19) and (21) leads to

$$\widehat{W}(\beta ,\stackrel{}{k})=\frac{1\widehat{\mathrm{\Psi }}(\beta )}{\beta }\frac{\widehat{S}_M(\beta ,\stackrel{}{k})}{1n\widehat{\mathrm{\Psi }}(\beta )\widehat{\mathrm{\Phi }}(\stackrel{}{k})},$$
(22)
In the following, we exploit (22) to obtain analytical solutions of the spatial ETAS model in different regimes, that provide specific predictions on the conditions necessary for observing aftershock diffusion. In addition, we provide specific predictions on the exponent $`H`$ of the diffusion law $`Rt^H`$ that are tested by numerical simulations.

## IV Critical regime $`n=1`$

### IV.1 Classification of the different regimes

Numerous works on the CTRW have investigated many possible forms for $`\mathrm{\Psi }(t)`$ and $`\mathrm{\Phi }(\stackrel{}{r})`$ and have provided the asymptotic long time and large scale dependence of $`W(t,\stackrel{}{r})`$ (see Weiss1 ; Hughes ; Meltzner2 ; Berko1  and references therein). Here, we restrict our discussion to the cases where both $`\mathrm{\Psi }(t)`$ and $`\mathrm{\Phi }(\stackrel{}{r})`$ have power law tails as given by (4) and (5). The long-time and large scale behavior of the ETAS and CTRW are controlled by the behavior of the Laplace-Fourier transforms for small $`\beta `$ and small $`|\stackrel{}{k}|`$.

Two cases must be distinguished depending on the exponent $`\mu `$ controlling the weight of the tail of $`\mathrm{\Phi }(\stackrel{}{r})`$.

* For $`\mu >2`$, the variance $`(\stackrel{}{r})^2=\sigma ^2`$ of the jump size distribution exists. To leading order in $`k=|\stackrel{}{k}|`$, $`\widehat{\mathrm{\Phi }}(\stackrel{}{k})`$ can be expanded as

$`\widehat{\mathrm{\Phi }}(\stackrel{}{k})=1\sigma ^2k^2+𝒪(k^o),\mathrm{with}o>2.`$ (23)
* For $`\mu 2`$, the variance $`(\stackrel{}{r})^2`$ is infinite. This regime of “long jumps” leads to so-called Lévy flights. In this case, to leading order in $`k=|\stackrel{}{k}|`$, $`\widehat{\mathrm{\Phi }}(\stackrel{}{k})`$ can be expanded as

$$\widehat{\mathrm{\Phi }}(\stackrel{}{k})=1\sigma ^\mu k^\mu +𝒪(k^o),\mathrm{where}0<\mu 2,\mathrm{with}o>\mu ,$$
(24)
where $`\sigma `$ is a characteristic distance defined by

$$\sigma =\{\begin{array}{ccc}d[\mathrm{\Gamma }(1\mu )]^{1/\mu },\hfill & & 0<\mu <1,\hfill \\ \frac{d\pi }{\mu \mathrm{\Gamma }(\mu 1)\mathrm{sin}(\pi \mu /2)},\hfill & & 1<\mu <2.\hfill \end{array}$$
(25)

For a distribution $`\mathrm{\Psi }(t)`$ of waiting times of the form of a local Omori law (4) with exponent $`\theta <1`$, $`\widehat{\mathrm{\Psi }}(\beta )`$ can be expanded for small $`\beta `$ as

$$\widehat{\mathrm{\Psi }}(\beta )=1(\beta c^{})^\theta +𝒪(\beta ^\omega ),\mathrm{with}\omega 1.$$
(26)
where $`c^{}`$ is proportional to $`c`$ up to a numerical constant $`c^{}=c\left(\mathrm{\Gamma }(1\theta )\right)^{1/\theta }`$ in the case $`\theta <1`$.

Putting the leading terms of the expansions of $`\widehat{\mathrm{\Phi }}(\stackrel{}{k})`$ for small $`|\stackrel{}{k}|`$ and of $`\widehat{\mathrm{\Psi }}(\beta )`$ for small $`\beta `$ in (21) gives

$$\widehat{N}(\beta ,\stackrel{}{k})=\frac{\widehat{S}_M(\beta ,\stackrel{}{k})}{1n+n(\beta c^{})^\theta +n\sigma ^\mu k^\mu }.$$
(27)
The corresponding $`\widehat{W}(\beta ,\stackrel{}{k})`$ is obtained from (22) by

$$\widehat{W}(\beta ,\stackrel{}{k})=\widehat{S}_M(\beta ,\stackrel{}{k})\frac{(\beta )^{\theta 1}c^\theta }{1n+n(\beta c^{})^\theta +n\sigma ^\mu k^\mu }.$$
(28)
The critical regime $`n=1`$ gets rid of the constant term $`1n`$ in the denominator of (27) and (28). This case is analyzed in details below.

The regime $`n1`$ introduces a characteristic time $`t^{}`$ given by (1). In the sub-critical regime, equation (27) can be rewritten as

$$\widehat{N}(\beta ,\stackrel{}{k})=\frac{\widehat{S}_M(\beta ,\stackrel{}{k})}{(1n)}\frac{1}{1+(\beta t^{})^\theta +(kr^{})^\mu }.$$
(29)
where $`r^{}`$ is defined by

$$r^{}=\sigma \left(\frac{n}{1n}\right)^{1/\mu }.$$
(30)
For $`t<t^{}`$ and $`r<r^{}`$, the dressed propagator is given by the same expression as for the critical case and all our results below hold. For large times $`t>t^{}`$ and large distances $`r>r^{}`$, we can factorize (29) as a product of a function of time and a function of space

$$\widehat{N}(\beta ,\stackrel{}{k})\frac{\widehat{S}_M(\beta ,\stackrel{}{k})}{(1n)}\frac{1}{(1+(\beta t^{})^\theta )}\frac{1}{(1+(kr^{})^\mu )}.$$
(31)
Thus, there is no diffusion in the sub-critical regime for $`t>t^{}`$ and $`r>r^{}`$. We shall not analyze further this trivial regime $`n<1`$ and $`t>t^{}`$ and will only analyze the case $`t<t^{}`$. If there is the need, the cross-over can be calculated explicitly using (27).

In order to get the leading behavior of $`N(t,\stackrel{}{r})`$ from that of $`W(t,\stackrel{}{r})`$, we see from (21) and (22) that $`\widehat{N}(\beta ,\stackrel{}{k})=\frac{\beta }{1\widehat{\mathrm{\Psi }}(\beta )}\widehat{W}(\beta ,\stackrel{}{k})\beta ^{1\theta }c^\theta \widehat{W}(\beta ,\stackrel{}{k})`$. The inverse Laplace transform of $`1/\beta ^\theta `$ is $`1/[\mathrm{\Gamma }(\theta )t^{1\theta }]`$. Using the fact that the Laplace transform of $`df/dt`$ is $`\beta `$ times the Laplace transform of $`f(t)`$ minus $`f(0)`$, we get $`N(t,\stackrel{}{r})`$ as the derivative of a convolution

$$N(t,\stackrel{}{r})=\frac{c^\theta }{\mathrm{\Gamma }(\theta )}\frac{d}{dt}_0^t𝑑t^{}\frac{W(t^{},\stackrel{}{r})}{(tt^{})^{1\theta }}=c^\theta {}_{0}{}^{}D_{t}^{1\theta }W(t,\stackrel{}{r}).$$
(32)
In (32), we have dropped the Dirac function coming from the inverse Laplace transform of the constant term $`f(0)`$, which provides a contribution only at the origin of time $`t=0`$. Note that the operator $`\frac{1}{\mathrm{\Gamma }(\theta )}\frac{d}{dt}_0^t𝑑t^{}\frac{W(t^{},\stackrel{}{r})}{(tt^{})^{1\theta }}`$ is nothing but the so-called fractional Riemann-Liouville derivative operator of order $`1\theta `$ applied to the function $`W(t,\stackrel{}{r})`$ of time $`t`$ and is usually denoted $`{}_{0}{}^{}D_{t}^{1\theta }W(t,\stackrel{}{r})`$.

### IV.2 The standard diffusion case $`\theta >1`$ and $`\mu >2`$

The standard diffusion process is recovered for $`\theta 1`$ (for which the average waiting time is finite) and for $`\mu 2`$ (for which the variance of the jump length is finite). In this case, $`\widehat{N}(\beta ,\stackrel{}{k})=\frac{\widehat{S}_M(\beta ,\stackrel{}{k})}{\beta c^{}+\sigma ^2k^2}`$. For an impulsive source leading to $`\widehat{S}_M(\beta ,\stackrel{}{k})=`$ constant, this is the Laplace-Fourier transform of the standard diffusion propagator

$$N(t,\stackrel{}{r})\frac{1}{(Dt)^{d/2}}\mathrm{exp}[(\stackrel{}{r})^2/Dt],\mathrm{where}D=\sigma ^2/c^{},$$
(33)
where $`d`$ is here the space dimension. This solution is valid for $`|\stackrel{}{r}|/\sqrt{Dt}`$ not too large. For larger values, large deviations lead to corrections with the power law tail of the input jump distribution $`\mathrm{\Phi }(\stackrel{}{r})1/|\stackrel{}{r}|^{1+\mu }`$ defined in (5), along the lines presented for instance in Sorbook  (section 3.5). This regime is not relevant to the aftershock problem for which usually $`0<\theta <1`$.

### IV.3 Long waiting times ($`\theta <1`$) and finite variance of the jump sizes ($`\mu >2`$)

Putting the leading terms of the expansions of $`\widehat{\mathrm{\Phi }}(\stackrel{}{k})`$ (23) and of $`\widehat{\mathrm{\Psi }}(\beta )`$ (26) in (21) gives

$$\widehat{N}(\beta ,\stackrel{}{k})=\frac{1}{(\beta c^{})^\theta +(\sigma k)^2}$$
(34)
The expression (34) can be inverted with respect to the Fourier transform, and then inverted with respect to the Laplace transform using Fox functions Meltzner2 ; barkai3 . The solution for $`W(t,\stackrel{}{r})`$ in one dimension is given for instance in Meltzner2  in terms of an infinite sum

$$W(t,\stackrel{}{r})=\frac{1}{2D}\frac{1}{t^{\frac{\theta }{2}}}\underset{k=0}{\overset{\mathrm{}}{}}\frac{(1)^kz^k}{k!\mathrm{\Gamma }(1\theta (k+1)/2)}$$
(35)
where

$$z=\frac{Dt^{\theta /2}}{|\stackrel{}{r}|}$$
(36)
and $`D=\sigma /c^{\theta /2}`$.

Expression (35) and many others below involve the Gamma function of negative arguments. We recall that the Gamma function $`\mathrm{\Gamma }(u)`$ can be analytically continued to the whole complex plane, except for the simple poles $`u=0,1,2,3,\mathrm{}`$ Thus, $`\mathrm{\Gamma }(u)`$ is defined everywhere but at these poles. In order to get the expression of the Gamma function for negative arguments, one can use two formulae: $`\mathrm{\Gamma }(1u)\times \mathrm{\Gamma }(u)=\pi /\mathrm{sin}(\pi u)`$ and $`\mathrm{\Gamma }(1+u)=u\mathrm{\Gamma }(u)`$. Both these formulae are valid for all points with the possible exception of the arguments at poles $`0,1,2,\mathrm{}`$ For instance, $`\mathrm{\Gamma }(\theta )=\mathrm{\Gamma }(1\theta )/(\theta )=[\pi /\theta \mathrm{sin}(\pi \theta )]/\mathrm{\Gamma }(\theta )`$, for $`0<\theta <1`$.

Expression (35) can be rewritten as a Fox-function Mathai

$$W(t,z)=\frac{1}{2D}\frac{1}{t^{\frac{\theta }{2}}}H_{1,1}^{1,0}\left[\frac{1}{z}|\begin{array}{c}(1\theta /2,\theta /2)\hfill \\ (0,1)\hfill \end{array}\right]$$
(37)
whose asymptotic dependence for large $`z`$, obtained from a standard theorem of the Fox function (equation (1.6.3) of Mathai ),

$$W(t,z)\frac{1}{Dt^{\frac{\theta }{2}}}\frac{1}{z^{\frac{1\theta }{2\theta }}}\mathrm{exp}\left(\left(1\frac{\theta }{2}\right)\left(\frac{\theta }{2}\right)^{\frac{\theta }{2\theta }}z^{\frac{2}{2\theta }}\right)$$
(38)
is in agreement with the result of Roman and Alemany Roman  and Barkai et al. barkai3  for a space dimension $`d_f=1`$, including the dependence in the power law prefactor to the exponential. The exponential dependence $`W(t,r)\mathrm{exp}\left(\mathrm{const}(r/Dt^{\theta /2})^{\frac{2}{2\theta }}\right)`$ in (38) holds in arbitrary dimensions $`d_f`$, the only modification occurring in the prefactor whose power of $`z`$ change with the space dimension $`d_f`$ as Roman ; barkai3

$$W_{d_f}(t,z)\frac{1}{Dt^{\frac{\theta }{2}}}\frac{1}{z^{\frac{d_f(1\theta )}{2\theta }}}\mathrm{exp}\left(\left(1\frac{\theta }{2}\right)\left(\frac{\theta }{2}\right)^{\frac{\theta }{2\theta }}z^{\frac{2}{2\theta }}\right).$$
(39)
The expression of $`N(t,\stackrel{}{r})`$ can be obtained from $`W(t,\stackrel{}{r})`$ using the fractional Riemann-Liouville derivation (32) of order $`1\theta `$. Inserting expression (35) in (32) and using the expression of the fractional Riemann-Liouville derivative operator $`{}_{0}{}^{}D_{t}^{\alpha }`$ applied to an arbitrary power $`t^\mu `$, i.e., $`{}_{0}{}^{}D_{t}^{\alpha }t^\mu =\frac{\mathrm{\Gamma }(1+\mu )}{\mathrm{\Gamma }(1+\mu \alpha )}t^{\mu \alpha }`$, we obtain

$$N(t,\stackrel{}{r})=\frac{c^\theta }{2Dt^{1\frac{\theta }{2}}}\underset{k=0}{\overset{\mathrm{}}{}}\frac{(1)^kz^k}{k!\mathrm{\Gamma }((1k)\theta /2)}.$$
(40)
Expression (40) can be used to evaluate $`N(t,\stackrel{}{r})`$ for small $`z`$, but the numerical evaluation of (40) is impossible for large $`z`$. In order to obtain the asymptotic behavior of $`N(t,\stackrel{}{r})`$, expression (40) can be rewritten as a Fox-function Mathai

$$N(t,\stackrel{}{r})=\frac{c^\theta }{2Dt^{1\frac{\theta }{2}}}H_{1,1}^{1,0}\left[\frac{1}{z}|\begin{array}{c}(\theta /2,\theta /2)\hfill \\ (0,1)\hfill \end{array}\right].$$
(41)
Employing again the standard theorem of the Fox function (equation (1.6.3) of Mathai ), the asymptotic behavior of $`N(t,r)`$ for large distances $`r`$ such that $`r>Dt^{\theta /2}`$ is given by

$$N(t,r)\frac{c^\theta }{Dt^{1\frac{\theta }{2}}}\left(\frac{|\stackrel{}{r}|}{Dt^{\theta /2}}\right)^{\frac{1\theta }{2\theta }}\mathrm{exp}\left(\left(1\frac{\theta }{2}\right)\left(\frac{\theta }{2}\right)^{\frac{\theta }{2\theta }}\left(\frac{|\stackrel{}{r}|}{Dt^{\theta /2}}\right)^{\frac{2}{2\theta }}\right).$$
(42)
The exponential dependence $`N(t,r)\mathrm{exp}\left(\mathrm{const}(r/Dt^{\theta /2})^{\frac{2}{2\theta }}\right)`$ in (42) holds in arbitrary dimensions.

This expression becomes incorrect for very large distances because it would predict an exponential or slightly super-exponential decay with $`r`$. This cannot be true as the global law cannot decay faster than the local law (5). The reason for (42) to become incorrect at large distances is that the expansion of $`\widehat{N}(\beta ,\stackrel{}{k})`$ for small $`|\stackrel{}{k}|`$ (large distances) given by (34) has been truncated at the order $`k^2`$. There is however a subdominant term $`k^\mu `$ that describes the power law tail of the local law (5) and also of the global law asymptotically. A similar situation occurs in the application of the central limit theorem for sums of $`N`$ random variables with power law distributions with exponents $`\mu >2`$ Sorbook : the distribution of the sum $`S`$ is a Gaussian in its bulk for $`|S|<\sqrt{N\mathrm{ln}N}`$ and crosses over to a power law with tail exponent $`\mu `$ for larger $`S`$. In a similar way, the cross-over of $`N(t,r)`$ to the asymptotic local power law (5) can be recovered by an analysis including the subleading correction $`k^\mu `$ to the expansion (34).

Expression (40) shows that the global rate of seismicity cannot be factorized as a product of a distribution of times and a distribution of distances. This space-time coupling implies that the seismic activity diffuses with time, and that the decay of the rate of aftershocks depends on the distance from the first mainshock. This coupling of space and time stems from the cascade of aftershocks, from the primary aftershocks to the secondary aftershocks to the tertiary aftershocks and so on.

Figure 4 presents the decay of the seismic activity $`N(r,t)`$ obtained using expression (40) for small $`z`$ and expression (42) for large $`z`$, as a function of the time from the mainshock and as a function of the distances $`r`$. Close to the mainshock epicenter, expression (40) predicts that the global seismicity rate decays with time as the renormalized Omori law

$$N(t,0)\frac{1}{t^{1\theta /2}}.$$
(43)
The same decay is found at any fixed point $`\stackrel{}{r}`$ for times $`t>(|\stackrel{}{r}|/D)^{2/\theta }`$. At all times, the same decay $`1/t^{1\theta /2}`$ is also obtained by measuring the aftershock seismicity in a local box at a distance from the main shock origin increasing with time as $`rt^{\frac{\theta }{2}}`$ (this is nothing but putting $`z=`$ constant in (40)). At large distances $`r>Dt^{\theta /2}`$, the global decay law is different from a power-law decay. Figure 4 shows that the rate of aftershocks presents a truncation at early times, which increases as the distance $`r`$ increases. At large times, the rate of aftershocks recovers the $`1/t^{1\theta /2}`$ power-law decay (43). We stress that a fit of the global law $`N(r,t)`$ over the whole time interval by an Omori law would yield an apparent exponent $`p<1\theta /2`$ that decreases with $`r`$.

Integrating (40) over the whole one-dimensional space, we recover the global Omori law

$$N(t)=𝑑rN(t,r)\frac{1}{t^{1\theta }}$$
(44)
found in SS ; HS1 . Thus, we have found an additional source of variability of the exponent $`p`$ of the Omori law: if measured over the whole catalog, we should measure $`p=1\theta `$ in the critical regime $`n=1`$ while $`p=1\theta /2`$ is slightly larger when measured in certain time- and space-windows, as described above. Thus, in this regime, pruning of catalogs may lead to continuous change from the value $`1\theta `$ to $`1\theta /2`$. In addition, as we have mentioned, the cross-over in time may lead to still smaller apparent exponents, thus enhancing the impression of variability of the exponent $`p`$. In reality, this range of $`p`$-values are seen to result from the complex spatio-temporal organization of the aftershock seismicity of the ETAS model. These results should lead us to be cautious when analyzing real catalogs with respect to the conditions and regimes under which the analysis is performed.

There is another observable that characterizes how an aftershock sequence invades space as a function of time. Expression (40) indeed predicts a sub-diffusion process quantified by

$$|\stackrel{}{r}|^2t^{2H},$$
(45)
with $`H=\theta /2`$ since the natural variable is $`z`$ given by (36). Indeed, expression (40) tells us that, up to a global rescaling function of time, the rate of aftershocks is identical for a fixed value of $`z`$. Thus, any aftershock structure diffuses according to (45).

This prediction is checked in Figure 5 by numerical simulations. 1000 synthetic catalogs have been generated with $`\mu =3`$, $`\theta =0.2`$ and $`n=1`$. The average distance between the first mainshock and its aftershocks as a function of the time from the mainshock has been averaged over these 1000 simulations. The theoretical diffusion exponent is $`H=\theta /2=0.1`$, in good agreement with the asymptotic behavior observed in the numerical simulation. In practice, in order to minimize the effect of fluctuations and optimize the speed of convergence, we estimate numerically $`\mathrm{exp}[\mathrm{ln}|\stackrel{}{r}|]`$ which is also expected to scale as $`\mathrm{exp}[\mathrm{ln}|\stackrel{}{r}|]t^{\theta /2}`$ due to the simple scaling form of (41).

This problem has also been solved exactly in barkai1  in the context of the so-called fractional Fokker-Planck equation, which amounts to replace the distribution $`\mathrm{\Phi }(\stackrel{}{r})`$ of jumps (5) by a Gaussian function. This fractional Fokker-Planck equation allows one to introduce the possibility of bias or drift in the CTRW and therefore in the aftershock sequence.

### IV.4 Exponential waiting time distribution and long jump size Lévy distribution ($`\mu <2`$)

This case with an exponential distribution

$$\mathrm{\Psi }(t)=\lambda e^{\lambda t}$$
(46)
of waiting times with a Lévy distribution $`\mathrm{\Phi }(\stackrel{}{r})=L_\mu (|\stackrel{}{r}|)`$ of jump sizes with tail exponent $`\mu <2`$ has been investigated by Budde et al. Budde . One finds

$$|\stackrel{}{r}|^2^{1/2}t^{1/\mu },$$
(47)
corresponding to a superdiffusion regime with Hurst exponent $`H=1/\mu >1/2`$. The full distribution function $`W(t,\stackrel{}{r})`$ corresponding to the critical regime $`n=1`$ is known for $`\lambda t>>1`$:

$$W(t,\stackrel{}{r})\frac{1}{(\lambda t)^{1/\mu }}L_\mu \left(\frac{|\stackrel{}{r}|}{(\lambda t)^{1/\mu }}\right).$$
(48)
The corresponding $`N(t,\stackrel{}{r})`$ is obtained from (20). The Laplace transform of the exponential distribution (46) is $`\widehat{\mathrm{\Psi }}(\beta )=\lambda /(\beta +\lambda )`$. We thus get

$$\widehat{N}(\beta ,\stackrel{}{k})=\left(\beta +\lambda \right)\widehat{W}(\beta ,\stackrel{}{k}),$$
(49)
and thus

$$N(t,\stackrel{}{r})=\frac{W(t,\stackrel{}{r})}{t}+\lambda W(t,\stackrel{}{r}).$$
(50)
Expression (50) together with (48) predicts a diffusion law $`rt^H`$ with $`H=1/\mu `$ which is in good agreement with our simulations. At large times $`|\stackrel{}{r}|(\lambda t)^{1/\mu }`$, $`N(t,\stackrel{}{r})\lambda W(t,\stackrel{}{r})1/t^{1/\mu }`$, given an apparent local Omori exponent $`\theta =11/\mu `$. This offers a new mechanism for generating Omori law for aftershocks from purely exponential local relaxation but with a heavy distribution of jump sizes. This power-law decay should be observed only at a fixed distance $`r`$ or over a limited domain from the mainshock in the regime of large times.

Integrating over the whole space, $`𝑑\stackrel{}{r}W(t,\stackrel{}{r})=1`$ which gives $`N(t)=\delta (t)+\lambda `$ equal to a constant seismic rate. This results from an initial mainshock at $`t=0`$ leading to the cascade of aftershocks adjusting delicately to this constant rate for the critical value $`n=1`$ of the branching parameter. In the sub-critical regime $`n<1`$, the Omori law integrated over space gives instead $`N(t)\mathrm{exp}[(1n)\lambda t]`$, showing that the characteristic decay time $`1/(1n)\lambda `$ of the dressed Omori propagator $`N(t)`$ becomes much larger (much longer memory) that the decay time $`1/\lambda `$ of the bare Omori propagator.

For $`\mu >2`$, we recover the standard diffusion corresponding to $`\theta >1`$ and $`\mu >2`$ discussed in section IV.2.

### IV.5 Long waiting times ($`\theta <1`$) and long jump sizes (Lévy flight regime for $`\mu 2`$)

Putting the leading terms of the expansions of $`\widehat{\mathrm{\Phi }}(\stackrel{}{k})`$ and of $`\widehat{\mathrm{\Psi }}(\beta )`$ in (21) gives

$$\widehat{N}(\beta ,\stackrel{}{k})=\widehat{S}_M(\beta ,\stackrel{}{k})\frac{1}{(\beta c^{})^\theta +(\sigma k)^\mu }.$$
(51)
The corresponding $`\widehat{W}(\beta ,\stackrel{}{k})`$ is given by

$$\widehat{W}(\beta ,\stackrel{}{k})=\widehat{S}_M(\beta ,\stackrel{}{k})\frac{(\beta )^{\theta 1}c^\theta }{(\beta c^{})^\theta +(\sigma k)^\mu }.$$
(52)
Equation (52) has been studied extensively in the context of the CTRW model as a long wavelength $`|\stackrel{}{k}|0`$ and long time $`\beta 0`$ approximation to investigate the long time behavior of the CTRW. Kotulski Kot1  has developed a rigorous approach, based on limit theorems, to classify the asymptotic behaviors of different type of CTRWs and justifies the approximation (52) for the long time behavior. Barkai barkai2  has studied the quality of the long wavelength $`|\stackrel{}{k}|0`$ and long time $`\beta 0`$ approximation (52) by solving the exact CTRW problem for the case when the waiting time distribution $`\mathrm{\Psi }(t)`$ is a one-sided stable Lévy law of index $`\theta `$ with the same tail as (4) and the distribution $`\mathrm{\Phi }(\stackrel{}{r})`$ of jumps is a symmetric stable Lévy of index $`\mu `$ with the same tail as (5). Their Laplace and Fourier transforms, that appear in the denominator of (22), are respectively $`\widehat{\mathrm{\Psi }}(\beta )=\mathrm{exp}[\beta ^\theta ]`$ and $`\widehat{\mathrm{\Phi }}(\stackrel{}{k})=\mathrm{exp}[|\stackrel{}{k}|^\mu /2]`$. Note that the long wavelength $`|\stackrel{}{k}|0`$ and long time $`\beta 0`$ approximation gives $`1\mathrm{exp}[(c^{}\beta )^\theta ]\mathrm{exp}[|\sigma \stackrel{}{k}|^\mu ]=(c^{}\beta )^\theta +|\sigma \stackrel{}{k}|^\mu `$, which recovers (51). By comparing the exact solution of (21) for $`\mathrm{\Psi }(t)`$ and $`\mathrm{\Phi }(\stackrel{}{r})`$ of the above Lévy form with that of the long wavelength $`|\stackrel{}{k}|0`$ and long time $`\beta 0`$ approximation (52), Barkai barkai2  finds that certain solutions of (52) diverge on the origin, a behavior not found for the corresponding solutions of (21). In addition, certain solutions of the full equation (21) converge only very slowly for $`\mu <1`$ to the solutions of the long-time approximation (52). These results validate our use of the asymptotic long time behavior with respect to the scaling laws but provide a note of caution if one needs more precise non-asymptotic information. In this case, such information can be obtained by a suitable analysis of the full equation (21).

Using power counting, expression (52) predicts a diffusion process (45) with exponent

$$H=\frac{\theta }{\mu }.$$
(53)
This prediction is checked by numerical simulation of the ETAS model in the critical regime $`n=1`$, with $`\theta =0.2`$, $`\mu =0.9`$, shown in figure 6. The average distance between the first mainshock and its aftershocks as a function of the time from the mainshock indeed increases according to (45) with an exponent $`H`$ in very good agreement with the prediction $`H=\theta /\mu =0.22`$. As the form of the denominator in (52) is independent of the space dimension, the prediction (53) is valid in any space dimension.

The natural variable for the expansions given below allowing to compute $`N(t,\stackrel{}{r})`$ is

$$z=\frac{Dt^{\theta /\mu }}{|\stackrel{}{r}|},$$
(54)
where $`D=\sigma /c^{\theta /\mu }`$ and $`c^{}=c\left(\mathrm{\Gamma }(1\theta )\right)^{1/\theta }`$.

#### IV.5.1 $`z`$-expansion of the solution

$`W(t,\stackrel{}{r})`$ can be obtained as the following sum (equation (5.10) of Saichev )

$$W(t,\stackrel{}{r})=\frac{1}{\pi |\stackrel{}{r}|}\underset{m=0}{\overset{+\mathrm{}}{}}(1)^mz^{m\mu }\frac{\mathrm{\Gamma }(m\mu +1)}{\mathrm{\Gamma }(m\theta +1)}\mathrm{cos}\left[\frac{\pi }{2}(m\mu +1)\right].$$
(55)
Applying (32) to (55) term by term in the sum, we get

$$N(t,\stackrel{}{r})=\frac{c^\theta }{D\pi t^{1\theta +\theta /\mu }}\underset{m=0}{\overset{+\mathrm{}}{}}(1)^mz^{1+m\mu }\frac{\mathrm{\Gamma }(m\mu +1)}{\mathrm{\Gamma }((m+1)\theta )}\mathrm{cos}\left[\frac{\pi }{2}(m\mu +1)\right],$$
(56)
The asymptotics

$$\frac{\mathrm{\Gamma }(m\mu +\mu +1)\mathrm{\Gamma }(m\theta +1)}{\mathrm{\Gamma }(m\theta +\theta +1)\mathrm{\Gamma }(m\mu +1)}\frac{\mathrm{\Gamma }(m\mu +\mu +1)\mathrm{\Gamma }((m+1)\theta )}{\mathrm{\Gamma }((m+2)\theta )\mathrm{\Gamma }(m\mu +1)}m^{\mu \theta }$$
(57)
show that the series (55) and (56) exist only for $`\mu <\theta `$. It can be shown that these series exist for all $`z`$ in this case. This series converges very slowly for large $`z`$ but the Padé summation method BeOr  can be used to improve the convergence of (56) in the case $`\mu <\theta `$, and can also be used to estimate (56) in the case $`\mu >\theta `$ for which the series diverges.

The space integral $`𝑑rN(t,r)`$ over the whole one-dimensional volume $`V`$, with $`N(t,r)`$ given by (56), recovers the global Omori law

$$_V𝑑rN(t,r)\frac{1}{t^{1\theta }}.$$
(58)
Note the non-trivial phenomenon in which the superposition of all aftershock activities transforms the local Omori law or “bare propagator” (4) $`\mathrm{\Psi }(t)\frac{1}{t^{1+\theta }}`$ into the global Omori law or “dressed propagator” $`\frac{1}{t^{1\theta }}`$. This effects was predicted in SS ; HS1  in the version of the ETAS model without space dependence. These results are consistent with the claim of section II.4 according to which all results reported previously for the version of the ETAS model without space dependence hold also for the version of the space-dependent ETAS model studied here, when averaging over the whole space.

The asymptotic behavior for $`|\stackrel{}{r}|Dt^{\frac{\theta }{\mu }}`$ (i.e., $`z1`$) and $`\mu <\theta `$ is obtained by keeping only the first non-zero term ($`m=1`$) in (56) which is convergent for all $`z`$ in the case $`\mu <\theta `$

$$N(t,\stackrel{}{r})=\frac{\mathrm{sin}\left(\frac{\pi \mu }{2}\right)}{\sigma c^{}\pi }\frac{\mathrm{\Gamma }(1+\mu )}{\mathrm{\Gamma }(2\theta )}\left(\frac{c^{}}{t}\right)^{12\theta }\left(\frac{\sigma }{|\stackrel{}{r}|}\right)^{1+\mu },\mathrm{for}|\stackrel{}{r}|Dt^{\frac{\theta }{\mu }}.$$
(59)
At fixed large $`|\stackrel{}{r}|`$ and for $`t<|\stackrel{}{r}/D|^{\frac{\mu }{\theta }}`$, this predicts a local Omori law with exponent $`p=12\theta `$.

#### IV.5.2 $`1/z`$-expansion of the solution

We use the theory of Fox functions Mathai  to obtain $`N(t,\stackrel{}{r})`$ as an infinite series in $`1/z`$. For this, we first rewrite expression (56) as a Fox function Mathai

$$N(t,\stackrel{}{r})=\frac{c^\theta }{D\mu \pi t^{1\theta +\theta /\mu }}R\left(H_{2,2}^{1,2}\left[ze^{i\pi /2}|\begin{array}{c}(1/\mu ,1/\mu ),(1,1)\hfill \\ (1/\mu ,1/\mu ),(\theta /\mu \theta +1,\theta /\mu )\hfill \end{array}\right]\right),$$
(60)
where $`R(z)`$ indicates the real part of $`z`$.

The $`1/z`$ expansion of $`N(t,\stackrel{}{r})`$ can be obtained using the dual expansion of the Fox function (60) (expression (3.7.2) of Mathai )

$$N(t,\stackrel{}{r})=\frac{c^\theta }{D\pi \mu t^{1\theta +\theta /\mu }}\underset{m=0}{\overset{+\mathrm{}}{}}(1)^m[\mu z^{1\mu m\mu }\frac{\mathrm{\Gamma }(1(m+1)\mu )\mathrm{sin}((m+1)\mu \pi /2)}{\mathrm{\Gamma }(m\theta )}$$

$$+\frac{z^m}{m!}\frac{\pi \mathrm{cos}(m\pi /2)}{\mathrm{sin}((m+1)\pi /\mu )\mathrm{\Gamma }(\theta (m+1)\theta /\mu )}].$$
(61)
This expansion exists only for $`\mu >\theta `$ (conditions of page 71 below eq. (3.7.2) of Mathai ). This is easily checked by the behavior of an asymptotics similar to (57). Note that the series (61) is not defined in the special case $`\mu =1`$ due to the presence of the ill-defined ratio $`\mathrm{\Gamma }(0)/\mathrm{\Gamma }(0)`$ and a different approach is required, such as the integral representation of $`W(t,\stackrel{}{r})`$ developed in Saichev . The global Omori law obtained by integrating over the whole space (61) is again $`N(t)1/t^{1\theta }`$ as expected from the analysis of the ETAS model without space dependence HS1 .

Keeping only the largest term of (61) for large $`z`$, we obtain the asymptotic behavior for small distances $`r<Dt^{\theta /\mu }`$

$`N(t,r)`$ $``$ $`{\displaystyle \frac{\mathrm{\Gamma }(12\mu )\mathrm{sin}(\pi \mu )\mathrm{sin}(\pi \theta )}{c^{}\sigma \pi ^2}}{\displaystyle \frac{\mathrm{\Gamma }(1+\theta )}{(r/\sigma )^{12\mu }}}{\displaystyle \frac{1}{(t/c^{})^{1+\theta }}}\mathrm{for}\mu <0.5`$
$`N(t,r)`$ $``$ $`{\displaystyle \frac{c^\theta }{c^{}\sigma \mu \mathrm{\Gamma }(\theta \theta /\mu )\mathrm{sin}(\pi /\mu )}}{\displaystyle \frac{1}{(t/c^{})^{1\theta +\theta /\mu }}}\mathrm{for}0.5<\mu <2.`$ (62)
Note that for $`r<Dt^{\theta /\mu }`$ and $`0.5<\mu <2`$, the leading behavior of $`N(t,r)`$ is independent of $`r`$.

Equation (62) thus predicts an apparent exponent

$`p`$ $`=`$ $`1+\theta \mathrm{for}\mu <0.5`$
$`p`$ $`=`$ $`1\theta +\theta /\mu \mathrm{for}0.5<\mu <2`$ (63)
for small distances $`r<Dt^{\theta /\mu }`$. This prediction is valid only in the case $`\mu >\theta `$ for which the series (61) is convergent. However, the same asymptotic results are also obtained by different methods in the case $`\mu <\theta `$, for instance expression (63) is recovered for all $`\mu <2`$ using the integral representation of Saichev  \[A. Saichev, private communication\]. The numerical evaluation of (56), which converges for $`\mu <\theta `$, also recovers the asymptotic results (62). The two regimes $`\mu <0.5`$ and $`0.5<\mu <2`$ are illustrated in Figures 7 and 8 respectively. The seismicity rate $`N(t,\stackrel{}{r})`$ is evaluated from expression (56) for small $`z`$ and from expression (61) for large $`z`$.

We also performed numerical simulations of the ETAS and CTRW models and the results are in good agreement with expression (56) and (61) for $`N(\stackrel{}{r},t)`$ for $`tc`$ and $`rd`$. For very small times $`tc`$, or for very small distances $`rd`$, expressions (56) and (61) are not valid because they are based on a long wavelength $`|\stackrel{}{k}|0`$ and long time $`\beta 0`$ approximation. Numerical simulations of the ETAS model in the case $`\theta =0.2`$ and $`\mu =0.9`$ are presented in Figure 9, and are in good agreement with the analytical solutions (56) and (61) shown in Figure 8 for the same parameters, except from the truncation of $`N(t,r)`$ for times $`tc`$ and distances $`rd`$ that are not reproduced by the analytical solution.

### IV.6 A simple non-separable joint distribution of waiting times and jump sizes: coupled spatial diffusion and long waiting time distribution

Consider the choice for $`\varphi _{m_i}(tt_i,\stackrel{}{r}\stackrel{}{r}_i)`$ replacing (2) by

$$\varphi _{m_i}(tt_i,\stackrel{}{r}\stackrel{}{r}_i)=\rho (m_i)\mathrm{\Psi }(tt_i)\mathrm{\Phi }(|\stackrel{}{r}\stackrel{}{r}_i|/\sqrt{Dt}),$$
(64)
where $`\rho (m_i)`$ and $`\mathrm{\Psi }(t)`$ are again given by (3) and (4) while (5) is changed into

$$\mathrm{\Phi }(|\stackrel{}{r}\stackrel{}{r}_i|/\sqrt{Dt})=\frac{1}{\sqrt{2Dt}}\mathrm{exp}\left(|\stackrel{}{r}\stackrel{}{r}_i|^2/Dt\right).$$
(65)
The spatial diffusion of seismic activity is now coupled to the waiting time distribution. Expression (65) captures the effect that, in order for aftershocks to spread over large distances by the underlying physical process, they need time. In fact, returning to the discussion in the introduction on the various proposed mechanisms for aftershocks, expression (65) embodies a microscopic diffusion process.

In this case, (21) must be replaced by

$$\widehat{N}(\beta ,\stackrel{}{k})=\frac{\widehat{S}_M(\beta ,\stackrel{}{k})}{1n\widehat{\varphi }(\beta ,\stackrel{}{k})},$$
(66)
where $`\widehat{\varphi }(\beta ,\stackrel{}{k})`$ is the Laplace-Fourier transform of the product $`\mathrm{\Psi }(t)\mathrm{\Phi }(|\stackrel{}{r}|/\sqrt{Dt})`$. For large times and long distances for which the first terms in the expansion in $`\beta `$ and $`k`$ are sufficient, and for $`n=1`$, we obtain

$$\widehat{\varphi }(\beta ,\stackrel{}{k})\frac{\widehat{S}_M(\beta ,\stackrel{}{k})}{(\beta +Dk^2)^\theta }.$$
(67)
The inverse Laplace-Fourier transform of (66) is

$$N(t,\stackrel{}{r})\frac{1}{t^{1\theta }}\frac{1}{\sqrt{2\pi Dt}}\mathrm{exp}\left(|\stackrel{}{r}|^2/Dt\right).$$
(68)
As expected, expression (68) recovers the dressed Omori propagator in the case of absence of space dependence HS1 . At finite $`r`$ and long times, the dressed Omori law also decay as $`1/t^{1\theta }`$. The diffusion of aftershocks is normal with the standard diffusion exponent $`H=1/2`$.

## V New questions on aftershocks derived from the CTRW analogy

We list a series of comments and questions suggested from the analogy between the ETAS model and the CTRW model. In particular, we discuss the possibility of defining new observables for earthquake aftershocks, that could be worthwhile to investigate in future empirical studies of earthquake aftershocks.

### V.1 Recurrence of aftershock activity in the proximity of the main shock

A quantity often investigated in studies of random walks is the probability $`W(t,\stackrel{}{0})`$ to find the random walker at its starting point (the origin) at time $`t`$. In the earthquake framework, this is the seismic aftershock rate close to the main shock.

### V.2 First-passage times

The first passage time of a random walk is the first arrival time of the random walk at a given point $`\stackrel{}{r}`$. In the earthquake context, this translates into the study of the waiting time for a given region to have its own first aftershock after the main shock occurs. The distribution of such first passage waiting times gives the distribution of times with no nearby seismic activity. See for instance barkai1  in the case of a power law distribution of waiting times and Gaussian distribution of jump sizes. Margolin and Berkowitz Margolin  give the distribution of first-passage times in the case where the jump distribution is narrow and the waiting distribution is long-tailed $`1/t^{1+\theta }`$. They analyze the three different regimes $`\theta <1`$, $`1<\theta <2`$ and $`\theta 2`$.

### V.3 Occupation time of seismic activity

Weiss and Calabrese Weiss2  have studied the total amount of time spent by a lattice CTRW on a subset of points. In the seismic language, this amounts to study the probability distribution of the durations of aftershock sequences that are localized in a specific subset of the space. In other words, how probable are aftershock sequences that are found only within a given spatial subset over a certain duration?

### V.4 Transience and recurrence of seismic activity

Another question that has been studied in some details in the CTRW framework is whether random walks are transient or recurrent. A transient random walk visits any point $`\stackrel{}{r}`$ at most a finite number of times before escaping to infinity. For earthquakes, the transient regime corresponds to the activation of at most a finite number of aftershocks in any given point $`\stackrel{}{r}`$. In contrast, a recurrent random walk may return a growing number of times to all or a subset of points at time increases. In the aftershock language, this means that these points will have a never-ending (decaying) aftershock activity. We stress here the difference between the global Omori law giving a never-ending power law decay of the aftershock activity (in the sub-critical regime $`n<1`$) and its spatial dependence which must exhibit important variations. In particular, in the recurrent regime, an Omori law can be documented by counting aftershocks in those limited regions of space which are activated again and again.

### V.5 Probability for the cumulative number of aftershocks

Let us define a basic quantity in the CTRW formalism, namely the probability $`\chi _m(t)`$ to make exactly $`m`$ steps up to time $`t`$. In the earthquake context, $`\chi _m(t)`$ is the probability to have exactly $`m`$ aftershocks after the main shock. In the case in which the spatial transition probability $`\mathrm{\Phi }(\stackrel{}{r})`$ between different positions is independent of the waiting times (corresponding to factorizing $`\varphi _{m_i}(tt_i,\stackrel{}{r}\stackrel{}{r}_i)`$ as in (2)), the probability density $`W(t,\stackrel{}{r})`$ to find the walker at position $`\stackrel{}{r}`$ at time $`t`$ can be written

$$W(t,\stackrel{}{r})=\underset{m=0}{\overset{+\mathrm{}}{}}W_m(\stackrel{}{r})\chi _m(t),$$
(69)
where $`W_m(\stackrel{}{r})`$ is the probability to reach $`\stackrel{}{r}`$ from $`\stackrel{}{0}`$ in $`m`$ steps. In the earthquake context, $`W_m(\stackrel{}{r})`$ is the probability that there has been exactly $`m`$ events in the time interval $`[0,t]`$ and that the last one occurred at $`\stackrel{}{r}`$. Equation (69) states that the CTRW is a random process subordinated to simple random walks described by $`W_m(\stackrel{}{r})`$ under the operational time given by the $`\chi _m(t)`$ distribution Sokolov1 ; Sokolov2 .

### V.6 Random walk models with birth and death and background seismicity from localized sources

Bender et al. Bender  have studied models of random walks in which walkers are born in proportion to the population at one specific site (for instance the origin) with probability $`a1`$ (with $`a>1`$) and die at all other sites with probability $`1n`$ (with $`n1`$). In the earthquake context, this consists in assuming that the aftershock activity is fed by a localized region in space, which is itself activated by the aftershocks returning to this region, furthering the overall activity. This may be considered to describe the seismic activity close to a plate boundary, in which the plate boundary is the constant self-consistent source of a seismic activity which may spread over a significant region away from the boundary. The excursion of the random walkers quantify the spread of the seismic activity away from the main fault structure. The rate of death of the walkers correspond exactly to the distance $`1n`$ from the critical value $`n=1`$. Bender et al. Bender  find a phase diagram in the $`(a1,1n)`$ parameter space in which a boundary separates two possible asymptotic regimes:

1. for small $`a1`$ and large $`1n`$, the seismic activity at the origin and everywhere eventually dies off;
2. for large $`a1`$ and small $`1n`$, the average seismic activity at the origin approaches a positive constant at long times. In this regime, there is a transition as $`a1`$ is decreased or as $`1n`$ is increased, between a case where the global seismic activity outside the origin goes to zero and a case where it diverges at long times. On the boundary between these two regimes in the $`(a1,1n)`$ parameter space, the distribution of seismic activity approaches a steady state at long times. There is a critical point (for space dimensions different from $`2`$) at a certain value $`(a_c1,1n_c)`$, for which the long-time seismic activity away from the source is given by $`(aa_c)^\nu `$ where $`\nu `$ is a critical exponent equal to $`2`$ in three dimensions.

Note that the results of Bender  are obtained for random walks on a lattice. This can easily be converted into a CTRW by the fact that a CTRW is nothing by a process subordinated to discrete random walks under the operational time defined by the process $`\{t_i\}`$ of the time of just arrival to a given site, as given by (69).

## VI Discussion

Using the analogy between the ETAS model and the CTRW model established here, we have derived the relation between the average distance between aftershocks and the mainshock as a function of the time from the mainshock, and the joint probability distribution of the times and locations of aftershocks.

We have assumed that each earthquake triggers aftershocks at a distance $`r`$ and time $`t`$ according to the bare propagator $`\varphi (r,t)`$, which can be factorized as $`\mathrm{\Psi }(t)\mathrm{\Phi }(r)`$. This means that the distribution $`\mathrm{\Phi }(r)`$ of the distances between an event and its direct aftershocks is decoupled from the distribution $`\mathrm{\Psi }(t)`$ of waiting time. Hence, the direct aftershocks triggered by a single mainshock do not diffuse in space with time. Notwithstanding this decoupling in space and time of the bare propagator $`\varphi (r,t)`$, we have shown that the global law or dressed propagator $`N(t,\stackrel{}{r})`$ defined as the global rate of events at time $`t`$ and at position $`\stackrel{}{r}`$, cannot be factorized into two distributions of waiting times and space jumps. This joint distribution of waiting times and positions of the whole sequence of aftershocks cascading from a mainshock is different from the product of the bare time and space propagators.

The mean distance between the mainshock and its aftershocks, including secondary aftershocks, increases with the time from the mainshock, due to the cascade process of aftershocks triggering aftershocks triggering aftershocks, and so on. In the critical case $`n=1`$, this diffusion takes the form of a power-law relation $`Rt^H`$ of the average distance $`R`$ between aftershocks and the mainshock, as a function of the time $`t`$ from the mainshock. If the local Omori law is characterized by an exponent $`0<\theta <1`$, and if the space jumps follow a power law $`\mathrm{\Phi }(r)1/(r+d)^{1+\mu }`$, the diffusion exponent is given by $`H=\theta /\mu `$ in the case $`\mu <2`$ and $`H=\theta /2`$ in the case $`\mu >2`$. Depending on the $`\theta `$ and $`\mu `$ values, we can thus observe either sub-diffusion ($`H<1/2`$) or super-diffusion ($`H>1/2`$), as summarized in Figure 10. In the sub-critical ($`n<1`$) and super-critical ($`n>1`$) regimes, this relation is still valid up to the characteristic time $`t^{}`$ given by (1) and for distances smaller than $`r^{}Dt^H`$ given by (30). For $`t>t^{}`$ and $`r>r^{}`$ in the sub-critical regime, the global distributions of times and distances between the mainshock and its aftershocks are decoupled and there is therefore no diffusion. In the super-critical regime, the aftershock rate increases exponentially for $`t>t^{}`$ and the aftershocks diffuses more rapidly than before $`t^{}`$.

In the critical regime, the cascade of secondary aftershocks introduces a variation of the apparent Omori exponent as a function of the distance from the mainshock. The asymptotic values of the Omori exponent in the different regimes are summarized in Table 2. In the regime $`\mu <2`$, we observe a transition from an Omori law decay with an exponent $`p=12\theta `$ at early times $`t^Hr/D`$ to a larger exponent at large times. This provides another mechanism to explain the observed variability of the Omori exponent. In the regime $`\mu >2`$, a power-law decay of the seismicity with time is observed only at large times $`t^Hr/D`$. At early times, or at large distances $`rDt^H`$, the seismicity rate is very small, because the seismicity as not yet diffused up to the distance $`r`$.

We should emphasize that our theoretical analysis of aftershock diffusion predicts the behavior of the ensemble average of aftershock sequences. Individual sequences may depart from this ensemble average, especially for sequences with few earthquakes and limited durations. For long sequences (20,000 events say), we have verified that the exponent $`H`$ measured on individual sequences does not deviate from the ensemble average value by more than about 20%. As already discussed, the impact of fluctuations become however more effective as the parameter $`\alpha `$ increases above $`b/2`$.

The diffusion of the seismicity also renormalizes the spatial distribution of the seismicity, which is very different from the local distribution $`\mathrm{\Phi }(r)`$ of distances between a triggering event and its direct aftershocks. In the regime $`\mu >2`$, the global seismicity rate $`N(t,\stackrel{}{r})`$ decays exponentially with the distance from the mainshock, whereas the local distribution of distances $`\mathrm{\Phi }(r)`$ is a power-law distribution. In the regime $`\mu <2`$, the local law $`\mathrm{\Phi }(r)r^{1\mu }`$ is recovered at large distances, but a slower decay for $`0.5<\mu <2`$ or a constant rate for $`\mu <0.5`$ is observed at small distances $`rDt^H`$. These predictions on the decrease of the Omori exponent with $`r`$ have not yet been observed in earthquake catalogs, but an expansion of the aftershock zone has been reported in many studies Mogi ; Imoto ; Chatelain ; Tajima1 ; Tajima2 ; Wesson ; Ouchi ; Noir ; Jacques . However, very few studies have quantified the diffusion law. Noir et al. Noir  show that the earthquake Dobi sequence (central Afar, August 1989) composed of 22 $`M>4.6`$ earthquakes presented a migration that was in agreement with a diffusion process due to fluid transfer in the crust, characterized by a normal diffusion process with exponent $`H=0.5`$. Tajima and Kanamori Tajima1 ; Tajima2  studied several aftershock sequences in subduction zone and observed a much slower logarithmic diffusion, which is compatible with a low exponent $`H`$ close to $`0.1`$. In some cases, the aftershock sequence displays no expansion with time. For instance, Shaw Shaw  studied several aftershock sequences in California and concluded that the distribution of distances between the mainshock and its aftershocks is independent of time. This can be explained by the fact that the 0mori exponent measured in Shaw  is very close to $`1`$, thus $`\theta `$ is very small and our prediction is that the exponent $`H`$ should be very small.

In fact, the ETAS model predicts that diffusion should be observed only for aftershock sequences with a measured Omori exponent $`p`$ significantly smaller than $`1`$, which can only occur according to our model when the bare Omori propagator with exponent $`1+\theta `$ is renormalized into the dressed propagator with global exponent $`1\theta `$. We have shown that this renormalization of the exponent only occurs at times less than $`t^{}`$, while for longer times in the sub-critical regime $`n<1`$ the dressed Omori propagator recovers the value of the bare exponent $`1+\theta >1`$ (see figure 2). Therefore, identifying an empirical observation of $`p<1`$ with our prediction $`p=1\theta `$ indicates that the aftershock sequence falls in the “good” time window $`t<t^{}`$ in which the renormalization operates. We have also shown that the dressed propagator gives a diffusion only for $`t<t^{}`$. We can thus conclude that, according to the ETAS model, the observation of an empirical Omori exponent larger than $`1`$ is indicative of the large time $`t>t^{}`$ behavior in the sub-critical regime $`n<1`$, for which there is no diffusion. This provides a possible explanation for why many sequences studied by Tajima1 ; Tajima2 ; Shaw  do not show a diffusion of the aftershock epicenters. Reciprocally, a prerequisite for observing diffusion in a given aftershock sequence is that the empirical $`p`$-value be less than $`1`$ in order to qualify the regime $`t<t^{}`$.

An alternative model has been discussed by Dieterich Diete  who showed that the spatial variability of the stress induced by a mainshock, coupled with a rate and state friction law, results in an expansion of the aftershock zone with time. This expansion does not take the form of a diffusion law as observed in the ETAS model, the relation between the characteristic size of the aftershock zone does not grow as a power law of the time from the mainshock (equation (22) and Figure 6 of Diete ).

Marsan et al. Marsan1 ; Marsan2  and Marsan and Bean Marsan3  studied several catalogs at different scales, from the scale of a deep mine to the world-wide seismicity, and observed that the average distance between two earthquakes increases as a power-law of the time between them, with an exponent often close to 0.2, indicative of a sub-diffusion process. They interpreted their results as a mechanism of stress diffusion, that may be due to fluid transfer with heterogeneous permeability leading to sub-diffusion. Their analysis is quite different from those used in other studies, because they consider all pairs of events, without distinction between aftershocks and mainshocks. This analysis can however lead to spurious diffusion, and in some cases this method does not detect diffusion in synthetic data set with genuine diffusion. We have tested their analysis on a synthetic catalog generated by superposing a background seismicity with uniform spatial and temporal distribution, and 10 mainshocks with poissonian distribution in time and space, and with a power-law distribution of energies. Each of these mainshocks generates only direct aftershocks, without secondary cascades of aftershocks, and the number of aftershocks increases exponentially with the magnitude of the mainshock. This way, we generate a synthetic catalog without any physical process of diffusion, and which includes all the other well-established characteristics of real seismicity: clustering in space and time superposed to a seismicity background. Applying the analysis of Marsan1 ; Marsan2 ; Marsan3  to this synthetic data set leads to an apparent diffusion process with a well-defined exponent $`H=0.5`$. However, this apparent diffusion does not reflect a genuine diffusion but simply describes the crossover from the characteristic size of an aftershock zone at early times to the larger average distance between uncorrelated events at large times. In plain words, the apparent power law $`Rt^H`$ is nothing but a cross-over and is not real. Furthermore, applying this analysis to a synthetic catalog generated using the ETAS model, without seismicity background, and with a theoretical diffusion exponent $`H=0.2`$, the method yields $`H=0.01`$ if we use all the events of the catalog. If we select only events up to a maximum distance $`r_{max}`$ to apply the same procedure as in Marsan1 ; Marsan2 ; Marsan3 , we obtain larger values of $`H`$ which are more in agreement with the theoretical exponent $`H=0.2`$ but with large fluctuations that are function of $`r_{max}`$. Therefore, it is probable that the diffusion reported in Marsan1 ; Marsan2 ; Marsan3  is not real and results from a cross-over between two characteristic scales of the spatial earthquake distribution. It may be attributed to the analyzing methodology which mixes up uncorrelated events. We are thus reluctant to compare the results of Marsan et al. Marsan1 ; Marsan2 ; Marsan3  with the predictions obtained with the ETAS model.

One can similarly question the results on anomalous diffusion of seismicity obtained by Sotolongo-Costa et al. Sotolongo , who considered 7500 micro-earthquakes recorded by a local spanish network from 1985 to 1995. They interpret the sequence of earthquakes as a random walk process, in which the walker jumps from an earthquake epicenter to the next in sequential order. The time between two successive events is seen as a waiting time between two jumps and the distance between these events is taken to correspond to the jump size. Since the distributions of time intervals and of distances between successive earthquakes are both heavy-tailed (approximately power laws), their model is a CTRW. We cannot stress enough that their CTRW model of seismicity has nothing to do with our results on the mapping of the ETAS model onto a CTRW. Their procedure is ad-hoc and their results depend obviously strongly on the space domain of the analysis since distant earthquakes that are completely unrelated can be almost simultaneous! We also stress that our mapping of the ETAS model onto the CTRW model does not correspond to identifying an earthquake sequence as a single realization of a CTRW, as assumed arbitrarily by Sotolongo-Costa et al. Sotolongo .

Our predictions obtained here are thus difficult to test on seismicity data, due to the small number of events available and the restricted time periods and distance ranges, and because the seismicity background can strongly affect the results. New methods should hence be developed to investigate if there is a real physical process of diffusion in seismic activity and to compare the observations of real seismicity with the quantitative predictions of the ETAS model. Preliminary study of aftershock sequences in California leads to the conclusion that most aftershock sequences are characterized by an Omori exponent $`p>1`$, indicative of the sub-critical regime with $`t>t^{}`$. As expected from our predictions in this regime, we do not observe an expansion of the aftershock zone. However, a few sequences give a value $`p<1`$ and also exhibit an increase of the average distance between the mainshock and its aftershocks consistent with our predictions. A detailed report of this analysis will be reported elsewhere.

## VII Conclusion

We have studied analytically and numerically the ETAS (epidemic-type aftershock) model, which is a simple stochastic process modeling seismicity, based on the two best-established empirical laws for earthquakes, the power law decay of seismicity after an earthquake and a power law distribution of earthquake energies. This model assumes that each earthquake can trigger aftershocks, with a rate increasing with its magnitude. In this model, the seismicity rate is the result of the whole cascade of direct and secondary aftershocks.

We have first established an exact correspondence between the ETAS model and the CTRW (continuous-time random walk) model. We have then used this analogy to derive the joint probability of times and distances of the seismicity following a large earthquake and we have characterized the different regimes of diffusion.

We have shown that the diffusion of the seismicity should be observed only for times $`t<t^{}`$, where $`t^{}`$ is a characteristic time depending on the model parameters, corresponding to an observed Omori exponent smaller than one. Most aftershock sequences have an observed Omori exponent larger than one, corresponding to the subcritical regime of the ETAS model, for which there is no diffusion. The diffusion of the seismicity produces a decrease of the Omori exponent as a function of the distance from the mainshock, the decay of aftershocks being faster close to the mainshock than at large distances. The spatial distribution of seismicity is also renormalized by the cascade process, so that the observed distribution of distances between the mainshock and its aftershocks can be fundamentally different from the bare propagator $`\mathrm{\Phi }(r)`$ which gives the distribution of the distances between triggered and triggering earthquakes. We have also noted that the ETAS model generates apparent but realistical fractal spatial patterns.

Assuming that the distances between triggering and triggered events are independent of the time between them, this model generates a diffusion of the whole sequence of aftershocks with the time from the mainshock, which is induced by the cascade of aftershocks triggering aftershocks, and so on. Our results thus provides a simple explanation of the diffusion of aftershock sequences reported by several studies, which was often interpreted as a mechanism of anomalous stress diffusion. We see that no such “anomalous stress diffusion” is needed and our theory provides a parsimonious account of aftershock diffusion resulting from the minimum physical ingredients of the ETAS model. As Einstein once said, “A theory is more impressive the greater the simplicity of its premises, the more different the kinds of things it relates and the more extended its range of applicability.”

###### Acknowledgements.

We are very grateful to B. Berkowitz, S. Gluzman, J.-R. Grasso, Y. Klafter, L. Margerin, A. Saichev and G. Zaslavsky for useful suggestions and discussions.