{'fraction_non_alphanumeric': 0.04107753705815279, 'fraction_numerical': 0.014709236031927024, 'mean_word_length': 4.5011761016151794, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Universal mean moment rate profiles of earthquake ruptures

Earthquake phenomenology exhibits a number of power law distributions including the Gutenberg-Richter frequency-size statistics and the Omori law for aftershock decay rates. In search for a basic model that renders correct predictions on long spatio-temporal scales, we discuss results associated with a heterogeneous fault with long range stress-transfer interactions. To better understand earthquake dynamics we focus on faults with Gutenberg-Richter like earthquake statistics and develop two universal scaling functions as a stronger test of the theory against observations than mere scaling exponents that have large error bars. Universal shape profiles contain crucial information on the underlying dynamics in a variety of systems. As in magnetic systems, we find that our analysis for earthquakes provides a good overall agreement between theory and observations, but with a potential discrepancy in one particular universal scaling function for moment-rates. The results reveal interesting connections between the physics of vastly different systems with avalanche noise.


INTRODUCTION
Earthquake phenomenology is characterized by several power law distributions. The most famous of these are the frequency-size distributions (i.e. histograms) of regional and global earthquakes [1,6], and the modified Omori law for the aftershock decay rate around large rupture zones [2,3]. Using the seismic moment M 0 for the earthquake size, the frequency-size distributions (or moment histogram) has the form (e.g., [4]) where M 0 ∼ i ∆u i ∆A i with ∆u i and ∆A i being the local slip and rupture area during an earthquake, respectively. A related, more commonly used form in terms of earthquake magnitude M is log(n(M )) = a − bM (2) where n(M 0 )dM 0 = n(M )dM , the constant a characterizes the overall rate of activity in a region, and the "bvalue" gives the relative rates of events in different magnitude ranges. Using the observed moment-magnitude scaling relation M ∼ 2/3log(M 0 ) for large earthquakes [4,5,6], the exponent β of (1) is related to the b-value of (2) as b = 1.5β. The modified Omori law for aftershock decay rates is: where N is the cumulative number of aftershocks, t is the time after the mainshock, and K, c, and p are empirical constants. The exponents in (1) and (3) are stable for data collected over large space-time domains, with some clear deviations from global averages related to faulting type and regional properties [6]. For example, the b-values of strike-slip, thrust, and normal earthquakes with depth ≤ 50 km in the global Harvard catalogue are about 0.75, 0.85, and 1.05, respectively [7]. As another example, regions with high heat flow often have short aftershock sequences with relatively large exponent (e.g., p > 1. 25), while regions with low heat flow have long aftershock sequences with low exponent (e.g., p < 0.9) [2,3]. The association of earthquake statistics with power law relations like Eq. (1) and Eq. (3) led some to suggest that earthquake dynamics is associated with an underlying critical point [8,9,10,11,12]. However, power law distributions can be generated by many other mechanisms [15,16] and it is important to develop criteria that can provide stronger evidence for or against the association of earthquakes with criticality.
Recently, enough data have been collected to extract statistics of earthquakes on individual fault zones occupying long (order 100 km) and narrow (order 10 km) regions of space. Wesnousky and collaborators [13,14] found that the frequency-size statistics of earthquakes on highly irregular fault zones, with many offsets and branches, as the San Jacinto fault in California, are also described by the Gutenberg-Richter power law relation up to the largest events. However, relatively regular fault zones (presumably generated progressively with increasing accumulated slip over time), such as the San Andreas fault in California, display power-law frequency-size statistics only for small events. These occur in the time intervals between roughly quasi-periodic earthquakes of a much larger "characteristic" size that is related to large-scale dimensions of the fault zone [4,13,17,18]. (If the ratio of the mean divided by the standard deviation of the dis-tribution of time intervals between characteristic earthquakes is larger than 1, the distribution is referred to as quasi-periodic [18]). Earthquakes of intermediate magnitude are typically not observed on these faults (other than, perhaps, during aftershock sequences). The corresponding frequency size statistics are called the "characteristic earthquake" distribution [4,13].
Previously these two types of behavior on individual fault zones have been modeled as statistics close-to and far-from an underlying critical point [8,19], using a model for a strike-slip fault that incorporates long-range interactions and strong heterogeneities [17,18]. The different dynamic regimes were associated with a competition between failure-promoting effects of elastic stresstransfer or dynamic weakening, and the opposing effect of strength inhomogeneities in the fault structure. Fisher et al. [8] found that near the critical point the frequencysize statistics follow a power law distribution (with a cutoff at large magnitude), with the same scaling exponent of observed data for strike-slip faults (i.e., a b-value of 0.75). A similar form of frequency-size statistics and predicted b-value were obtained also for a critical parameter value in a stochastic branching model [20].
To provide an improved understanding of earthquake dynamics that can suggest additional observables, we focus on faults with Gutenberg-Richter like earthquake statistics (i.e. near-critical behavior) and develop two universal scaling functions associated with mean moment-rate time profiles at either fixed total moment or fixed total earthquake duration. Universal scaling functions (or shape profiles) give important information on the underlying dynamics, and may be found in solar flares in astrophysics [21], price fluctuations in financial markets [22], Barkhausen noise in magnets [23], and, as shown here, also in earthquake phenomenology. If the behavior of fault zones with earthquakes following the Gutenberg-Richter statistics is indeed critical, then the shapes of these functions should be as universal as the exponent β in Eq. (1). Comparing the scaling functions in our earthquake model to observations constitutes a much stronger test of the theory than merely comparing a discrete, finite set of critical exponents.
In this study we compute the scaling functions for both model predictions and observational data and compare the results. In section of the paper we review the model. In sections to , we extend the model, while also introducing an extended phase diagram for the model dynamics, and the scaling behavior on long length scales. In section we introduce the universal scaling functions and their scaling forms, and in section we extract the functions from both simulation and observational data. Finally, in section we discuss the results and emerging new questions.

EARTHQUAKE MODEL
The model we use was developed originally by Ben-Zion and Rice [17,18], who suggested that a narrow irregular strike-slip fault zone of horizontal length L and vertical depth W may be represented by an array of N ∼ LW cells in a two dimensional plane, with constitutive parameters that vary from cell to cell to model the disorder (offsets etc.) of the fault zone structure (FIG. 1). The cells represent brittle patches on the interface between two tectonic blocks that move with slow transverse velocity v in the x direction at a great distance from the fault. The interaction between cells during slip events is governed by 3-D elasticity and falls off with a distance r from the failure zone as 1 r 3 . These interactions are sufficiently long range that scaling in mean field theory (where the interaction range is set to infinity) becomes exact, up to logarithmic corrections, in the physical fault dimension (d = 2) [8,17,18].
In mean field theory, the local stress τ i on a given cell i is [17]: where u i is the total offset of cell in the horizontal x direction,ū = j u j /N is the average displacement, J/N is the mean-field elastic coupling strength between cells, and K L ∼ 1/ √ N is the loading stiffness [19] of the tectonic blocks that move far away from the fault with relative velocity v. Initially the stresses τ i are randomly distributed with τ a,i ≤ τ i ≤ τ s,i , where τ s,i is a fixed local static failure threshold stress and τ a,i is the fixed local arrest stress. The distributions of static failure stresses and arrest stresses represent the heterogeneity or geometrical disorder in the fault system. The differences between the failure and arrest stresses give the local distribution of stress drops during brittle failures; the earthquake dynamics depend only on the stress drop distribution (no stochasticity). In addition, the scaling behavior of the system is not sensitive to the exact form of the distributions as long as they are bounded and τ a,i < τ s,i . We choose a compact distribution for τ a,i , such that p(τ a,i ) = 3(W 2 − 4τ 2 a,i )/(2W 3 ) for −W/2 ≤ τ a,i ≤ W/2 and 0 outside of these bounds. Also, we look at the low disorder limit where W ≪ τ s,i , where we choose τ s,i = 1, so all cells will fail at this point.
The fault is stuck while the stress on each cell is increased uniformly as dτ i /dt = K L v as a result of the external loading which is increased adiabatically (that is, we take the limit v → 0). When the stress on a cell reaches its failure threshold τ s,i , the cell slips by the amount: This stress drop is uniformly redistributed to all other cells (employed in the mean field approximation) by an amount: where c ≡ J/(K L + J) is the conservation parameter which gives the fraction of the stress drop of a cell that is retained in the system after it slips [19]. The resulting stress increase on the other cells can cause some of them to slip as well, leading to an avalanche of cell slips, or an earthquake.

DYNAMICAL WEAKENING AND STRENGTHENING
The model includes dynamic weakening effects during the failure process [17,18]: after an initial slip in an earthquake, the strength of a failed cell is reduced to a dynamical value: with 0 ≤ ǫ ≤ 1 parameterizing the relative importance of dynamical weakening effects in the system. This weakening represents the transition from static friction to dynamic friction during the rupture and the strength of a failed cell remains at the dynamic value throughout the remainder of the earthquake. In the time intervals between earthquakes all failure thresholds heal back to their static value τ s,i . Fisher et al. [8] found that at exactly ǫ = 0 the model produces a power law distribution of earthquake moments M 0 following equation (1), cutoff by the finite fault size, with an analytical exponent β = 1/2 (FIG. 2). This corresponds to a b-value of 0.75, close to that associated with observed earthquakes on strike-slip faults [7]. The power law scaling of the frequency-size statistics and other variables [8] indicates that the model with ǫ = 0 operates at a critical point. In contrast, for a finite weakening ǫ > 0 the model produces the characteristic earthquake distribution, with power law statistics for the small events up to a cutoff moment that scales like: and quasi-periodically recurring large characteristic events that scale with the fault size (M 0 ∼ (LW ) 3/2 ). The model can be expanded further to include dynamic strengthening represented by ǫ < 0. Ben-Zion and Sammis [24] summarized multidisciplinary observations which indicate that brittle failure of rock has an initial transient phase associated with strengthening, distributed deformation, and creation of new structures. Detailed frictional studies also show an initial strengthening phase associated with the creation of a new population of asperity contacts [4,25]. In our model (FIG. 1) we associate ǫ < 0 with regions off the main fault segments that are in an early deformation stage. To capture basic aspects of brittle deformation on such regions in the three-dimensional volume around the main fault (FIG. 1), we change the model as follows: when any cell i slips during an earthquake, and thereby reduces its stress by ∆τ i ≡ τ f,i − τ a,i , the failure stress τ f,j of every cell j = 1, ..., N is strengthened by an amount |ǫ|∆τ i /N . Once the earthquake is complete, the failure stress of each cell is slowly lowered back to its original value. This represents in a simple way the brittle deformation that occurs during an earthquake in the off-fault regions, which are first in a strengthening regime, compared to the main fault, and then have a weakening process. The events that are triggered as the failure stresses are lowered in the weakening period are referred to as aftershocks. The occurrence of aftershocks in this version of the model for off-fault regions is in agreement with the observation that a large fraction of observed aftershocks typically occur in off fault regions [26]. For this version of the model with ǫ < 0, both the primary earthquakes (i.e., main shocks) and the triggered aftershocks are distributed according to the Gutenberg-Richter distribution, up to a cutoff moment scaling as 1/ǫ 2 . Assuming that the increased failure stress thresholds τ f,i are slowly lowered with time as log(t) towards their earlier static values τ s,i , and that the stresses are distributed over a wide range of values, we show analytically in Appendix that the temporal decay of aftershock rates at long times is proportional to 1/t, as in the modified Omori law of Eq. (3) with p = 1 [2,3,4].
Remarkably, the long length scale behavior of this model can be shown to be the same as the behavior of the model given in Eq. (5) with an added "antiferroelastic" term (−|ǫ|Jū): In Eq. (10) every time a cell fails, it slips by an amount ∆u i that leads to stress loading of the other cells, lessened by |ǫ|J∆u i /N compared to our original model (Eq. (5)).
On the other hand, in the global strengthening model (described above) when a cell slips the failure stresses of all cells are strengthened by |ǫ|J∆u i /N . On long length scales the global strengthening of the failure stress has equivalent effects on the earthquake statistics as the dissipation of the redistributed stress, up to corrections of order O(1/N ), so the scaling behavior for large events of both models are the same. Moreover, Eq. (10) can be rewritten as: We can now absorb |ǫ| by defining J ′ = J(1 − |ǫ|) and K ′ L = K L + J|ǫ|. Rewriting Eq. (11) with the new definitions, and dropping the |ǫ| contribution in [K ′ L − J|ǫ|]vt since v → 0, we find: Therefore we recover Eq. (5) with J → J ′ and K L → K ′ L . This amounts to changing the stress conservation parameter c (from reference [19]). For Eq. (12): where K L → 0 since we are concerned with the adiabatic limit. We also know (from reference [19]) that the cutoff S cf for the Gutenberg-Richter distribution scales as Thus, from Eq. (13) we find that the cutoff for Eq. (10) will scale as ∼ 1/|ǫ| 2 .

MAPPING TO SINGLE INTERFACE MAGNET MODEL
The mean field version of the single interface magnet model with infinite range antiferromagnetic interactions is given by [27,28]: where h i (t) is the position of the domain wall, H(t) is the external driving field, k is the coefficient of the antiferromagnetic term, and η i (h) is the pinning field. In the paper by Fisher et al. [8] it has been shown that the scaling behavior on long length scales resulting from Eq. (10), without the −|ǫ|Jū term, is same as that of Eq. (14) without the antiferromagnetic term −kh. Furthermore, upon inspection we see the following correspondence between the single interface magnet model (Eq. (14)), and the mean field earthquake model (Eq. (10)): In other words, the coefficient of the antiferromagnetic term k plays the same role in the magnet model (Eq. (14)), as the coefficient of strengthening |ǫ|J does in the earthquake model (Eq. (10)).

PHASE DIAGRAM
The regimes with various statistics produced by the model are summarized by the phase diagram given in FIG. 2. The range ǫ > 0 corresponds to "mature" localized faults with a weakening rheology and characteristic earthquake statistics. The value ǫ = 0 corresponds to "immature" strongly inhomogeneous fault zones with Gutenberg-Richter statistics. Finally, the range ǫ < 0 corresponds to the fracture network away from the main fault, characterized by strengthening due to the creation of new structures and associated emerging aftershocks. It may be surprising that the discussed simple model can capture many of the essential general features of earthquake statistics (or other systems with avalanches, such as driven magnetic domain walls). This can be understood through the renormalization group [29,30], a powerful mathematical tool to coarse grain a system and extract its effective behavior on long space-time scales. Many microscopic details of a system are averaged out under coarse graining, and universal aspects of the behavior on long scales depend only on a few basic properties such as symmetries, dimensions, range of interactions, weakening/strengthening, etc. When a model correctly captures those basic features, the results provide proper predictions for statistics, critical exponents, and universal scaling functions near the critical point. Consequently, many models that are in the same universality class lead to the same statistics and exponents [8,19,29,30]. The universal scaling functions around the critical point, discussed in the next section, provide additional information that can be used to distinguish between different models and universality classes.

MOMENT RATE SHAPES
In this section we focus on fault zones with Gutenberg-Richter power law statistics, modeled by systems at or close to the ǫ = 0 critical point. Recent analysis allowed researchers to obtain the moment rate dm 0 (t)/dt, which gives the slip on a fault per unit time during the propagation of earthquake rupture, for hundreds of large seismic events recorded on global networks [31,36]. The moment rates are derived from inversions of teleseismically recorded seismograms on a global seismic network [38]. Motivated by works on statistical physics of magnetic systems discussed in the previous chapter (see also [23,30]), we are interested in studying the event-averaged moment rate time profile (FIG. 3) for earthquakes with given total moment M 0 , denoted with dm 0 (t|M 0 )/dt and the event-averaged moment rate time profile (FIG.  4) for earthquakes with given duration T , denoted with dm 0 (t|T )/dt . Here m 0 (t|T ) is the (cumulative) moment at time t of the propagating earthquake of total duration T , and m 0 (t|M 0 ) is the cumulative moment at time t of the earthquake of total moment M 0 . Theoretical analysis of phase diagrams similar to that shown in FIG. 2 implies that near the critical point there should be, in addition to scaling exponents, also universal scaling functions (up to a rescaling of the ordinate and abscissa) [23]. In our model the two scalable functions of interest, dm 0 (t|M 0 )/dt and dm 0 (t|T )/dt , obey respectively the following scaling relations [32,33]: and We determined these scaling functions from corresponding results for magnets, using the fact that our mean field version of the Ben-Zion and Rice model of Eq. (5) [8,19] is in the same universality class (i.e. has the same universal behavior on long length scales) as the above mentioned model for domain wall motion in magnets.

EXPONENTS AND DATA COLLAPSES
We compare the observation results with our model and find remarkable agreement in most cases. The frequency-moment distribution, D(M 0 ) ∼ M −1−β 0 of the observed data [31] has (inset of FIG. 2) three decades of scaling and an exponent of β = 1/2 ± 0.05, in close agreement with the model near ǫ = 0. The deviation from a power law distribution at the low moment range is associated with the reduced resolution of the observational network for small events. In mean field theory the universal scaling function f (x) in Eq. (16) is of the exact form [33]:f mf (x) = Axe −Bx 2 /2 with non-universal constants A = B = 1. In FIG. 3 we present a collapse of the observational data of dm 0 (t|M 0 )/dt for four different values of M 0 to obtain the corresponding function f exp (x) for observations with x = t/M 1/2 0 . The observational curves not only collapse, and are therefore universal, the mean field exponent 1/2 in the scaling variable x is in excellent agreement with observations. We fit the functional form f mf (x) with A = 4 and B = 4.9 to the collapse of the observed data; f exp (x) deviates from the f mf (x) for small values of the ordinate.
In mean field theory, the function g(x) of Eq. (17) is of the symmetric form: g mf (x) = Ax(1 − x), where A is a non-universal constant. In FIG. 4 we collapse observational data for dm 0 (t|T )/dt for three values of T to obtain the function g exp (x) with x = t/T . Again we find that the curve collapses well, even though only small data sets were available, and the exponent of 1 obtained from the collapse is in excellent agreement with mean field theory.
We also plot the mean field scaling function g mf (x) with A = 80. The results show that while the scaling exponents agree, there are notable differences between the observational function g exp (x) and the mean field function g mf (x), especially for small values of the ordinate. We checked that finite size effects do not play a role in g mf (x) (or in f mf (x) for that matter). Also, we find that the mean skewness coefficient for the g exp (x) curves is 0.878 and the mean standard error of skewness is 0.705: since twice the standard error is greater than the absolute value of the skewness, the asymmetry is not statistically significant. Therefore more observational work with a larger data set is required to verify the moment rate shape asymmetry and clarify its origin. An asymme-try may result from a rupture process that begins with a failure of a large asperity, from finite fault size effects if the rupture process slows down once the rupture has traversed the fault in one direction, or from contributions of early aftershocks.

DISCUSSION
The employed earthquake model [17,18] was shown in the past to have a critical point at ǫ = 0 and additional dynamic regimes for ǫ > 0 [8,19] compatible with observed frequency-size statistics of earthquakes on individual fault zones [13,14]. We have generalized the theoretical analysis to include a strengthening regime ǫ < 0 with aftershocks, and derived universal scaling functions around the critical point ǫ = 0. The results provide new tools for data analysis that may be used to obtain an improved understanding of earthquake dynamics. The analysis indicates that near ǫ = 0 the model is in the same universality class as a recent model for domain wall motion in magnets and we can match the phase diagram regions ǫ < 0 and ǫ > 0 to those of the corresponding magnet model [8,28,30,32]. In other words, the two systems are marked by the exact same universal scaling exponents, universal scaling functions, and similar phase diagrams. The model predictions for frequency-size statistics and moment-rates of earthquakes near ǫ = 0 are overall in close agreement with observational data of relatively large earthquakes recorded on the global seismic network [31]. However, the observed mean moment-rate of earthquakes with a given duration T apparently increases with time more rapidly than it drops off, contrary to the corresponding symmetric model function. The potential asymmetry in the observed data may result from a rupture process that begins with a failure of a large asperity. This is compatible with observations that hypocenter locations tend to be located close to an area on a fault that produces large moment release, (e.g. [34]). As we explain below, an analogy to magnetic systems suggests that an asymmetry could potentially also stem from momentary initial threshold strengthening associated with the creation of a new population of asperity contacts upon local failure and followup aftershocks. Further study will be necessary to clarify this issue.
Theoretical analysis of such a potential asymmetric rupture process requires corrections to the mean field earthquake model results. Recently it has been shown that the corresponding magnetic domain-wall model [23,30] predicts well the critical scaling exponents for Barkhausen noise experiments in magnets. Significantly, the experimental scaling function for magnetization avalanches or Barkhausen "pulses" [23,30], that is the analogue of the moment-rate time profile for fixed earthquake duration of Eq. (17), shows the same type of asymmetry that is apparently observed for earthquakes (FIG. 4). It has been suggested that the asymmetry in the function is due to eddy currents in the magnet [35]. Eddy currents have a similar effect in magnets as transient threshold strengthening would have on earthquakes. In this paper we have shown how long-term threshold strengthening (on time scales longer than individual earthquakes) leads to aftershocks, which effectively represent an asymmetry on time scales longer than individual earthquakes. This raises the possibility that the origin of this asymmetry may be similar in both magnets and earthquakes: in both cases it may be due to a transient force (due to eddy currents or threshold strengthening respectively) that counteracts the propagation of an event and thus leads to asymmetric event profiles that taper off more slowly than they began. Our study shows that there are important theoretical and observational connections between processes in earthquake and magnetic systems.
A related study of earthquake moment rate shapes was done by Houston [36] who used a data set similar to the data of Bilek [31] used here. In order to compare the average moment rate profiles of these earthquakes, the profiles were rescaled by two methods: moment scaling and duration scaling. Moment scaling rescales the time axis of the profiles using assumptions about cracklike scaling of the events (different from our mean-field result), and then rescales the height of the profiles so that the area (or moment) under all profiles are the same. The duration scaling of [36] rescales the time axis such that all duration-scaled profiles end at the same reference duration, and then rescales the moment rates so that all scaled time profiles have the same area underneath. Fig. 5 shows average moment rate profiles that were obtained by Houston from moment scaling (top) and duration scaling (bottom) for data from several subduction zones. The top plots correspond to dm 0 (t|M 0 )/dt collapsed, and the bottom plots correspond to dm 0 (t|T )/dt collapsed. We see from the top part in Fig. 5 that overall the moment rate shapes seem to agree rather well with our predicted mean field shapes of Fig. 3. Likewise the bottom plots of Fig. 5 appear to agree quite well the mean field curve of Fig. 4, except for an additional slight asymmetry, with positive skewness, similar to our result from observations.
In [36] the skewness of the duration scaled moment rate profiles, i.e. the skewness of the dm 0 (t|T )/dt collapses, is calculated and found to range from 0.12 to 0.36 depending on the depth of the earthquakes (bigger skewness for greater depth) and the method used to extract moment rates from seismograms. Since the statistical error of the skewness is not given in [36], we cannot determine if the above skewness values are statistically significant. It is interesting, however, that just like in the case of magnets, and in our analysis of the observational data, the skewness is slightly positive, indicating that on average earthquakes tend to grow faster than they die down. The skewness values quoted in [36] are smaller than the 0.878 value obtained here, though one of them (the 0.36 value for earthquakes at greater depth) is within errorbars.
While overall the shapes of the moment rate profiles look rather similar for both Houston's analysis and our mean field theory results, there are significant differences in the analysis that can be summarized as follows: (1) The exponents used in [36] to create the momentscaled average profile dm 0 (t|M 0 )/dt differ from the exponents we used. In both cases the x axis is divided by M αt 0 and the y axis is divided by M 1−αt 0 in order to ensure normalization of the curves to a fixed reference moment. In [36] the exponent value is α t = 1/3, while in this paper we use the mean field value α t = 1/2 (see Fig.  3). The value α t = 1/3 in [36] is associated with cracklike scaling; for cracks the moment M 0 scales with the rupture area A as M 0 ∼ A 3/2 [4]. Since the rupture area scales with the diameter L of the crack as A ∼ L 2 and the duration T of the earthquake scales with its diameter as T ∼ L one arrives at the proposed scaling T ∼ M 1/3 0 . In [36] it is mentioned, however, that the data would also be compatible with rescaling the time (or x) axis with an exponent of α t ≃ 0.41, i.e. that T ∼ M 0.41 0 for the data. This scaling is much closer to the mean field prediction T ∼ M 1/2 0 , i.e. α t = 1/2, that is used in this work. Our exponent α t = 1/2 is derived from the mean field prediction [8] that the duration T of earthquakes in the critical power law region of Fig. 2, scale as T ∼ M 1/2 0 . In [8] it is shown analytically that the moment M 0 of earthquakes in the critical regime scales with the rupture area as M 0 ∼ A. Since the rupture area scales, as before, with the earthquake diameter as A ∼ L 2 , and the earthquake duration scales as T ∼ L, one obtains for the mean field prediction T ∼ M 1/2 0 for critical earthquakes. Note that for earthquakes which are so large that their horizontal diameter L is much larger than the vertical width of the fault, one expects the moment to scale as M 0 ∼ L. Since again the duration scales as T ∼ L, one obtains T ∼ M 0 , i.e. α t = 1. The data used here and in [36] apparently does not fall within this scaling regime.
(2) Our theoretical analysis leads to predictions of both scaling exponents and entire scaling functions, which can be compared to data. In contrast, in [36] the shape is obtained only empirically by averaging rescaled data, and it is dependent on the chosen scaling exponents. Our procedure of comparing the theoretical results to data is designed to provide simultaneously, through a scaling collapse that minimizes deviations of individual averaged curves from the collapsed one, the scaling exponent α t and the entire scaling function. On the other hand, in [36] the value of α t is obtained from a separate analysis of duration versus moment of observed data, and the shape is then obtained empirically, by averaging rescaled data.
(3) For our simulated moment rates, the exponents used in [36] would not lead to a good collapse, i.e. the deviation between various rescaled curves would be larger for α t = 1/3 compared to α t = 1/2.
Unfortunately, since the available observational data have large error bars, we cannot determine at present unequivocally, which set of exponents works better. More moment rate data, especially for many small earthquakes would allow us to reduce statistical errorbars and obtain a more precise comparison of the scaling functions in theory and experiment. The smaller earthquakes are also more likely to fall into the critical scaling regime of our mean field prediction. However, moment rate data are more difficult to obtain for small earthquakes, as they require greater spacial resolution. We hope that this study will motivate more observational work and analysis to answer these questions.

DETERMINATION OF OMORI LAW OF AFTERSHOCK DECAY
We derive the Omori law for aftershock decay from our mean field earthquake model in the case of dynamic strengthening. After a primary earthquake, the static failure stress, τ s , is increased by an amount ǫ s /N , where s is the mean event size. Given in terms of moment, the mean avalanche size in mean field theory is simply given by [19]: where 1/s 3/2 ∼ D(s), is the scaling behavior of the mean field event size distribution. Now since s ∼ 1/ǫ, the failure stress, τ s , is then shifted by an amount 1/N . In order to produce aftershocks, we assume that the new larger failure stress, which we call τ 0 f , is lowered slowly, logarithmically in time until it returns to its original value, τ s ≃ τ 0 f − 1/N . More precisely, using t = 0 for the time at the end of the previous earthquake, τ 0 f decays as follows: We define η(t) as the amount by which τ 0 f has been lowered at time t: Eq. (19) is only accurate to first order since lowering the threshold triggers on average ∼ η(t) aftershocks (assuming nonsingular stress distributions). Since aftershocks occur quickly compared to the weakening time scale, the aftershocks that are triggered due to lowering the threshold will cause τ f (t) to shift up by: η(t) × ǫ s /N ∼ η(t)/N . Thus the expression for τ f (t) becomes: and Eq. (20) gives: If we assume that the number of aftershocks triggered by time t is A(t) ∼ η(t), we obtain the modified Omori law (Eq. (3)) for aftershock decay [2,3,4]: For large t and large N we have the original Omori law: We are grateful to Susan Bilek for giving us the observational data. We thank James P. Sethna 3: A collapse of averaged earthquake pulse shapes, dm0(t|M0)/dt , with the size of the moment M0 in Newton meters within 10% of each size given in the legend respectively. In order to obtain each collapsed moment rate shape, five to ten earthquakes were averaged for each value of M0. The collapse was obtained using the mean field scaling relation [8]: dm0(t|M0)/dt /M  4: A collapse of averaged earthquake pulse shapes, dm0(t|M0)/dt with a duration of T (seconds) within 10% (given in legend), is shown. The collapse was obtained using the mean field scaling relation [32]: dm0(t|T )/dt ∼ g(t/T ) . In order to obtain each collapsed pulse shape, two to ten earthquakes were averaged for each value of T . In our mean field theory the universal scaling function is g mf (x) = Ax(1 − x) with x = t/T . We plot this functional form (bold curve) with A = 80. Note the apparent asymmetry to the left in the observed data while the theoretical curve is symmetric around its maximum. Inset: The raw data and the averaged data (before collapsed). Central America (11) Aleutians (6) Kuriles (11) Vanuatu (6) Tonga-Kermadec (4) FIG. 5: Rescaled moment rate versus time profiles of seismic events determined with moment scaling (top) and duration scaling (bottom) obtained from [36] (see text). The data shown are obtained from several subduction zones in the indicated geographical locations and the numbers of earthquakes used for each region are given in parentheses.