Mehdi Zakroum

Abdellah Houmz

Mounir Ghogho

Ghita Mezzour
email: ghita.mezzour@uir.ac.maabdelkader.lahmadi

Abdelkader Lahmadi

Jérôme Franc ¸ois § Mohammed El Koutbi

Exploratory Data Analysis of a Network Telescope Traffic and Prediction of Port Probing Rates

Keywords: Cyber Intelligence, Cyber Security, Network Telescope, Darknet, Probing Patterns, Transition Graphs, Prediction of Probing Rates, Non-stationary Autoregressive Model, Non-stationary Vector Autoregressive Model, Machine Learning

Understanding the properties exhibited by large scale network probing traffic would improve cyber threat intelligence. In addition, the prediction of probing rates is a key feature for security practitioners in their endeavors for making better operational decisions and for enhancing their defense strategy skills. In this work, we study different aspects of the traffic captured by a /20 network telescope. First, we perform an exploratory data analysis of the collected probing activities. The investigation includes probing rates at the port level, services interesting top network probers and the distribution of probing rates by geolocation. Second, we extract the network probers exploration patterns. We model these behaviors using transition graphs decorated with probabilities of switching from a port to another. Finally, we assess the capacity of Non-stationary Autoregressive and Vector Autoregressive models in predicting port probing rates as a first step towards using more robust models for better forecasting performance.

I. INTRODUCTION

New cyber threat vectors and vulnerabilities are constantly emerging with the evolution of technology. Attackers commonly scan networks to find vulnerable devices which can be used for malicious intents. One of the major attacks happened in 2016 is the Dyn DDoS attack. The attackers used botnets of vulnerable devices as a primary source of their DDoS traffic generation, making leading internet platforms unavailable for a large number of users.

Improving our knowledge on scan activities will help to prevent cyber attacks through early detection, and in general, to enhance security policies. Many causes can trigger scan campaigns such as vulnerability disclosure, worm spread and zero days. Generally, such malicious traffic is hidden by a large amount of legitimate traffic, making it complex to be identified by internet service providers and network security operators to protect target users.

A passive approach for identifying network probing activities are network telescopes, also known as darknets. A network telescope is a sensor logging the traffic received by a set of passive unallocated network addresses. Therefore, the traffic received by the network telescope is considered suspicious, requiring thus to be examined.

To collect such traffic, we use a network telescope hosted at INRIA Nancy-Grand Est consisting of nearly 4096 IPv4 addresses. By analyzing the collected data, we aim to answer the following questions:

• What are the most targeted services? What are the services targeted by the top network probers? • How network probers are exploring the target network?

How to model these probing activities? • Can we predict probing rates of the targeted services? The remainder of this paper is structured as follows. The next section provides a review of the related work to this study. Section III presents an exploratory analysis of the darknet traffic. In section IV, we identify the attackers probing patterns. Finally, in section V, we explore the capacity of Nonstationary Autoregressive and Vector Autoregressive models to forecast the probing rates at the port level.

II. RELATED WORK

Reconnaissance is the first phase in the cyber kill chain, where the attacker scans the target infrastructure looking for vulnerabilities. A more generic approach for finding vulnerable devices consists of scanning the whole IPv4 address space, including network telescopes. Many studies leverage the traffic captured by the latter to study different aspects of probing activities. Durumeric et al. [START_REF] Durumeric | An Internet-Wide View of Internet-Wide Scanning[END_REF] studied the traffic acquired by a large network telescope consisting of 5.5M IP addresses. The study includes the origin of scans, the targeted services by network probers and the effect of vulnerability publication on probing activities. Bou-harb et al. [START_REF] Bou-Harb | A Statistical Approach for Fingerprinting Probing Activities[END_REF] used a probabilistic and statistical approach to identify the origin of the probing activities: whether they are generated by scanning tools or by worms and botnets. They also studied whether probing activities are random or they exhibit specific patterns. Eto et al. [START_REF] Eto | A proposal of malware distinction method based on scan patterns using spectrum analysis[END_REF] proposed a method to extract the features of scanning malwares based on the oscillation of destination IP addresses in the captured scan packets. Li et al. [START_REF] Li | Honeynet-based botnet scan traffic analysis[END_REF] proposed a general framework that identifies scanning events and analyzing methods used by botnets in probing compaigns. They applied their framework to extract the scanning characteristics of a set of 6 botnets. Papale et al. [START_REF] Dainotti | Analysis of a "/0" Stealth Scan from a Botnet[END_REF] analyzed a 12-day world-wide cyber scanning campaign targeting VoIP (SIP) servers caught by a /8 darknet. They found that the origin is the Sality botnet which generated about 20 million packet from roughly 3 million IP addresses.

Few studies explored the dependencies between targeted ports. McNutt and Markus [START_REF] Mcnutt | Correlations Between Quiescent Ports in Network Flows[END_REF] presented a method for detecting the start of anomalous port-specific activity by recognizing deviation from correlated activities. They found a high correlation between time series of flow counts on unassigned or obsolete ports that do not have active services. Therefore, they can detect ports receiving anomalous activities. In contrast to our work, they used in their study a traffic of an organization network (not a darknet traffic), where the amount of benign traffic is large, hiding thus malicious traffic. Lagraa and Franc ¸ois [START_REF] Lagraa | Knowledge discovery of port scans from darknet[END_REF] inferred the dependency between services using graph analysis. They proposed a graph-based model to discover port scanning behavior patterns. They applied methods utilized for community structure discovery in large graphs in order to identify clusters of ports. Our work generalizes their approach: instead of constructing graphs for each pair of source and destination IP addresses, our graphs aggregate the probing activity by source IP address. This approach shows the general exploration pattern followed by a network prober, regardless of the target host.

III. EXPLORATORY DATA ANALYSIS

A. Data Set

The data we use is collected by a /20 network telescope. The traffic was recorded from November 2014 until October 2017 and has a size of 2 TB. The collected traffic consists of timestamped packet headers. We record for each packet the source and the destination IP addresses, the source and the destination ports and the packet's flags. Our study focuses on stateful connections established from the source. Hence, we consider only packets with a TCP SYN flag which count for approximately 4.5 billion packets.

B. Traffic by Port

We begin by extracting the traffic received by each port. That is, we aggregate the received traffic by destination port and we count the number of TCP SYN packets for each destination port. Figure 1 shows the 30 most targeted ports. We observe that the most targeted services are remote access services, web servers, database management systems and some Microsoft services. The network probers tend to use alternative ports in addition to the official ones. The port 23 (telnet) generates more than 50% of the traffic. Figure 2

C. Top Network Probers' Interests

The intent of a network prober might be manifested in the services he targets. Knowing the ports that interest the top network probers determines services requiring particular security efforts. We consider as a top network prober one maintaining an average probing rate higher than 150 TCP SYN packets per day. It is noteworthy that our definition of top network probers does not include probers performing distributed probing activities.

First, we count TCP SYN packets sent by more than 64 million source IP addresses included in our data set. Then, for each top network prober (they are nearly 1500), we extract the probing rates received by each port and we aggregate the counts by port. Figure 3 shows the 30 most targeted ports by top network probers. In contrast to the results in Figure 1, top network probers focuse their probing activities on the port 22 (SSH) rather than the port 23 (telnet).

D. Traffic by Country

The distribution of the received traffic by geolocation helps determining how likely an occurred probing campaign is originating from a specific country. We extract the total traffic caught by the network telescope by country. We infer the country code from the source IP address using the DB-IP Fig. 4: Top 30 countries and their corresponding traffic percentages database 1 . Figure 4 summarizes the received traffic by country code.

IV. PROBING PATTERNS

In this section, we model the behavior of network probers using transition graphs. We assess the relationship between targeted services by determining the probability of transition from a port to another, then, we identify different network exploration patterns.

A. Graph Modeling

A network communication can be identified by a 5-tuple: source and destination IP addresses, source and destination ports, and the transport protocol. Our aim is to analyze for each network prober the exploration behavior of the whole darknet. Hence, we take into account only two features: the source IP address and the destination port.

To extract the graphs, we begin by aggregating the traffic generated by each network prober. Then, we count the number of transitions from a destination port to another by sequentially browsing the extracted traffic. These counts are finally normalized in order to get the transition probabilities. It is to note that the time dimension is omitted during this process.

Formally, we extract for each source IP address i a transition graph G i (V i , E i ), in which V i is the set of destination ports targeted by the network prober i and E i represents the transition probabilities between destination ports that are elements of E i . The association between two elements p a and p b of E i represents the probability that the network prober i switches from p a to p b .

B. Extracted Graphs

Figure 5 shows a sample of transition graphs corresponding to 3 network probers. Figure 5a represents a network prober sequentially targeting services typically deployed in web servers: SSH (22), RDP (3389), MySQL (3306) and FTP (21), while focusing on the HTTP (80) server. Figure 5b shows a network prober targeting only the MySQL server port in addition to two of its alternatives. Figure 5c corresponds to a network prober targeting remote access services such as SSH (22) and its alternative (2222), and telnet (23) and its alternative (2323). Many other probing patterns were identified but they are too large to fit in this paper.

The extracted transition graphs differ from each other by two main components: the number of vertices that corresponds to the destination ports and the number of edges describing the Fig. 7: Transition matrix of all network probers exploration behavior of a network prober. Figure 6 represents the cumulative distribution function of the number of target ports by individual network probers. We observe that more than 80% of network probers target less than five ports in the whole darknet space. This means that most attackers are focusing their probing activities only on services of interest, which might be related to a vulnerability disclosure for example.

C. Relationship Between Ports

In this section, we aim to identify the relationship between ports in terms of transition probabilities. We begin by aggregating by network prober the number of transitions from a port to another. Then, normalize these counts by the total number of transitions. We repeat this process for all network probers combined together and for top network probers (as defined in section III-C).

Figure 7 represents the transition matrix of the 30 most targeted ports in the whole darknet (see section III-B). Unsurprisingly, the figure shows a high association between ports and their alternatives: 23 and 2323, 80 and 8080, and 22 and 2222. The figure also emphasizes a strong relationship between services of the same type such as MS-SQL SERVER (1433) and MySQL (3306).

Similarly, Figure 8 shows the transition probabilities of the 30 most targeted ports by top network probers. We observe fewer relationships compared to the previous transition matrix.

V. PREDICTION OF PORT PROBING RATES

Predicting probing activities at the service level is a key feature for making better security operational decisions. Observing a significant disparity between the predicted probing rate and the actual value may help detecting an imminent threat. In this section, we forecast the probing rate of a target port by measuring its previous probing rates as well as the traffic received by the other ports. The predictions are performed one step ahead of time using the non-stationary autoregressive model (AR) and the non-stationary vector autoregressive model (VAR), for each port of the 550 most targeted ports (see Section III-B), for different time resolutions.

A. Data Set

The probing rate is inferred for different time resolutions: 1 hour, 3 hours, 6 hours, 12 hours and 1 day. For each time resolution, we extract a data set consisting of the probing rate time series of the most targeted ports D = {X (i) } i∈ [START_REF] Durumeric | An Internet-Wide View of Internet-Wide Scanning[END_REF]550] in which each record is a vector of probing rates occurring in the same time interval.

B. Non-stationarity of Probing Rate Time Series and Design Parameters

The analysis of the port probing rate time series showed non-stationarity of the first and the second order statistics over the period of 3 years. However, we observed that when considering shorter time windows, the non-stationarity tends to be alleviated at least in terms of average. Therefore, we introduce a "rolling window" over the probing rate time series. We train the estimators using the data falling in the rolling window, then the prediction is performed one step ahead of time. The size of the rolling window is a design parameter that could be interpreted as follows: a short rolling window allows to track the trend in time while larger rolling windows has the effect of averaging the trend. Another design parameter to consider is the autoregressive order which allows the model to infer the linear short term dependencies.

C. Non-stationary AR and VAR Models

To forecast probing rates, we use the non-stationary autoregressive model of order p given by:

x (i) (t; p) = w (i) 0 (t) + p h=1 w (i) h (t)x (i) (t -h) + (i) (t) (1)
where x (i) (t) is the probing rate received by the i th port at time t. (i) (t) is the white noise at time t. W

(i) t = (w (i) 0 (t), w (i) 1 (t), . . . , w (i) p (t))
is the vector of the model parameters. These parameters are estimated using data falling in the rolling window (see Section V-B) and they vary in time.

The non-stationary vector autoregressive model of order p is given by:

x (i) (t; p) = w (i) 0 (t)+ j∈I k p h=1 w (ij) h (t)x (j) (t-h)+ (i) (t) (2)
where x (i) (t) is the probing rate received by the i th port (the target port) at time t, I k is the set of indexes of the k retained features (see Section V-F), x (j) (t) is the probing rate on the j th selected port at time t, and (i) (t) is the white noise at time t. w [1,p] are the parameters of the model varying in time.

(i) 0 (t) and {w (ij) h (t)} j∈I k ,h∈

D. Training Algorithm

After transforming the time series into a supervised learning problem, the algorithm used to train the non-stationary AR and VAR regressors is the straightforward normal equation given by:

W (i) = (X T X) -1 X T X (i) (3)
where W (i) is the vector of the trainable weights, X = (X (j) ) j∈I k is the probing rates feature matrix and X (i) is the probing rates response vector. In order to reduce the computation complexity, no regularization is used.

E. Design Parameter Selection

The design parameters, namely the size of the rolling window N and the autoregressive order p, are determined using a grid search strategy. We varied p in [START_REF] Durumeric | An Internet-Wide View of Internet-Wide Scanning[END_REF]10] for the 5 considered time resolutions. Then, we tried an exhaustive set of rolling window sizes for each autoregressive order p. The range of N starts with 10 × p time units and ends with 75% of the time series length (leaving thus at least 25% of data for validation) with an increment of 10 time units. The optimal design parameter values p and N are given by the estimator providing the best coefficient of determination R 2 .

F. Feature Selection for the Non-stationary VAR Model

To improve the performance of the non-stationary VAR estimators, we select features according to their individual effect on the response variable using the Pearson correlation coefficient. This process has the effect of reducing the noise introduced by uncorrelated features. First, we split the data falling in the rolling window (the one giving the best nonstationary AR performance) into two subsets: a feature selection set F including 75% of the data and a validation set V.

Second, we compute on F the univariate correlations in term of probing rates between the target time series and the time series serving as features to the non-stationary VAR model, including the autoregressive features. Third, we iteratively select the k most correlated features which we use to train nonstationary VAR model on F and we calculate the coefficient of determination R 2 on the validation set V. The optimal set of features given by our feature selection strategy is the one providing the best coefficient of determination. It is worth mentioning that the selected features may change over time based on the location of the rolling window in the time series. Finally, the selected features are scaled to zero mean and to unit standard deviation.

G. Results and Discussion

Table I summarizes the performance of the non-stationary AR and VAR models for 5 different time resolutions for a set of popular services. It also includes the optimal design parameters for the non-stationary AR estimators. We used the same design parameters for the non-stationary VAR estimators.

The performance of the regressors tends to increase for larger time resolutions, for all the ports except the telnet service. This is due to lowered stochasticity of probing rates when considering larger time resolutions. Also, we observe that the probing rates of remote access services are the most predictable. The reason is that such services are highly targeted by network probers and their probing rate time series are stationary when considering short time resolutions. Also, we observe that the non-stationary AR model produces satisfying results for services exhibiting low short term probing rates variability such as telnet (ports 23 and 2323). Figure 9 shows that non-stationary VAR model consistently produces better results for services exhibiting high probing rate variability such as the web services (ports 80 and 443) and the database management systems (ports 1433 and 3306).

It is noteworthy that the non-stationary autoregressive model fails in predicting abrupt probing rate changes because of its persistence property. More powerful and stable models such as FARIMA+GARCH could be used to predict these extreme values if the probing rate time series exhibit longrange dependence phenomenon [START_REF] Zhan | Characterizing honeypot-captured cyber attacks: Statistical framework and case study[END_REF], [START_REF]Predicting cyber attack rates with extreme values[END_REF]. Also, such nonstationary AR and VAR models, as defined in our paper, require constant update of their parameters (the trainable weights) and their hyperparameters (the selected features) due to the non-stationarity of the probing rate time series, which is computationally expensive.

VI. CONCLUSION

This work presented an exploratory data analysis performed on 2 TB of traffic collected by a network telescope during the period of 3 years. The investigation of the network telescope traffic showed that 90% of probing activities are targeting only 550 ports of the port space. The latter include remote access services which are the most sought by network probers, followed by database management systems, web services and The second task was about inferring network probers reconnaissance patterns. We modeled these exploration behaviors using transition graphs showing the probabilities of switching from a port to another. This would be exploited in different applications such as port clustering or classifying network probers based on their exploration behaviors.

Finally, we assessed to which extent the non-stationary autoregressive and vector autoregressive models could produce reliable short term probing rate predictions at the port level. Due to their short memory property, such models could be used to learn non-stationary probing rate processes with short term persistence. However, when probing rate processes exhibit long-range dependencies, more robust models could be utilized such as GRU and LSTM recurrent neural networks. ACKNOWLEDGMENT This research work is part of the ThreatPredict2 project partly funded by the NATO Science for Peace and Security (SPS) programme under research contract SPS G5319 "ThreatPredict: From Global Social and Technical Big Data to Cyber Threat Forecast". We acknowledge the support from the National Center of Scientific and Technical Research (CNRST), Rabat, for the grant of an excellence scholarship. The authors would like to thank Frédéric Beck from the High

Fig. 1 :Fig. 2 :

 12 Fig. 1: The 30 most targeted ports and their corresponding traffic percentage (the length of the bars are in log scale and the labels above the bars are the actual values)

Fig. 5 :

 5 Fig. 5: A sample of transition graphs of 3 network probers. The size of the vertices corresponds to the number of TCP SYN packets received by the targeted port.

  Cumulative Distribution of the Number of Targeted Ports by Network Probers

Fig. 6 :

 6 Fig. 6: Cumulative distribution of the number of ports targeted by network probers

Fig. 8 :

 8 Fig. 8: Transition matrix of top network probers

Fig. 9 :

 9 Fig. 9: Comparison of performances of the non-stationary AR and VAR models

TABLE I :

 I Performances of non-stationary AR and VAR estimators for different time resolutions for a set of popular services. p and w are the optimal design parameters for the non-stationary AR estimators.Performances of theNon-stationary AR and VAR Models for Different Time Resolutions AR, 1 hour time resolution AR, 3 hours time resolution AR, 6 hours time resolution AR, 12 hours time resolution AR, 1 day time resolution VAR gain in performance

	R 2	0.0 0.1 0.2 0.3 0.4 0.5 0.6 1.0 0.9 0.8 0.7	23	2323	445	22	1883	Port	2222	1433	443	80	3306

https://threatpredict.loria.fr Security Laboratory at INRIA Nancy-Grand Est, LORIA, for his efforts in managing data and computation servers.


Agnès Helmstetter

Didier Sornette

Foreshocks Explained by Cascades of Triggered Seismicity

The observation of foreshocks preceding large earthquakes and the suggestion that foreshocks have specific properties that may be used to distinguish them from other earthquakes have raised the hope that large earthquakes may be predictable. Among proposed anomalous properties are the larger proportion than normal of large versus small foreshocks, the power law acceleration of seismicity rate as a function of time to the mainshock and the spatial migration of foreshocks toward the mainshock, when averaging over many sequences. Using Southern California seismicity, we show that these properties and others arise naturally from the simple model that any earthquake may trigger other earthquakes, without arbitrary distinction between foreshocks, aftershocks and mainshocks. We find that foreshocks precursory properties are independent of the mainshock size. This implies that earthquakes (large or small) are predictable to the same degree as seismicity rate is predictable from past seismicity by taking into account cascades of triggering. The cascades of triggering give rise naturally to long-range and long-time interactions, which can explain the observations of correlations in seismicity over surprisingly large length scales.

Foreshocks, mainshocks and aftershocks: hypothesis and predictions

It has been recognized for a long time that large earthquakes are sometimes preceded by an acceleration of the seismic activity, known as foreshocks [START_REF] Jones | Frequency of foreshocks[END_REF][START_REF] Abercrombie | Occurrence patterns of foreshocks to large earthquakes in the western United States[END_REF]. In addition to the increase of the seismicity rate a few hours to months before large earthquakes, other properties of foreshocks have been reported, which suggest their usefulness (when present) as precursory patterns for earthquake prediction. Special physical mechanisms have been proposed for foreshocks with the hope of helping earthquake prediction hope of helping earthquake prediction [START_REF] Yamashita | A model of foreshock occurrence[END_REF][START_REF] Sornette | Statistical model of earthquake foreshocks[END_REF][START_REF] Dodge | Detailed observations of California foreshock sequence: implications for the earthquake initiation process[END_REF][START_REF] Dieterich | Implications of fault constitutive properties for earthquake prediction[END_REF]. In addition, anomalous precursory seismic activity extending years to decades before large earthquakes and at distances up to ten times the mainshock rupture size are often thought to require different physical mechanisms [START_REF] Keilis-Borok | One regularity in the occurrence of strong earthquakes[END_REF][START_REF] Knopoff | Increased long-range intermediate-magnitude earthquake activity prior to strong earthquakes in California[END_REF][START_REF] Bowman | An Observational test of the critical earthquake concept[END_REF][START_REF] Jaumé | Evolving towards a critical point: A review of accelerating seismic moment/energy release prior to large and great earthquakes[END_REF][START_REF] Sammis | Positive feedback, memory and the predictability of earthquakes[END_REF][START_REF] Keilis-Borok | Earthquake prediction: State-of-the-art and emerging possibilities[END_REF] than for foreshocks closer to the mainshock epicenters.

The division between foreshocks, mainshocks, and aftershocks has a long and distinguished history in seismology. Within a pre-specified space-time domain, foreshocks are usually defined as earthquakes (above the background rate) preceding a larger earthquake (mainshock), which is itself followed by an increase in seismicity of smaller earthquakes (aftershocks). However, recent empirical and theoretical scrutiny suggests that this division might be arbitrary and physically artificial [START_REF] Kagan | Stochastic synthesis of earthquake catalogs[END_REF][START_REF] Shaw | Generalized Omori law for aftershocks and foreshocks from a simple dynamics[END_REF][START_REF] Jones | Are foreshocks mainshocks whose aftershocks happen to be big? Evidence From California and Italy[END_REF][START_REF] Hough | Aftershocks; Are they earthquakes or afterthoughts? EOS[END_REF]Felzer et al., 2002]. Since the underlying physical processes are not fully understood, the qualifying time and space windows used to select aftershocks, mainshocks and aftershocks are more based on common sense than on hard science. If the spacetime window is extended and a new event not considered previously is found with a magnitude larger than the previously classified mainshock, it becomes the new mainshock and all preceding events are retrospectively called foreshocks. A clear identification of foreshocks, aftershocks and mainshocks is hindered by the fact that nothing distinguishes them in their seismic signatures: at the present level of resolution of seismic inversions, they are found to have the same double-couple structure and the same radiation patterns [START_REF] Hough | Aftershocks; Are they earthquakes or afterthoughts? EOS[END_REF]. Statistically, the aftershock magnitudes are distributed according to the Gutenberg-Richter (GR) distribution P (m) ∼ 10 -bm with a b-value similar to other earthquakes [START_REF] Ranalli | A statistical study of aftershock sequences[END_REF][START_REF] Utsu | The centenary of the Omori formula for a decay law of aftershock activity[END_REF]. However, some studies [START_REF] Knopoff | b-values for foreand aftershocks in real and simulated earthquakes sequences[END_REF][START_REF] Molchan | Statistical analysis of the results of earthquake prediction, based on bursts of aftershocks[END_REF][START_REF] Molchan | Immediate foreshocks: time variation of the b-value[END_REF] have suggested that foreshocks have a smaller b-value than other earthquakes, but the physical mechanisms are not yet understood.

The Omori law describes the power law decay of the aftershock rate ∼ 1/(t -t c ) p with time from a mainshock that occurred at t c [START_REF] Omori | On the aftershocks of earthquakes[END_REF][START_REF] Kagan | Statistical study of the occurrence of shallow earthquakes[END_REF][START_REF] Utsu | The centenary of the Omori formula for a decay law of aftershock activity[END_REF], and which may last from months up to decades. In contrast with the well-defined Omori law for aftershocks, there are huge fluctuations of the foreshock seismicity rate, if any, from one sequence of earthquakes to another one. Figure 1 shows foreshock sequences of all M ≥ 6.5 mainshocks in the catalog of Southern California seismicity provided by the Southern California Seismic Network for the period 1932-2000. By stacking many foreshock sequences, a well-defined acceleration of the seismicity preceding mainshocks emerges, quantified by the so-called inverse Omori law ∼ 1/(t c -t) p ′ , where t c is the time of the mainshock [START_REF] Papazachos | The time distribution of reservoirassociated foreshocks and its importance to the prediction of the principal shock[END_REF][START_REF] Kagan | Statistical study of the occurrence of shallow earthquakes[END_REF][START_REF] Jones | Some characteristics of foreshocks and their possible relationship to earthquake prediction and premonitory slip on fault[END_REF]. While we see clearly an acceleration for the averaged foreshock number (black line in Figure 1), there are huge fluctuations of the rate of foreshocks for individual sequences. Most foreshock sequences are characterized by the occurrence of a major earthquake before the mainshock, which has triggered the mainshock most probably indirectly due to a cascade of multiple triggering. For instance, 66 days before its occurrence, the M = 7.3 Landers earthquake (upper curve in Figure 1) was preceded by the M = 6.1 Joshua-Tree earthquake. The successive oscillations of the cumulative number of events after the Joshua-Tree earthquake correspond to secondary, tertiary, etc., bursts of triggered seismicity. The inverse Omori law is usually observed for time scales shorter than the direct Omori law, of the order up to weeks to a few months before the mainshock. However, there seems to be no way of identifying foreshocks from usual aftershocks and mainshocks in real time (see [START_REF] Jones | Are foreshocks mainshocks whose aftershocks happen to be big? Evidence From California and Italy[END_REF][START_REF] Hough | Aftershocks; Are they earthquakes or afterthoughts? EOS[END_REF] for a pioneering presentation of this view point). In other words, is the division between foreshocks, mainshocks and aftershocks falsifiable [START_REF] Kagan | Is earthquake seismology a hard, quantitative science?[END_REF]?

In order to address this question, we present a novel analysis of seismic catalogs, based on a parsimonious model of the foreshocks, mainshocks and aftershocks, in terms of earthquake triggering: earthquakes may trigger other earthquakes through a variety of physical mechanisms [START_REF] Harris | Earthquake stress triggers, stress shadows, and seismic hazard[END_REF] but this does not allow one to put a tag on them. Thus, rather than keeping the specific classification that foreshocks are precursors of mainshocks and mainshocks trigger aftershocks, we start from the following hypothesis:

Hypothesis Foreshocks, mainshocks, and aftershocks are physically indistinguishable.

In this paper, we study a simple model of seismicity based on this Hypothesis and demonstrate that it can explain many properties of foreshocks, such as the larger proportion than normal of large versus small foreshocks, the power law acceleration of seismicity rate as a function of time to the mainshock and the spatial migration of foreshocks toward the mainshock, when averaging over many sequences. These properties and others arise naturally from the simple model that any earthquake may trigger other earthquakes, without arbitrary distinction between foreshocks, aftershocks and mainshocks.

Definition of the model

The simplest construction that embodies the Hypothesis is the epidemic-type aftershock sequence (ETAS) model introduced in [START_REF] Kagan | Stochastic synthesis of earthquake catalogs[END_REF][START_REF] Kagan | Statistical short-term earthquake prediction[END_REF] (in a slightly different form) and in [START_REF] Ogata | Statistical models for earthquake occurrence and residual analysis for point processes[END_REF]. In this model, all earthquakes may be simultaneously mainshocks, aftershocks and possibly foreshocks. An observed "aftershock" sequence in the ETAS model is the sum of a cascade of events in which each event can trigger more events. The triggering process may be caused by various mechanisms that either compete or combine, such as 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, etc.. The epidemic-type aftershock (ETAS) model assumes that a given event (the "mother") of magnitude m i occurring at time t i and position r i gives birth to other events ("daughters") of any possible magnitude m at a later time between t and t + dt and at point r ± dr at the rate

φ mi (t -t i , r -r i ) = ρ(m i ) Ψ(t -t i ) Φ( r -r i ) . (1)
We will refer to φ mi (t -t i , r -r i ) as the "local"' Omori law, giving the seismic rate induced by a single mother. It is the product of three independent contributions: 1. ρ(m i ) gives the number of daughters born from a mother with magnitude m i . This term is in general chosen to account for the fact that large earthquakes have many more triggered events than small earthquakes. Specifically,

ρ(m i ) = K 10 αmi , (2)
which is in good agreement with the observations of aftershocks productivity (see Helmstetter [2003] and references therein).

2. Ψ(t -t i ) is a normalized waiting time distribution giving the rate of daughters born at time t -t i after the mother

Ψ(t) = θ c θ (t + c) 1+θ .
(3)

The "local" Omori law Ψ(t) is generally significantly different from the observed Omori law of aftershocks. We have shown in [Helmstetter and Sornette, 2002a] that the exponent p of the observed "global" Omori law is equal to or smaller than the exponent 1+θ of the "local" Omori law Ψ(t) due to the cascades of triggering.

3. Φ( r -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 | r -r i | from the mother. Specifically, we take

Φ( r) = µ d( | r| d + 1) 1+µ , ( 4
)
where µ is a positive exponent and d is a characteristic dimension of the earthquake source, as suggested by [START_REF] Ogata | Seismicity analysis through point-process modeling : a review[END_REF].

The last ingredient of the ETAS model is that the magnitude m of each daughter is chosen independently from that of the mother and of all other daughters according to the Gutenberg-Richter distribution

P (m) = b ln(10) 10 -b(m-m0) , (5)
with a b-value usually close to 1. m 0 is a lower bound magnitude below which no daughter is triggered.

The total seismicity rate is the sum of all aftershocks of past events and of a constant external source s used to model the constant tectonic loading. The full distribution of seismicity in time, space and magnitude is thus described by the 8 parameters of the ETAS model K, c, θ, µ, d, α, b and s. See [START_REF] Ogata | Seismicity analysis through point-process modeling : a review[END_REF], [Helmstetter and Sornette, 2002a] and references therein for a discussion of the values of these parameters in real seismicity.

The ETAS model has been used previously to give short-term probabilistic forecast of seismic activity [START_REF] Kagan | Statistical short-term earthquake prediction[END_REF][START_REF] Kagan | Probabilistic forecasting of earthquakes[END_REF][START_REF] Console | A simple and testable model for earthquake clustering[END_REF], and to describe the temporal and spatial clustering of seismic activity [START_REF] Ogata | Statistical models for earthquake occurrence and residual analysis for point processes[END_REF][START_REF] Kagan | Likelihood analysis of earthquake catalogues[END_REF][START_REF] Kagan | Probabilistic forecasting of earthquakes[END_REF][START_REF] Console | A simple and testable model for earthquake clustering[END_REF]Felzer et al., 2002]. The Hypothesis and the ETAS model allow us to study two classes of foreshocks, called of type I and of type II defined in the next section.

1.2. Definition of foreshocks of type I and of type II 1.2.1. Formal definitions. The usual definition of foreshocks, that we shall call "foreshock of type I," refers to any event of magnitude smaller than or equal to the magnitude of the following event, then identified as a "mainshock." This definition implies the choice of a space-time window R × T used to define both foreshocks and mainshocks. Mainshocks are large earthquakes that were not preceded by a larger event in this space-time window.

In contrast, the Hypothesis that the same physical mechanisms describe the triggering of a large earthquake by a smaller earthquake (mainshock triggered by a type I foreshock) and the triggering of a small earthquake by a larger event (aftershock triggered by a mainshock) makes it natural to remove the constraint that foreshocks must be smaller than the mainshock. We thus define "foreshock of type II" as any earthquake preceding a large earthquake which is defined as the mainshock, independently of the relative magnitude of the foreshock compared to that of the mainshock. That is, foreshocks of type II may actually be larger than their mainshock.

The advantage of this second definition is that foreshocks of type II are automatically defined as soon as one has identified the mainshocks, for instance, by calling mainshocks all events of magnitudes larger than some threshold of interest. Foreshocks of type II are thus all events preceding these large magnitude mainshocks in a space-time window R × T . In con-trast, foreshocks of type I need to obey a constraint on their magnitude, which may be artificial, as we shall see further down.

1.2.2. Practical implementation. In our analysis of the SCEC catalog, we construct foreshock and aftershock sequences as follows. A mainshock is defined as an earthquake in the magnitude range (M, M +∆M ) that was not preceded by a larger event in a space-time window (R 2 , T ) before the mainshock. The distance R 2 ∼ = 50 km is here chosen to be close to (but smaller than) the maximum size of the spatial clusters of seismicity in the California catalog, in order to minimize the influence of large earthquakes that occurred before the mainshock. Other choices between 20 km to 200 km have been tested and give essentially the same results. The aftershocks are all events that occurred in a space-time window (R, T ) after each mainshock. The foreshocks of type I are selected in a space-time window R, T before each mainshock. Foreshocks of type II are selected in the same space-time window (R, T ) before each mainshock, now defined without the constraint that they were not preceded by a larger event in a space-time window (R 2 , T ). The only difference between foreshocks of type I and of type II is that the selection of their respective mainshocks is different.

The minimum magnitude cutoff used for aftershocks and foreshocks is M 0 = 3. The fact that the SCEC catalog is not complete below M = 3.5 for the first part of the catalog 1932-1975 does not affect the results on the inverse and direct Omori laws, it simply under-estimates the number of foreshocks and aftershocks but does not change the temporal evolution of the rate of seismicity before and after mainshocks. The incompleteness of the catalog for M < 3.5 before 1975 explains the roll-off of the magnitude distribution for small magnitudes. Even for the most recent part, the SCEC catalog is not complete above M = 3 in the few hours or days after a large M ≥ 6 earthquake, due to the saturation of the seismic network. These missing events may induce a spurious roll-off of both the direct and inverse Omori laws at short times before or after a mainshock, but cannot induce a spurious acceleration of the rate of seismicity before mainshocks.

We stack all foreshocks and aftershocks sequences synchronized at the time of the "mainshocks" in different mainshock magnitude intervals [M ,M + ∆M ], and for different choices of the space-time window R, T used to define foreshocks and aftershocks. We use larger magnitude intervals ∆M for larger mainshock magnitudes M to compensate for the smaller number of large mainshocks. R has been tested between 10 km and up to 500 km with no essential change, except for an increasing sensitivity to the background seismicity for the largest R (see Figure 3). T has been tested between 0.5 year to 10 years with similar results. Tests have also been performed with the spatial window size R adjusted to scale with the mainshock magnitude with no significant difference. Our results presented below thus appear robust with respect to the (arbitrary) definitions of the space-time windows and the definition of mainshocks. In the main text, foreshocks of type I and of type II are treated separately. Previous studies of foreshocks using a stacking method [Papazachos, 1975a;Papazachos, 1975b;[START_REF] Jones | Frequency of foreshocks[END_REF][START_REF] Kagan | Statistical study of the occurrence of shallow earthquakes[END_REF][START_REF] Jones | Some characteristics of foreshocks and their possible relationship to earthquake prediction and premonitory slip on fault[END_REF][START_REF] Shaw | Generalized Omori law for aftershocks and foreshocks from a simple dynamics[END_REF][START_REF] Reasenberg | Foreshocks occurrence before large earthquakes[END_REF] have considered foreshocks of type I only.

Foreshock properties derived from the ETAS model

The simple embodiment of the Hypothesis in the ETAS model leads to the following consequences and predictions [START_REF] Helmstetter | Mainshocks are aftershocks of conditional foreshocks: how do foreshock statistical properties emerge from aftershock laws[END_REF], which are proposed as crucial tests of the Hypothesis.

1. The rate of foreshocks of type II is predicted to increase before the mainshock according to the inverse Omori law N (t) ∼ 1/(t c -t) p ′ with an exponent p ′ smaller than the exponent p of the direct Omori law. The exponent p ′ depends on the "local" Omori exponent 1 + θ describing the direct triggering rate between earthquakes (first-generation triggering), on the b-value of the GR distribution and on the exponent α quantifying the increase ∝ 10 αM in the number of aftershocks as a function of the magnitude M of the mainshock [Helmstetter, 2003]. The inverse Omori law also holds for foreshocks of type I preceding large mainshocks. The inverse Omori law results from the direct Omori law for aftershock and from the multiple cascades of triggering. The inverse Omori law emerges as the expected (in a statistical sense) trajectory of seismicity, conditioned on the fact that it leads to the burst of seismic activity accompanying the mainshock.

2. In contrast with the direct Omori law, which is clearly observed after all large earthquakes, the inverse Omori law is observed only when stacking many foreshock sequences. Even for small mainshocks, for which the number of aftershocks is similar to the number of foreshocks, there are much larger fluctuations in the rate of foreshocks before mainshocks than in the rate of aftershocks.

3. While the number of aftershocks increases as 10 αM with the magnitude M of the mainshock, the number of foreshocks of type II is predicted to be independent of M . Thus, the seismicity should increase on average according to the inverse Omori law before any earthquake, whatever its magnitude. For foreshocks of type I, the same results should hold for large mainshocks. For small and intermediate values of the mainshock magnitude M , the conditioning on foreshocks of type I to be smaller than their mainshock makes their number increase with M solely due to the constraining effect of their definition.

4. The GR distribution for foreshocks is predicted to change upon the approach of the mainshock, by developing a bump in its tail. Specifically, the modification of the GR law is predicted to take the shape of an additive correction to the standard power law, in which the new term is another power law with exponent b -α. The amplitude of this additive power law term is predicted to also exhibit a power law acceleration upon the approach to the mainshock.

5. The spatial distribution of foreshocks is predicted to migrate toward the mainshock as the time increases toward the time of the mainshock, by the mechanism of a cascade of seismic triggering leading to a succession of jumps like in a continuous-time random walk [Helmstetter and Sornette, 2002b].

We now proceed to test systematically these predictions on the catalog of the Southern California Data Center (SCEC) over the period 1932-2000, which is almost complete above M = 3 and contains more than 22000 M ≥ 3 earthquakes, using the methodology described in section 1.2.

Properties of foreshocks, mainshocks and aftershocks and comparison with the ETAS model

Observation of the direct and inverse Omori law

Figure 2 shows the rate of foreshocks of type II as a function of t c -t and of aftershocks as a function of t -t c , where t c is the time of the mainshocks, for different mainshock magnitude between 3 and 7. Aftershocks and foreshocks have been selected with the space-time window T = 1 yr and R = 50 km. Both rates follow an approximate power law (inverse Omori law for foreshocks with exponent p ′ and Omori law for aftershocks with exponent p).

The fluctuations of the rate of foreshocks are larger for large mainshocks, because the number of mainshocks (resp. foreshocks) decreases from 15584 (resp. 1656249) for the magnitude range 3 -3.5 down to 47 (resp. 1899) for M > 6 mainshocks. In contrast, the fluctuation of the rate of aftershocks are larger for small mainshocks magnitudes, due to the increase of the number of aftershocks per mainshock with the mainshock magnitude, and to the rules of mainshock selection which reject a large proportion of small earthquakes. The number of mainshocks is different for aftershocks and foreshocks due to their distinct definition. The number of mainshocks associated with aftershocks decreases from 677 for the magnitude range 3-3.5 down to 39 for M ≥ 6 mainshocks. The number of aftershocks increases from 2614 for 3 ≤ M < 3.5 up to 7797 for M ≥ 6 mainshocks. The fluctuations of the data makes it hard to exclude the hypothesis that the two power laws have the same exponent p = p ′ ∼ = 1., even if p ′ seems slightly smaller. Note that it can be shown theoretically for α < b/2 that p = 1 -θ and p ′ = 1 -2θ for a local Omori law with exponent 1 + θ and that the difference p -p ′ should get smaller as α increases above b/2 [START_REF] Helmstetter | Mainshocks are aftershocks of conditional foreshocks: how do foreshock statistical properties emerge from aftershock laws[END_REF]. Since α ∼ = 0.8 > b/2 ≈ 0.5, this limit is met which explains the smallness of the difference p -p ′ . The truncation of the seismicity rate for small times |t c -t| < 1 day, especially for aftershocks of large M > 6 mainshocks and for foreshocks, is due to the incompleteness of the catalog at very short times after mainshocks due to the saturation of the seismic network. At large times from the mainshock, the seismicity rate decreases to the level of the background seismicity, as seen clearly for the rate of aftershocks following small M = 3 mainshocks.

The second striking observation is the strong vari-ation of the amplitude of the rate N a (t) of aftershocks as a function of the magnitude M of the mainshock, which is well-captured by an exponential dependence N a (t) ∝ 10 αM /(t -t c ) p with α ∼ = 0.8 [Helmstetter, 2003]. In contrast, the rates of foreshocks of type II are completely independent of the magnitude M of the mainshocks: quite strikingly, all mainshocks independently of their magnitudes are preceded by the same statistical inverse Omori law, with the same power law increase and the same absolute amplitude! All these results are very well modeled by the ETAS model with the parameters α = 0.8, θ = 0.2 and b = 1 using the theoretical framework and numerical simulations developed in [START_REF] Helmstetter | Mainshocks are aftershocks of conditional foreshocks: how do foreshock statistical properties emerge from aftershock laws[END_REF].

Another remarkable observation is presented in Figure 3 which shows the rate of foreshocks of type II for mainshock magnitudes between 4 and 4.5, for different values of the distance R used to select aftershocks and foreshocks. The inverse Omori law is observed up to R ≈ 200 km, and the duration of the foreshock sequences increases as R decreases due to the decrease of the effect of the background seismicity. Restricting to the shortest distances R to minimize the impact of background seismicity, the inverse Omori laws can be observed up to 10 yrs before mainshocks, for foreshocks of type II (Figure 4). Thus, foreshocks are not immediate precursors of mainshocks but result from physical mechanisms of earthquake triggering acting over very long times and large distances.

An important question concerns the relative weight of coincidental shocks, i.e., early aftershocks triggered by a previous large earthquake, which appear as foreshocks of type II to subsequent aftershocks (seen as mainshocks of these foreshocks of type II). Such coincidental shocks can give rise to an apparent inverse Omori law [START_REF] Shaw | Generalized Omori law for aftershocks and foreshocks from a simple dynamics[END_REF] when averaging over all possible positions of "mainshocks" in the sequence, without any direct interaction between these mainshocks and preceding events viewed as their foreshocks. Actually, these coincidental shocks form a minority of the total set, because the fraction of shocks directly triggered by a mainshock decays to negligible values beyond a few days for the range of parameters of the ETAS that realistically fit the SCEC catalog [Felzer et al., 2002].

Figure 5 shows the rate of foreshocks of type I as a function of t c -t and of aftershocks as a function of t -t c , where t c is the time of the mainshocks. The total number of foreshocks of type I is much smaller that the number of type II foreshocks for small main-shocks because a significant fraction of foreshocks of type II are "aftershocks" of large M > 6 earthquakes according to the usual definition and are therefore rejected from the analysis of foreshocks of type I, which are constrained to be smaller than their mainshock. There are much larger fluctuations for foreshocks of type I than for foreshocks of type II due to the smaller number of the former. Nevertheless, type I and type II foreshocks clearly follow the same trajectory of increased activity before a mainshock, suggesting that type I foreshocks, like type II foreshocks, are triggers of the mainshock.

The exponent p ′ of the inverse Omori law for foreshocks of type I is approximately equal to the exponent of foreshocks of type II and to the exponent p of the direct Omori law for aftershocks. The rate of foreshocks of type I increases slowly with the mainshock magnitude but this increase is not due to a larger predictability of larger earthquakes, as expected for example in the critical point theory [START_REF] Sammis | Positive feedback, memory and the predictability of earthquakes[END_REF] and as observed in a numerical model of seismicity [START_REF] Huang | Precursors, aftershocks, criticality and self-organized criticality[END_REF]]. The increase of the number of type I foreshocks with the mainshock magnitude can be reproduced faithfully in synthetic catalogs generated with the ETAS model and is nothing but the consequence of the algorithmic rules used to define foreshocks of type I. Namely, by definition, bigger mainshocks are allowed to have bigger type I foreshocks and therefore they have more type I foreshocks than smaller mainshocks. In other words, there is no physics but only statistics in the weak increase of foreshocks of type I with the mainshock magnitude. Confirming this concept, the inverse Omori law for foreshocks of type I becomes independent of the mainshock magnitudes M for large M , for which the selection constraint has only a weak effect. The dependence of the inverse Omori law for foreshocks of type I as a function of distance R used to select foreshocks is very similar (not shown) to that shown for foreshocks of type II in Figure 3. However, the duration of foreshock sequences is shorter for type I foreshocks because the number of type I foreshocks is smaller than the number of type II foreshocks due to the rules of foreshock selection.

Modification of the magnitude distribution before a mainshock

We have shown in [START_REF] Helmstetter | Mainshocks are aftershocks of conditional foreshocks: how do foreshock statistical properties emerge from aftershock laws[END_REF]] that selecting and stacking foreshock sequences in the ETAS model leads to a modification of the magnitude distribution compared to the theoretical distri-bution. We refer to [START_REF] Helmstetter | Mainshocks are aftershocks of conditional foreshocks: how do foreshock statistical properties emerge from aftershock laws[END_REF] for the derivation of the results stated below, which is too involved to be reported here. The distribution P (m) of foreshock magnitudes is predicted to get an additive power law contribution q(t)dP (m) with an exponent b ′ smaller than b and with an amplitude q(t) growing as a power law of the time to the mainshock:

P (m) = (1 -q(t))P 0 (m) + q(t) dP (m) , (6)
where P 0 (m) is the standard GR distribution P 0 (m) ∼ 10 -bm and dP (m) ∼ 10 -b ′ m with b ′ = b -α. The amplitude q(t) of the additive distribution dP (m) in ( 6) should increase as a power-law of the time to the mainshock according to

q(t) ∼ 1/(t c -t) θ b ′ α . ( 7
)
This analytical prediction has been checked with extensive numerical simulations of the ETAS model. This change of the magnitude distribution for foreshocks does not mean that foreshocks belong to a different population, but simply results from the definition of foreshocks, which are only defined as foreshocks after the mainshock occurred. The mechanism giving rise to the change of distribution for foreshocks is explained in [START_REF] Helmstetter | Mainshocks are aftershocks of conditional foreshocks: how do foreshock statistical properties emerge from aftershock laws[END_REF]. Intuitively, the modification of the magnitude distribution of foreshocks results from the increase of the number of triggered events with the mainshock magnitude. There are few large earthquakes, but they trigger many more earthquakes than smaller earthquakes. As a consequence, a large fraction of mainshocks are triggered, directly or indirectly, by large foreshocks. The proportion of large foreshocks is thus larger than the proportion of large earthquakes in the whole population, and gives an apparent lower b-value for foreshocks than for other earthquakes. The magnitude distribution of triggering events is given by the product ρ(m)P 0 (m) ∼ 10 αm 10 -bm ∼ dP (m), which gives the distribution of foreshock magnitudes for large m. We now test this prediction using the SCEC catalog on foreshocks of type II of M > 3 mainshocks, selected using R = 20 km and T = 1 yr. The magnitude distribution P (m) sampled at different times before mainshocks is shown in panel (a) of Figure 6, where the black to gray curves correspond to times preceding mainshocks decreasing from 1 year to 0.01 day with a logarithmic binning. As time approaches that of the mainshocks, one can clearly observe that the tails depart more and more from the standard GR power law P 0 (m) with b = 1.0 ± 0.1 estimated using the whole catalog and shown as the dashed line. The perturbation q(t)dP (m) in ( 6) of the foreshock magnitude distribution, shown in panel b) of Figure 6, can be estimated by fitting the prediction (6) to the observed magnitude distribution of foreshocks, by inverting the parameters b ′ and q(t) in ( 6) for different times before mainshocks. The obtained additive GR laws q(t)dP (m) are compatible with pure power laws with an approximately constant exponent b ′ = 0.6 ± 0.1 shown in panel (d) of Figure 6, except at very long times before the mainshock where it drops to 0 when the amplitude q(t) of the perturbation becomes too small. The amplitude q(t) is shown in panel (c) and is compatible with a power law (7) with a fitted exponent 0.3 ± 0.2. These observations are in good qualitative agreement with the predictions ( 6) and ( 7) on the nature of the modification of the GR law for foreshocks in terms of a pure additive power law perturbation with an amplitude growing as a power law of the time to the mainshocks. Quantitatively, b ′ is marginally outside the 2σ-confidence interval for the prediction b ′ = b -α = 0.2 ± 0.2 using the estimation α = 0.8 ± 0.1 given by [Helmstetter, 2003]. We attribute this discrepancy to the dual impact of the incompleteness of the catalog for small magnitudes after a large earthquake and to the smallness of the statistics. We stress that the prediction (6) with b ′ = b -α has been verified with good precision in synthetic catalogs which do not have these limitations [START_REF] Helmstetter | Mainshocks are aftershocks of conditional foreshocks: how do foreshock statistical properties emerge from aftershock laws[END_REF]. Using the best fitted value b ′ = 0.6, we obtain a reasonable agreement for the predicted exponent θb ′ /α and the fitted value 0.3±0.2 for the power law behavior of q(t) using θ in the range 0.2 -0.4.

Note that expression (6) contains as a special case the model in which the modification of the GR law occurs solely by a progressive decrease of the b-value as the time of the mainshock is approached (by putting q(t) = 1 and allowing b ′ to adjust itself as a function of time), as proposed in [START_REF] Berg | Relation between earthquake foreshocks, stress and mainshocks[END_REF][START_REF] Kagan | Statistical study of the occurrence of shallow earthquakes[END_REF][START_REF] Molchan | Statistical analysis of the results of earthquake prediction, based on bursts of aftershocks[END_REF][START_REF] Molchan | Immediate foreshocks: time variation of the b-value[END_REF]. Our quantitative analysis clearly excludes this possibility while being completely consistent with the mechanism embodied by the concept of triggered seismicity. Although the foreshock magnitude distribution is not a pure power-law but rather the sum of two power laws, our results rationalize the reported decrease of b-value before mainshocks [START_REF] Berg | Relation between earthquake foreshocks, stress and mainshocks[END_REF][START_REF] Kagan | Statistical study of the occurrence of shallow earthquakes[END_REF][START_REF] Molchan | Statistical analysis of the results of earthquake prediction, based on bursts of aftershocks[END_REF][START_REF] Molchan | Immediate foreshocks: time variation of the b-value[END_REF]. Indeed, with a limited number of events, the sum of two power laws predicted by ( 6) with an increasing weight of the additive law dP (m) as the time of the mainshock is approached will be seen as a decreasing b-value when fitted with a single GR power law.

Migration of foreshocks

The last prediction discussed here resulting from the Hypothesis is that foreshocks should migrate slowly toward the mainshock. Note that the specification (1) of the ETAS model defined in section 1.1 predicts no diffusion or migration if seismicity results solely from direct triggering (first generation from mother to daughter). Technically, this results from the separability of the space and time dependence of φ Mi (t -t i , r -r i ). In the ETAS model, diffusion and migration can be shown to result from the cascade of secondary, tertiary (and so on) triggered seismicity, akin to a (continuous-time) random walk with multiple steps [Helmstetter and Sornette, 2002b], which couples the space and time dependence of the resulting global seismicity rate. This migration or antidiffusion of the seismic activity toward the mainshock is quantified by the characteristic size R of the cluster of foreshocks which is predicted to decrease before the mainshock according to [Helmstetter and Sornette, 2002b;[START_REF] Helmstetter | Mainshocks are aftershocks of conditional foreshocks: how do foreshock statistical properties emerge from aftershock laws[END_REF]]

R ∼ (t c -t) H , (8)
with H = θ/µ for µ < 2 where θ and µ are defined in section 1.1. It is natural that the (sub-)diffusive exponent H combines the exponent θ (respectively µ) of the time-(resp. space-) dependent local processes (3) and ( 4). This law (8) describes the localization of the seismicity as the mainshock approaches, which is also observed in real seismicity [START_REF] Kagan | Statistical study of the occurrence of shallow earthquakes[END_REF][START_REF] Von Seggern | Seismicity parameters preceding moderate to major earthquakes[END_REF].

We use a superposed epoch analysis and stack all sequences of foreshocks of type II synchronized at the time of the mainshock and with a common origin of space at the location of each mainshock. The analysis of the California seismicity presented in inset of Figure 4 shows clearly a migration of the seismicity toward the mainshock, confirmed by the significant diffusion exponent H = 0.3 ± 0.1. This value is compatible with the estimates θ ∼ = 0.2 and µ ∼ = 1. We obtain the same pattern for type I and type II foreshocks, but the figures for type I foreshocks have more noise due to the smaller number of events. However, this migration is likely to be an artifact of the background activity, which dominates the catalog at long times and distances from the mainshocks. Indeed, the shift in time from the dominance of the background activity at large times before the mainshock to that of the foreshock activity clustered around the mainshock at times just before it may be taken as an apparent inverse diffusion of the seismicity rate when using standard quantifiers of diffusion processes (see [Helmstetter and Sornette, 2002b] for a discussion of a similar effect for the apparent diffusion of aftershocks).

Discussions and conclusions

By defining the foreshocks of type II and by comparing them with standard foreshocks of type I, we have revisited the phenomenology of earthquake foreshocks using the point of view of triggered seismicity formulated in our Hypothesis. We have found that the most salient properties of foreshock sequences are explained solely by the mechanism of earthquake triggering. This validates the Hypothesis.

An important result is that the precursory modification of the seismic activity before a mainshock is independent of its magnitude, as expected by the triggering model with a constant magnitude distribution. Therefore, large earthquakes are not more predictable than smaller earthquakes on the basis of the powerlaw acceleration of the seismicity before a mainshock or by using the modification of the magnitude distribution.

All these results taken together stress the importance of the multiple cascades of earthquake triggering in order to make sense of the complex spatiotemporal seismicity. In particular, our results do not use any of the specific physical mechanisms proposed earlier to account for some of the observations analyzed here. For instance, the ETAS model is different from the receding stress shadow model [START_REF] Bowman | Accelerating seismicity and stress accumulation before large earthquakes[END_REF][START_REF] Sammis | Positive feedback, memory and the predictability of earthquakes[END_REF], from the critical earthquake model [START_REF] Bowman | An Observational test of the critical earthquake concept[END_REF][START_REF] Jaumé | Evolving towards a critical point: A review of accelerating seismic moment/energy release prior to large and great earthquakes[END_REF][START_REF] Sammis | Positive feedback, memory and the predictability of earthquakes[END_REF] and from the pre-slip model [START_REF] Dodge | Detailed observations of California foreshock sequence: implications for the earthquake initiation process[END_REF], which in addition each addresses only a specific part of the seismic phenomenology. Our demonstration and/or confirmations of (i) the increase of rate of foreshocks before mainshocks (ii) at large distances and (iii) up to decades before mainshocks, (iv) a change of the Gutenberg-Richter law from a concave to a convex shape for foreshocks, and (v) the migration of foreshocks toward mainshocks are reminiscent of, if not identical to, the precursory patterns documented in particular by the Russian [Keilis-Borok and Ma-linovskaya, 1964;[START_REF] Keilis-Borok | Earthquake prediction: State-of-the-art and emerging possibilities[END_REF] and Japanese [START_REF] Mogi | Earthquake prediction research in Japan[END_REF] schools, whose physical origin has remained elusive an/or controversial. The present work suggests that triggered seismicity is sufficient to explain them. The concepts and techniques and their variations developed here could be applied to a variety of problems, such as to determine the origin of financial crashes [Sornette et al., 2003;[START_REF] Johansen | Endogenous versus Exogenous Crashes in Financial Markets[END_REF], of major biological extinctions, of change of weather regimes and of the climate, and in tracing the source of social upheaval and wars [START_REF] Sornette | Endogeneous versus exogeneous shocks in systems with memory[END_REF]].

The cascade model described here is sufficient to explain the properties of foreshocks in time, space and magnitude. There may be however other properties of foreshocks not explained by this model, such has an unusual waveform or focal mechanism, that may provide ways to distinguish foreshocks from other earthquakes in real time [START_REF] Dodge | Detailed observations of California foreshock sequence: implications for the earthquake initiation process[END_REF][START_REF] Kilb | The initial subevent of the 1994 Northridge, California earthquake: Is earthquake size predictable?[END_REF]. For instance, [START_REF] Dodge | Detailed observations of California foreshock sequence: implications for the earthquake initiation process[END_REF] found that some foreshock sequences in California were inconsistent with static stress triggering of the mainshock, and concluded that these foreshocks were more likely a by-product of an aseismic nucleation process. Several points may reconcile these observations with the cascade model. First, there are other mechanisms such as dynamical effects not taken into account in this study that may explain earthquake triggering. Indeed, a large fraction of aftershocks are inconsistent with static triggering by the mainshock. Second, some foreshocks may be missing in the catalog used and may have triggered the mainshock. Recall that, for α < b, the mainshock is more likely to be triggered by a small, possibly undetected, earthquake [Helmstetter, 2003]. 2) before M ≥ 4.5 mainshocks as a function of the distance from the mainshock for different values of the time before the mainshock ranging from 0.01 day (black line at the top) to 10 yrs (gray line at the bottom). We use logarithmic bins for the time windows, with a bin size increasing from 0.01 day up to 10 yrs as a geometric series with multiplicative factor 3.2. The number of events in each time window increases from N = 936 for 0.01 -0.03 days up to 2096633 for 1000 -3650 days. We evaluate the seismicity rate for different distances from the mainshock by counting the number of events in each shell (r, r + ∆r). The seismicity rate is normalized by the number of mainshocks, the duration of the time window and the widths of the space window ∆r (controlling the discretization of the curves) used to estimate the seismicity rate. The inset shows the characteristic size of the cluster of foreshocks, measured by the median of the distance between all foreshock-mainshock pairs, as a function of the time before the mainshock. The solid line is a fit by a power-law R ∼ t H with H = 0.3. Due to the large space-time window T = 10 yrs and R = 500 km used to select foreshocks, a large proportion of the seismicity are background events, which induces a spurious migration of seismicity toward the mainshock (see section 2.3). Figure 6. (a) Magnitude distribution P (m) of foreshocks of type II using a space window R = 20 km, for different time windows before the mainshock, ranging from black to gray as the time t c -t from the mainshock decreases from 1 yr to 0.01 day. Note the progressive increase in the proportion of large earthquakes by comparison to the normal distribution (dashed line, b = 1). Each curve contains the same number of events. (b) difference q(t) dP (m) between the foreshock magnitude. The foreshock magnitude distribution is well fitted in the magnitude range 4 ≤ m ≤ 7 by the sum of two power-laws (6), with an exponent b ′ ≈ 0.5 independently of the time from the mainshock. The amplitude q(t) of the perturbation is shown in panel (c) and the exponent b ′ of dP (m) is shown in panel (d). q(t) is fitted with a power law fit as a function of t c -t with exponent 0.3 ± 0.2 shown as the dashed line in panel (c).

Figure 1 .

 1 Figure1. Cumulative number of foreshocks of type II (see definition in section 1.2) for all M ≥ 6.5 mainshocks (thin gray lines), and average number of foreshocks per mainshock (heavy black line), obtained by stacking all (3700) foreshock sequences of M ≥ 4 mainshocks in the SCEC catalog. The foreshocks have been selected in a time-space window with T = 200 days and R = 30 km.

Figure 2 .Figure 3 .

 23 Figure2. Rate of seismic activity per mainshock for foreshocks of type II (continuous lines) defined in section 1.2 and for aftershocks (dashed lines) measured as a function of the time |t -t c | from the mainshock occurring at t c , obtained by stacking many earthquake sequences for different mainshock magnitude intervals given in the legend. This figure illustrates the power-law increase of seismicity before the mainshock (inverse Omori law) and the power-law decrease after the mainshock (direct Omori law), and show how the rate of foreshocks and aftershocks changes with mainshock magnitude.

Figure 4 .

 4 Figure4. Rate of foreshocks of type II (defined in section 1.2) before M ≥ 4.5 mainshocks as a function of the distance from the mainshock for different values of the time before the mainshock ranging from 0.01 day (black line at the top) to 10 yrs (gray line at the bottom). We use logarithmic bins for the time windows, with a bin size increasing from 0.01 day up to 10 yrs as a geometric series with multiplicative factor 3.2. The number of events in each time window increases from N = 936 for 0.01 -0.03 days up to 2096633 for 1000 -3650 days. We evaluate the seismicity rate for different distances from the mainshock by counting the number of events in each shell (r, r + ∆r). The seismicity rate is normalized by the number of mainshocks, the duration of the time window and the widths of the space window ∆r (controlling the discretization of the curves) used to estimate the seismicity rate. The inset shows the characteristic size of the cluster of foreshocks, measured by the median of the distance between all foreshock-mainshock pairs, as a function of the time before the mainshock. The solid line is a fit by a power-law R ∼ t H with H = 0.3. Due to the large space-time window T = 10 yrs and R = 500 km used to select foreshocks, a large proportion of the seismicity are background events, which induces a spurious migration of seismicity toward the mainshock (see section 2.3).

Figure 5 .

 5 Figure5. Same as Figure2for foreshocks of type I (defined in section 1.2) which have been selected using a space-time window R = R 2 = 50 km and T = 1 yr. The presented data and the statistics for aftershocks are the same as in Figure2. The total number of foreshocks of type I ranges from 1050 to 5462 depending on the mainshock magnitude. The same mainshocks are used for the selection of aftershocks and of type I foreshocks.

This preprint was prepared with AGU's L A T E X macros v4. File OminvDataRev2 formattedMay 29, 2006.

Acknowledgments. We are grateful to E. Brodsky, K. Felzer, J.-R. Grasso, H. Houston, Y.Y. Kagan, G. Ouillon, J. Vidale and an anonymous referee for stimulating discussions and constructive remarks that helped improve the manuscript. We are also grateful for the earthquake catalog data made available by the Southern California Earthquake Data Center. This work was partially supported by NSF-EAR02-30429, by the Southern California Earthquake Center (SCEC) and by the James S. Mc Donnell Foundation 21st century scientist award/studying complex system.