Agnès Helmstetter

Yan Y Kagan

David D Jackson

Comparison of short-term and long-term earthquake forecast models for southern California

. The forecasted seismicity rate is the sum of a constant background seismicity and of the aftershocks of all past earthquakes. The background rate is estimated by smoothing past seismicity. Each earthquake triggers aftershocks with a rate that increases exponentially with its magnitude and decreases with time following Omori's law. We use an isotropic kernel to model the spatial distribution of aftershocks for small (m ≤ 5.5) mainshocks, and by smoothing the location of early aftershocks for larger mainshocks. The model also assumes that all earthquake magnitudes follow the Gutenberg-Richter law with a uniform b-value. We use a maximum likelihood method to estimate the model parameters and test the short-term and long-term forecasts. A retrospective test using a daily update of the forecasts between 1985/1/1 and 2004/3/10 shows that the short-term model increases the average probability of an earthquake occurrence by a factor 11.5 compared to the long-term time-independent forecast.

Introduction

The best indicator of earthquake probability is the occurrence of recent earthquakes. Several studies show that many earthquakes are triggered in part by preceding events. Aftershocks are the most obvious examples, but many large earthquakes are preceded by smaller ones. Also, many events are apparently triggered through a cascade process in which triggered quakes trigger others in turn. Typically, the seismicity rate just after and close to a large m ≥ 7 earthquake can increase by a factor 10 4 and more, and stay above the background level for several decades especially in regions of low tectonic deformation. These large fluctuations of seismicity rate provide a way to improve earthquake forecasting significantly compared to time-independent forecasts. Small earthquakes also have a significant contribution in earthquake triggering because they are much more numerous than larger ones [START_REF] Helmstetter | Is earthquake triggering driven by small earthquakes?[END_REF][START_REF] Helmstetter | Importance of small earthquakes for stress transfers and earthquake triggering[END_REF]. As a consequence, many large earthquakes are triggered by previous smaller earthquakes ("foreshocks"), which can be useful to improve the forecasting of large earthquakes [START_REF] Kagan | Statistical shortterm earthquake prediction[END_REF][START_REF] Reasenberg | Earthquake hazard after a mainshock in California[END_REF][START_REF] Kagan | Likelihood analysis of earthquake catalogues[END_REF][START_REF] Reasenberg | Foreshock occurrence before large earthquakes[END_REF].

Short-term clustering has been recognized for some time, and several quantitative models of clustering have been proposed [START_REF] Kagan | Statistical shortterm earthquake prediction[END_REF][START_REF] Ogata | Statistical models for earthquake occurrence and residual analysis for point processes[END_REF][START_REF] Reasenberg | Earthquake hazard after a mainshock in California[END_REF][START_REF] Kagan | Likelihood analysis of earthquake catalogues[END_REF][START_REF] Kagan | Probabilistic forecasting of earthquakes[END_REF]Helmstetter and Sornette, 2003a;[START_REF] Ogata | Space-time model for regional seismicity and detection of crustal stress changes[END_REF]. Short-term forecasts based on earthquake clustering have already been developed. [START_REF] Kagan | Statistical shortterm earthquake prediction[END_REF] performed retrospective tests using the California seismicity. [START_REF] Jackson | Testable earthquake forecasts for 1999[END_REF] and [START_REF] Kagan | Probabilistic forecasting of earthquakes[END_REF] calculate in real time short-term hazard estimates for the northwest and southwest Pacific regions since 1999 (see http://scec.ess.ucla.edu/∼ykagan/predictions index.html). More recently, [START_REF] Gerstenberger | Real-time forecasts of tomorrow's earthquakes in California: a new mapping tool[END_REF] developed a method to provide daily forecasts of seismic hazard in southern California. However, these forecasts have not yet been fully optimized and tested.

We have devised and implemented a new method for issuing daily earthquake forecasts for southern California. Our forecasts are updated daily and available on-line (http://moho.ess.ucla.edu/∼helmstet/forecast.html). Short-term effects may be viewed as temporary perturbations to long-term earthquake potential. This long-term forecast could be a Poissonian process or a time-dependent process including, for example, stress shadows. It can include any geologic information based on fault geometry and slip rate, as well as data from geodesy or paleoseismicity. As a first step, we have measured the long-term seismic activity using historical seismicity only. We show that this simple model performs better than a more sophisticated model which incorporates geology data and characteristic earthquakes [START_REF] Frankel | Seismic hazard maps for California, Nevada, and Western Arizona/Utah[END_REF].

We use the Epidemic Type Earthquake Sequence (ETES) model [START_REF] Kagan | Statistical shortterm earthquake prediction[END_REF][START_REF] Ogata | Statistical models for earthquake occurrence and residual analysis for point processes[END_REF] to obtain short-term earthquake forecast including foreshocks and aftershocks. This model is usually called "Epidemic Type Aftershock Sequence" (ETAS), but in addition to aftershocks this model also describes background seismicity, mainshocks and foreshocks, using the same laws for all earthquakes. We use a maximum likelihood approach to estimate the parameters by maximizing the forecasting skills of the model [START_REF] Kagan | Statistical shortterm earthquake prediction[END_REF][START_REF] Schorlemmer | Earthquake likelihood model testing[END_REF].

Time-independent forecasts

Definition of the model

We have developed a method to estimate the longterm probability of an earthquake as a function of space and magnitude, from historical seismicity. We estimate the density of seismicity µ( r ) by declustering and smoothing past seismicity. We use the composite seismicity catalog from the Advanced National Seismic System (ANSS), available at http://quake.geo.berkeley.edu/anss/catalog-search.html. We selected earthquakes above m 0 = 3 for the time period and above m 0 = 2 since 1980. We computed the background density of earthquakes on a grid which covers southern California with a resolution of 0.05 • × 0.05 • . The boundary of the grid is [32.45 • N -36.65 • N ] in latitude and [121.55 • W -114.45 • W ] in longitude. Before computing the density of seismicity, we need to decluster the catalog, to remove the largest clusters of seismicity, which would give large peaks of seismicity which do not represent the long-term average. We used [START_REF] Reasenberg | Second-order moment of central California seismicity, 1969-82[END_REF] declustering algorithm with parameters r fact = 20, x meff = 2.00, p 1 = 0.99, τ min = 1.0 day, τ max = 10.0 days, and with a minimum cluster size of 5 events.

The declustered catalog is shown in Figure 1. Note that there is a better method of declustering, which do not need to specify a space-time window to define aftershocks, and uses an ETES-type model to estimate the probability that each earthquake is an aftershock [START_REF] Kagan | Statistical search for non-random features of the seismicity of strong earthquakes[END_REF][START_REF] Kagan | Universality of the seismic moment-frequency relation[END_REF][START_REF] Zhuang | Analyzing earthquake clustering features by using stochastic reconstruction[END_REF]. This method is however more complex to use and time-consuming for a large number of earthquakes.

The background density µ( r ) at point r is given by the sum over all past earthquakes 1 ≤ i ≤ N of the kernel K d ( r ) of each earthquake

µ( r ) = i=1 K di ( r -r i ) .
(1)

We use an isotropic kernel K d (r) given by

K d (r) = C(d) (r 2 + d 2 ) 1.5 , (2)
where C(d) is a normalizing constant, and d is the bandwidth of the kernel. We use adaptive kernels with a bandwidth d which depends on the density of earthquakes, so that we have a better resolution (smaller d) where the density is higher, but a larger d were the density of seismicity is smaller. Specifically, the bandwidth of the kernel for earthquake i is given by [START_REF] Izenman | Recent developments in nonparametric density estimation[END_REF]

d i = d 0 µ( r i ) and d i > d min . (3)
The background density defined by ( 1) and (3) has spatial variations at scales smaller than the resolution of the grid (0.05 • ) used for the forecasts. Therefore, we need to integrate the value of µ( r ) defined by (1) over each cell to obtain the background rate in this cell. The advantage of the function ( 2) is that we can compute analytically the integral of K d (x, y) over one dimension x or y, and then compute numerically the integral in the other dimension.

The function (3) has two adjustable parameters d 0 and the minimum value d min . We estimate these parameters by optimizing the likelihood of the model. We use the data from 1932 to 1995 to compute the density µ(r) on each cell, and the data since 1996 to evaluate the model. We estimate d 0 and d min by maximizing the log likelihood of the model given by the sum over all cells

LL = Nx ix=1 Ny iy=1 log p µ (i x , i y , n) , ( 4
)
where n is the number of events that occurred in the cell (i x , i y ), and the probability p µ (i x , i y , n) of having n events in the cell (i x , i y ) is given by (assuming a Poisson process)

p µ (i x , i y , n) = [T µ(i x , i x )] n exp[-µ (i x , i y )T ] n! , (5)
where µ(i x , i y ) represents the density of seismicity in the cell i x , i y per unit time and T is the duration of the catalog used to test the model. For each value of the parameters d 0 and d min we compute iteratively µ and d using expressions (1), ( 2) and (3), starting with d = d min for all earthquakes. After each iteration we get a smoother density µ. We stop the iteration when the likelihood LL of the model starts to decrease. The optimization gives d 0 = 0.0225 and d min = 0.5 km (minimum allowed value equal to location accuracy). We used 3 iterations to compute iteratively µ and d before the likelihood decreased.

The background density, estimated using the above parameters and the declustered catalog up to 2003, is shown in Figure 2.

Comparison with other long-term models

We have compared this model with other long-term forecasts for southern California, the model of [START_REF] Kagan | Long-term probabilistic forecasting of earthquakes[END_REF] (see also http://moho.ess.ucla.edu/-∼kagan/s cal tbl new.dat) and the model of [START_REF] Frankel | Seismic hazard maps for California, Nevada, and Western Arizona/Utah[END_REF]. [START_REF] Kagan | Long-term probabilistic forecasting of earthquakes[END_REF][START_REF] Wiemer | Introducing probabilistic aftershock hazard mapping[END_REF] forecasts . [START_REF] Kagan | Long-term probabilistic forecasting of earthquakes[END_REF][START_REF] Wiemer | Introducing probabilistic aftershock hazard mapping[END_REF] (KJ94) used a similar smoothing method to estimate the long-term seismicity rate in southern California (see also Kagan et al. [2003b]). The main differences between their algorithm and the present work are that • KJ94 use larger m ≥ 5.5 earthquakes since 1850, without declustering the catalog.

• KJ94 introduce a weight proportional to the logarithm of the moment of each earthquake, while we use the same weight for all earthquakes.

• KJ94 use an anisotropic power-law kernel (maximum density in the direction of the earthquake rupture), with a slower decay with distance K(r) ≈ 1/(r + R min ). Because the integral of K( r ) over the space does not converge, they truncate the kernel function after a distance R max = 200 km. Their kernel has a larger bandwidth R min = 5 km than the present model (smoother density).

• KJ94 add a constant uniform value to take into account "surprises": earthquakes that occur were no past earthquake occurred.

• KJ94 use extended sources instead of a point source: each m ≥ 6.5 earthquake is replaced by the sum of smaller scale ruptures (see http://moho.ess.ucla.edu/∼kagan/cal fps2.txt).

We modified our model to compare with KJ94 model, in order to use the same grid with the same resolution of 0.1 • . Both models were developed by using only data before 1990 to estimate the parameters and the long-term density of seismicity µ. We estimated the parameters d 0 and d min (defined in (3)) of our long-term model by using the data until 1986/1/1 to compute µ and the data from 1986/1/1 to 1990/1/1 to estimate the likelihood LL of the model. The optimization of the LL gave the parameters values d 0 = 0.02 and d min = 0.5 km. We then used those parameters and the data until 1990/1/1 to estimate the long-term density µ( r).

We use the log likelihood defined in (4) to compare KJ94 model with our long-term model. Because we want to test only the spatial distribution of earthquakes, not the predicted total number, we normalized both models by the observed number of earthquakes (N = 56). We obtain LL = -463 for KJ94 model and LL = -389 for our model. Both models are shown in Figure 3 together with the observed m ≥ 5 earthquakes since 1990. The present work thus improves the prediction of KJ94 by a factor (ratio of probabilities) exp((-389+463)/56) = 3.7 per earthquake, despite being much simpler (isotropic and point-source model). This suggests that including small earthquakes (m ≥ 2) to predict larger ones (m ≥ 5) considerably improves the predictions, because large earthquakes occur generally at the same place as smaller ones [START_REF] Kafka | Does the spatial distribution of smaller earthquakes delineate areas where larger earthquakes are likely to occur?[END_REF]. [START_REF] Frankel | Seismic hazard maps for California, Nevada, and Western Arizona/Utah[END_REF] forecasts. [START_REF] Frankel | Seismic hazard maps for California, Nevada, and Western Arizona/Utah[END_REF] (F97) model is a more complex model which includes both a smoothed historical seismicity (using m ≥ 4 earthquakes since 1933 and m ≥ 6 since 1850) and characteristic earthquakes on known faults, with a seismicity rate constrained by the geologic slip rate, and with a rupture length controlled by the fault length. The magnitude distribution follows the Gutenberg-Richter (GR) law with b = 0.9 for small magnitudes (m ≤ 6) and a bump for m > 6.2 due to characteristic events. We adjusted our model to use only data before 1996 to build the model and the same grid as F97 with a resolution of 0.1 • . We assumed a GR distribution with b = 1 and with an upper magnitude cut-off at m = 8. We used the ANSS catalog for the time period 1932-1995 to estimate the long-term rate of m ≥ 4 earthquakes (without declustering the catalog). We then assumed a GR law with b = 1 to estimate the long-term rate of m ≥ 5 earthquakes from the number of m ≥ 4 events. We use m ≥ 5 earthquakes in the ANSS catalog that occurred since 1996 to compare the models. Both models are illustrated on Figure 4, which represents the density of seismicity above m = 5.

We test how each model explains the number of observed events, as well as their location and magnitude, by comparing the likelihood of each model. The log likelihood is defined by

LL = Nx ix=1 Ny iy =1 Nm im=1 log p µ (i x , i y , i m , n) , ( 6
)
where n is the number of events that occurred in the cell (i x , i y ) and in the magnitude bin i m (the magnitude range [5.0 -8.0] is divided in bins of 0.1 unit). The log likelihood is LL = -155 for our model and LL = -161 for F97 model. We thus obtain a probability gain of 1.5 per earthquake for our model compared to F97. Our model thus better predicts the observed earthquake occurrence since 1996 than F97 model. F97 however better predicts the observed number than our model, because the number of m ≥ 5 earthquakes in the time period 1996-2004 was smaller than the long-term rate (predicted number N = 14.6 for F97 and N = 26.6 for our model, compared to the observed number N = 15). The difference in likelihood between the two models in mainly due to the choice of the kernel and of the minimum magnitude used to estimate the rate of seismicity from historical seismicity. F97 use a smoother kernel, with a fixed characteristic smoothing distance of 10 km and with a ∼ 1/r decay, and use only m ≥ 4 earthquakes for smoothing.

Focal mechanisms

In addition to the location and magnitude of earthquakes, seismic hazard estimation also requires the specification of earthquake focal mechanisms. [START_REF] Kagan | Long-term probabilistic forecasting of earthquakes[END_REF] have developed a method to estimate focal mechanisms based on catalogs of focal mechanisms (see http://moho.ess.ucla.edu/∼kagan/s cal tbl new.dat for the forecasts for southern California). They used the same catalog (m ≥ 5) and the same kernel to estimate the rate of seismicity and the average focal mechanisms in each cell (see section ). We would like to redo this analysis by applying the same smoothing kernel than for estimating the rate, including smaller m ≥ 2 earthquakes, and optimizing the forecast using a maximum likelihood method. This analysis will be presented in another paper.

Time dependent forecasts

Definition of the ETES model

The ETES model is based on two empirical laws of seismicity, which can also be reproduced by a multitude of physical mechanisms, the Gutenberg-Richter law to model the magnitude distribution, and Omori's law to characterize the decay of triggered seismicity with the time since the mainshock [START_REF] Kagan | Statistical shortterm earthquake prediction[END_REF][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]Helmstetter and Sornette, 2003a;[START_REF] Rhoades | Long-range earthquake forecasting with every earthquake a precursor according to scale[END_REF]. This model assumes that all earthquakes may be simultaneously mainshocks, aftershocks and possibly foreshocks. Each earthquake triggers direct aftershocks with a rate increasing exponentially ∼ 10 αm with the earthquake magnitude m, and that decays with time according to Omori's law ∼ 1/(t + c) p . We also assume that all earthquakes have the same magnitude distribution, which is independent of the past seismicity. Each earthquake has thus a finite probability of triggering a larger earthquake. An observed "aftershock" sequence in the ETES model is the sum of a cascade of events in which each event can trigger more events.

The global seismicity rate λ(t, r, m) is the sum of a background rate µ( r ), usually taken as a spatially non-homogeneous Poisson process, and the sum of aftershocks of all past earthquakes

λ(t, r, m) = µ( r) P m (m)+P m (m) ti<t φ mi ( r-r i , t-t i ) , (7)
where P m (m) is a time-independent magnitude distribution. The function φ m ( r, t) gives the spatiotemporal distribution of triggered events at point r and at time t after an earthquake of magnitude m

φ m ( r, t) = ρ(m) ψ(t) f ( r, m) , (8)
where ρ(m) is the average number of earthquakes triggered directly by an earthquake of magnitude m > m 0

ρ(m) = K 10 α(m-m0) , (9)
the function ψ(t) is Omori's law normalized to 1

( ∞ 0 ψ(t)dt = 1) ψ(t) = (p -1) c p-1 (t + c) p , (10)
and f ( r, m) is the normalized distribution of horizontal distances r relative to the mainshock of magnitude m. We have tested different choices for f ( r, m), which are described in section .

The exponent α has been found equal or close to 1.0 for the southern California seismicity [START_REF] Felzer | A common origin for aftershocks, foreshocks, and multiplets[END_REF][START_REF] Helmstetter | Importance of small earthquakes for stress transfers and earthquake triggering[END_REF], equal to the GR bvalue, showing that small earthquakes are collectively as important as larger ones for seismicity triggering. Note that we consider in the sum in (7) only earthquakes above the detection magnitude m 0 . Smaller undetected earthquakes may also have an important contribution to the rate of triggered seismicity. These undetected earthquakes may thus bias the parameters of the model, i.e., the parameters estimated by optimizing the likelihood of the models are "effective parameters", which more or less account for the influence of undetected small earthquakes with m < m 0 .

The ETES model assumes that each primary aftershock may trigger its own aftershocks (secondary events). The secondary aftershocks may themselves trigger tertiary aftershocks and so on, creating a cascade process. The exponent p, which describes the time distribution of direct aftershocks, is larger than the observed Omori exponent, which characterizes the whole cascade of direct and secondary aftershocks [START_REF] Helmstetter | Sub-critical and super-critical regimes in epidemic models of earthquake aftershocks[END_REF].

As a first step, we use a simple GR magnitude distribution

P m (m) ∼ 10 -b(m-m0) , (11)
with a uniform b-value equal to 1.0, a upper cut-off at m max = 8 [START_REF] Bird | Plate-tectonic analysis of shallow seismicity: apparent boundary width, beta, corner magnitude, coupled lithosphere thickness, and coupling in seven tectonic settings[END_REF], and a minimum magnitude m 0 = 2. If we want to predict relatively small m ≤ 4 earthquakes, we must take into account the fact that small earthquakes are missing in the catalog after a large mainshock. The procedure described for correcting for undetected small earthquakes is described in section .

The background rate µ in ( 7) is estimated by smoothing the historical seismicity, using the longterm forecast described by equations ( 1) and (2) in section . Alternatively, our short-term forecast model can be included in a more complex interdisciplinary model, which can include information from geodesy, including GPS, paleo-seismicity and fault slip rates.

We can build a composite model by adding the shortterm component to any long-term forecast. For this, we simply need to replace the background rate µ in ( 7) by the long-term forecast, multiplied by the fraction of time-independent events so that the total predicted seismicity rate (time-independent and triggered events) will be equal to the observed seismicity rate.

We use the QDDS (ftp://clover.wr.usgs.gov/pub/-QDDS/QDDS.html) and the QDM (ftp://clover.wr.usgs.gov/pub/QDM/QDM.html) java applications to obtain the data (time, locations and magnitude) in real time from several regional networks (southern and northern California, Nevada) and to create a composite catalog. We automatically update our forecast each day and we put our results on our web site (http://moho.ess.ucla.edu/∼helmstet/forecast.html). The model parameters are estimated by optimizing the prediction (maximizing the likelihood of the model) using retrospective tests. The inversion method and the results are presented in section .

The basic idea of our work is the same as the STEP model of [START_REF] Gerstenberger | Real-time forecasts of tomorrow's earthquakes in California: a new mapping tool[END_REF]. Both models are based on spatio-temporal clustering of earthquakes, but there are several significant differences between the two approaches:

• STEP model uses [START_REF] Frankel | Seismic hazard maps for California, Nevada, and Western Arizona/Utah[END_REF] forecasts for the background rate, while we estimate the background rate by declustering and smoothing the seismicity;

• In STEP model, the seismicity rate on each cell is the sum of the background rate plus the rate of aftershocks of only one mainshock, the previous earthquake which has the largest influence on this cell. If the rate of aftershocks of the mainshock is smaller than the background rate, then the seismicity rate is equal to the background rate. In our model, the seismicity rate on each cell is given by the sum of the background rate and of the aftershock rate of all past earthquakes.

• STEP model assumes than the rate of aftershocks scales as ∼ 10 αm with α = b, while we use α < b;

• STEP model (in their more complex version) adjusts the parameters p, K (aftershocks productivity) and c of Omori's law for each aftershock sequence and for each cell (if the number of aftershocks is sufficient), using a fit by a maximum likelihood method. If the number of events is too small to adjust the parameters, the generic aftershock model of [START_REF] Reasenberg | Earthquake hazard after a mainshock in California[END_REF] is used. In our model, we use the same parameters p, K and c for all earthquakes, which are estimated by optimizing the model using the whole catalog.

• STEP model introduces a variability of the GR b-value in time and in space, while we use a constant value b = 1.0. Bird and [START_REF] Kagan | Short-term properties of earthquake catalogs and models of earthquake source[END_REF] argue that the b-value has a uniform value 0.9 -1.0 for global M ≥ 5.6 earthquakes in various tectonic provinces.

We will compare our model with STEP and other time-dependent models as part of the Regional Earthquake Likelihood Model (RELM) project to test in real time daily and long-term models for California [Kagan et al. 2003a;[START_REF] Jackson | Prospective Tests of Southern California Earthquake Forecasts[END_REF][START_REF] Schorlemmer | Earthquake likelihood model testing[END_REF].

Application of ETES model for time-dependent forecasts

By definition, the ETES model provides the average instantaneous seismicity rate λ(t) at time t given by ( 7), if we know all earthquakes that occurred up to time t. To forecast the seismicity between the present time t p and a future time t p + T , we cannot use directly expression (7), because a significant fraction of earthquakes that will occur between time t p and time t p + T will be triggered by earthquakes that will occur between time t p and time t p + T (see Figure 5). Therefore using expression (7) to provide short-term seismicity forecasts, with a time window T of 1 day, may significantly underestimate the number of earthquakes. Helmstetter and Sornette [2003a] have tested several methods of prediction on synthetic catalogs generated by the ETES model. They have shown that the use of ( 7) to predict the number of earthquakes in the next day estimated by N p = tp+T tp λ(t) dt, using the exact parameters of the ETES model to compute λ(t), may underestimate the number of earthquakes by a factor 2, because it does not take into account the contribution of yet unobserved seismicity in the seismicity rate λ(t).

To solve this problem, Helmstetter and Sornette [2003a] proposed to generate synthetic catalogs with the ETES model in order to predict the seismicity for the next day, by averaging the number of earthquakes over all scenarios. This methods provides a much better estimation of the number of earthquakes than the direct use of ( 7), but is much more complex to apply. It is practically not useful if we want to predict not only the number of earthquakes, but also their location and magnitude, because this would require to generate a huge number of scenarios. Helmstetter and Sornette [2003a] have shown that, for synthetic ETES catalogs, the use of N p = tp+T tp λ(t) dt to predict the number of earthquakes between t p and t p +T underestimates the number of actually occurred earthquakes by an approximate constant factor, independent of the future events number. This means that the effect of yet unobserved seismicity is to amplify the aftershock rate of past earthquakes by a constant factor.

This result suggests a simple solution to take into account the effect of yet unobserved earthquakes. We can use N p = tp+T tp λ(t) dt to predict the number of earthquakes between t p and t p + T but using effective parameters K, µ s and α, which may be different from the true (i.e., t p = t) parameters of the ETES model. Instead of using the likelihood of the ETES model to estimate these parameters, as done by [START_REF] Kagan | Likelihood analysis of earthquake catalogues[END_REF], we will estimate the parameters of the model by optimizing the likelihood of the forecasts, defined in section . These effective parameters will depend on the duration (horizon) T of the forecasts.

Threshold Magnitude

An important problem when modeling the occurrence of relatively small earthquakes m ≤ 4 in California is that the catalog is not complete after large earthquakes [START_REF] Kagan | Short-term properties of earthquake catalogs and models of earthquake source[END_REF]. One effect of missing earthquakes is that the model will over-estimate the observed number of earthquakes, because small earthquakes are not detected. But another effect of missing early aftershocks is to underestimate the predicted seismicity rate, because we miss the contribution from these undetected small earthquakes in the future seismicity rate estimated from the ETES model ( 7). Indeed, secondary aftershocks (triggered by a previous aftershock) represent an important fraction of aftershocks [START_REF] Felzer | Secondary aftershocks and their importance for aftershock forecasting[END_REF]Helmstetter and Sornette, 2003b].

We have developed a method to correct from both effects of undetected small aftershocks. We first need to estimate the threshold magnitude as a function of the time from the mainshock and of the mainshock magnitude. We analyzed all aftershock sequences of m ≥ 6 earthquakes in southern California since 1985. We propose the following relation between the threshold magnitude m c (t, m) at time t (in days) after a mainshock of magnitude m m c (t, m) = m -4.5 -0.76 log 10 (t)

and m c (t, m) ≥ m 0 = 2 . (12)
Of course, there are some fluctuations between one sequence and another one, but the above relation is correct within ≈ 0.2 magnitude units. The above relation is illustrated on Figure 6 for Joshua-Tree m = 6.1, San-Simeon m = 6.5 and Landers m = 7.3 aftershock sequences.

We use expression ( 12) to estimate the detection magnitude magnitude m c (t) at the time t of each earthquake. The time-dependent detection threshold m c (t) is larger than the usual threshold m 0 for earthquakes that occurred at short times after a large m ≥ 5 earthquake. We select only earthquakes with m i > m c (t i ) to estimate the likelihood of the forecasts (18).

We can also correct the forecasts for the second effect: missing contribution from undetected aftershocks in the sum (7). We can take into account the effect of earthquakes below the detection threshold m c (t) and above the minimum magnitude m 0 by adding a contribution to the number ρ(m) of aftershocks of detected earthquakes m > m c (t), i.e., by replacing ρ(m) in ( 9) by

ρ * (m) = ρ(m)+ Kb b -α 10 b(mc(t)-m0) 1 -10 -(b-α)(mc(t)-m0) .
(13) where m c (t) is the detection threshold at the time t of the earthquake, estimated by ( 12), due to the effect of all previous m ≥ 5 earthquakes. The second contribution corresponds to the effect of all earthquakes with m 0 < m < m c (t) that occur on average for each detected earthquake. Practically, for a reasonable value of α ≈ 0.8, this correction ( 13) is of the same order as the contribution from observed earthquakes, because a large fraction of aftershocks are secondary aftershocks [START_REF] Felzer | Secondary aftershocks and their importance for aftershock forecasting[END_REF], and because small earthquakes are collectively as important as larger ones for earthquake triggering if α = b.

Spatial distribution of aftershocks

We have tested different choices for the spatial kernel f ( r, m), which models the distribution of distances between a mainshock of magnitude m and its aftershocks. We used a power-law function

f pl ( r, m) = C pl [ | r | 2 + d(m) 2 ] 1.5 , (14)
and a Gaussian distribution

f gs ( r, m) = C gs exp - | r | 2 2 d(m) 2 , ( 15
)
where C pl and C gs are normalizing constants such that the integral of f ( r, m) over an infinite surface is equal to 1. The spatial regularization distance d(m) accounts for the finite rupture size and for the location errors. We assume that d(m) is given by

d(m) = 0.5 + f d × 0.01 × 10 0.5m km, ( 16
)
where the first term accounts for location accuracy and the second term represents the aftershock zone length of an earthquake of magnitude m. The parameter f d is adjusted by optimizing the prediction and should be close to 1 if the aftershock zone size is equal to the rupture length as estimated by [START_REF] Wells | New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement[END_REF].

The Gaussian kernel (15), which describes the density of earthquakes at point r, is equivalent to the Rayleigh distribution ∼ exp[-(r/d) 2 /2]/r of distances | r| used by [START_REF] Kagan | Probabilistic forecasting of earthquakes[END_REF]. The choice of an exponent 1.5 in ( 14) is motivated by recent studies [START_REF] Ogata | Space-time model for regional seismicity and detection of crustal stress changes[END_REF][START_REF] Console | Refining earthquake clustering models[END_REF][START_REF] Zhuang | Analyzing earthquake clustering features by using stochastic reconstruction[END_REF] who inverted this parameter in earthquake catalogs by maximizing the likelihood of the ETES model, and who all found an exponent close to 1.5. This choice is also convenient because the function ( 14) is integrable analytically. It predicts that the aftershock density decreases with the distance r from the mainshock as 1/r 3 in the far field, proportionally to the static stress change.

For large earthquakes, which have a rupture length larger than the grid resolution of 0.05 • (≈ 5 km) and a large number of aftershocks, we can improve the model by using a more complex anisotropic kernel, as done previously by [START_REF] Wiemer | Spatial variability of seismicity parameters in aftershock zones[END_REF], [START_REF] Wiemer | Introducing probabilistic aftershock hazard mapping[END_REF] and [START_REF] Gerstenberger | Real-time forecasts of tomorrow's earthquakes in California: a new mapping tool[END_REF]. We use the location of early aftershocks as a witness for estimating the mainshock fault plane, and the other active faults in the vicinity of the mainshock. We compute the distribution of aftershocks of large m ≥ 5.5 mainshocks by smoothing the location of early aftershocks

f ( r, m) = N i=1 f (| r -r i |, m i ) , (17)
where the sum is on the mainshock and on all earthquakes that occurred within a distance D aft (m) from the mainshock, and at a time smaller than the present time t p and not larger than T aft from the mainshock. We took D aft (m) = 0.02 × 10 0.5m km (approximately 2 rupture lengths) and T aft = 2 days. The kernel f (r, m) in ( 17) used to smooth the location of early aftershocks, is either a power-law ( 14) or a Gaussian distribution (15), with an aftershock zone length given by ( 16) for the mainshock, but fixed to d = 2 km for the aftershocks. This way we recover the same expression for f ( r, m) for a large m ≥ 5.5 earthquake as for a smaller earthquake if we compute this term just after the mainshock has occurred (N = 1 in ( 17)). At larger times after the mainshock, when a sufficient number of aftershocks have already occurred, f ( r, m) is given by the smoothed density of aftershocks. The density of aftershocks estimated using ( 17) is shown in Figure 7 for the Landers earthquake, using a power-law kernel (a) or a Gaussian kernel (b). The distribution of aftershocks that occurred after more than two hours after Landers (black dots) is in good agreement with the prediction based on aftershocks that occurred in the first two hours (white circles). In particular, the largest Big-Bear aftershock (m = 6.4, latitude=34.2 • , longitude =-116.8 • ), which occurred about three hours after Landers, was preceded by other earthquakes in the first two hours after Landers, and is well predicted by our method. The Gaussian kernel (15) produces a density of aftershocks which is much more localized than with a power-law kernel.

The advantage of using the observed aftershocks to predict the spatial distribution of future aftershocks is that this method is completely automatic, fast, and uses only information from the time and location of aftershocks which are available a few minutes after the earthquake. It can provide an accurate prediction of the spatial distribution of future aftershocks after less than one hour after the mainshock when enough aftershocks have occurred. Our method has also the advantage of taking into account the geometry of the active fault network close to the mainshock, which is reflected by the spatial distribution of aftershocks.

Therefore, even if the spatial distribution of aftershock is controlled by the Coulomb stress change, it may be more accurate, much simpler and faster to use the method described above rather than to compute the Coulomb stress change. Indeed, Coulomb stress change calculation requires to know the mainshock fault plane and the slip distribution, which are available only several hours or days after a large earthquake [START_REF] Scotti | Coulomb stress modelling as a practical tool in real-time aftershock hazard assessment: the example of the PRESAP blind test[END_REF][START_REF] Steacy | Sensitivity of static stress calculations to the earthquake slip distribution[END_REF]. [START_REF] Felzer | Secondary aftershocks and their importance for aftershock forecasting[END_REF] have already shown that a simple fore-casting model (simplified ETES model), based on the time, location, and magnitudes of all previous aftershocks better predicts the location of future aftershocks than Coulomb stress change calculations.

Definition of the likelihood and estimation of the model parameters

We use a maximum likelihood method to test the forecasts and to estimate the parameters. For the temporal distribution of aftershocks, we choose to fix the c-value in Omori's law (10) equal to 0.0035 day (5 minutes). This parameter is not important as long as it is much smaller than the time window T = 1 day of the forecast. We estimate the parameter p of Omori's law, the exponent α and the constant K of the aftershock productivity ρ(m) defined by ( 9). We write the background rate µ( r) as the product µ s µ 0 ( r), where µ 0 ( r) is the background density normalized to 1 (see section ) and µ s gives the total number of background events per day. We have thus 5 parameters to estimate: p (Omori exponent defined in (10)), K and α (see eq. ( 9)), µ s , and f d (parameter defined by ( 16), which describes the size of the aftershock zone).

The log likelihood (LL) of the forecasts is defined by [START_REF] Kagan | Probabilistic forecasting of earthquakes[END_REF], Kagan et al. [2003b] where n is the number of events that occurred in the bin (i t , i x , i y , i m ). We take a step of T = 1 day in time, 0.05 degree in space and 0.1 in magnitude. The forecasts are updated each day at midnight Los Angeles time. We assume a Poisson process to estimate the probability p Np (i t , i x , i y , i m , n) of having exactly n events in each bin from the predicted number N p (i t , i x , i y , i m ) of events in this bin:

p λ (i t , i x , i y , i m , n) = (N p (i t , i x , i y , i m )) n exp[-N p (i t , i x , i y , i m )] n! , (19
) where N p p(i t , i x , i y , i m ) is given by the integral over each space-time-magnitude bin of the predicted seismicity rate λ( r, t, m) N p (i t , i x , i y , i m ) = We can simplify the expression of LL, by noting that we need to compute the seismicity rate only in the bins (i x , i y , i m ) that have a non-zero number of observed events n. We can rewrite ( 18) and ( 19) as

LL = Nt it=1 -N p (i t ) + n>0 ix n>0 iy n>0 im n log[N p (i t , i x , i y , i m )] -log(n!) , (21)
where N p (i t ) is the total predicted number of events between t p and t

p + T N p (i t ) = µ s + it+1 it ti<tp f i ρ(m i )ψ(t p -t i ) dt . (22)
The factor f i in ( 22) is the integral of the spatial kernel f i ( r -r i ) over the grid, which is smaller than 1 due to the finite size of the grid.

We maximize the log likelihood LL defined by (18) using a simplex algorithm [Press et al., 1992, p. 402], and using all earthquakes with m ≥ 2 since 1985/1/1 and until 2004/3/10 to test the forecasts. We take into account in the seismicity rate (7) the aftershocks of all earthquakes with m ≥ 2 since 1980/1/1 that occurred within the grid ([32.45 • N -36.65 • N ] in latitude and [121.55 • W -114.45 • W ] in longitude) or at less than 1 • outside the grid. There are 65,664 target earthquakes above the threshold magnitude m c in the time and space window used to compute the LL. We test different models for the spatial distribution of aftershocks, a power-law kernel (14) or a Gaussian (15).

We use the probability gain per earthquake G to quantify the performance of the short-term prediction by comparison to the long-term forecasts

G = exp LL -LL LT N , (23)
where LL LT is the log-likelihood of the long-term model, LL is the likelihood of the ETES model and N is the total number of events. The gain defined by ( 23) is related to the information per earthquake I defined by [START_REF] Kagan | Earthquake risk prediction as a stochastic process[END_REF] (see also [START_REF] Daley | Scoring probability forecasts for point processes: The entropy score and information gain[END_REF]) by G = 2 I . A certain caution is needed in interpreting the probability gain for the ETES model. Earthquake temporal occurrence is controlled by Omori's law which diverges to infinity for time approaching zero. Calculating the likelihood function for aftershock sequences illustrates this point: the rate of aftershock occurrence after a strong earthquake increases by a factor of thousands. Since log(1000) = 6.9, one early aftershock yields a contribution to the likelihood function analogous to about 7 additional free parameters. This means that the likelihood optimization procedure as well as the probability gain value is strongly dependent on early aftershocks. As Figure 6 demonstrates, many early aftershocks are missing from earthquake catalogs [START_REF] Kagan | Short-term properties of earthquake catalogs and models of earthquake source[END_REF], therefore the likelihood function value substantially depends on poor quality data in the beginning of aftershock sequence.

Similarly, earthquake hypocenters are concentrated on a fractal set with a correlation dimension slightly above 2.0 [START_REF] Helmstetter | Importance of small earthquakes for stress transfers and earthquake triggering[END_REF]. Due to random location errors for small inter-earthquake distances the dimension increases close to 3.0. This signifies that the likelihood function would substantially depend on location uncertainty, since kernel width K d ( r) (Equations 1 and 2) can be made smaller if a catalog with higher location accuracy is used.

Results and discussion

Model parameters and likelihood

The model parameters are obtained by maximizing the LL. The optimization converges after about 100 iterations. We have tried different values of the initial parameters and we have checked that the final values do not depend on the initial values. The results are given in Table 1 and in Figure 8. We have tested different versions of the model (spatial kernel, unconstrained or fixed α value, and different values of the minimum magnitude).

An example of our daily forecasts (using model 3 in Table 1) is shown in Figure 9, for the day of 2004/10/23. All 6 earthquakes which occurred during that day are located in areas of high predicted seismicity rate (large values of N p ). All except one occurred close enough in time and space from a recent earthquake, so that the short-term predicted number N p ( r ) is larger than the long-term average rate µ LT ( r ). The probability gain (23) per earthquake for this day is 26.

Figure 8 shows the LL of the daily forecasts, for the time period from 1985/1/1 to 2004/3/10, and for each iteration of the optimization as a function of the model parameters. The variation of the LL with each model parameter gives an idea of the resolution of this parameter. The unconstrained inversion gives a probability gain G = 11.7, and an exponent α = 0.43, much smaller than the direct estimation α = 0.8 ± 0.1 [START_REF] Helmstetter | Is earthquake triggering driven by small earthquakes?[END_REF] or α = 1 [START_REF] Felzer | A common origin for aftershocks, foreshocks, and multiplets[END_REF][START_REF] Helmstetter | Importance of small earthquakes for stress transfers and earthquake triggering[END_REF] using a fit of the number of aftershocks as a function of the mainshock magnitude. The optimization with α fixed to 0.8, closer to the observed value, provides a probability gain G = 11.1 slightly smaller than the best model. Note that there is a negative correlation between the parameters K and α in Table 1: K is larger for a smaller α in order to keep the number of forecasted earthquakes constant.

Comparison of predicted and observed aftershock rate

Figure 10 compares the predicted number of events following Landers mainshock, for the unconstrained model 2 (see Table 1), and for models 3 and 5 with α fixed to 0.8. Model 3 underestimates the number of aftershocks, but predicts the correct variation of the seismicity rate with time. In contrast, model 2 (with α = 0.43) highly underestimates the number of aftershocks until 10 days after Landers, because the low value of α yields a relatively small increase of seismicity at the time of the mainshock. Model 2 then provides a good fit to the end of the aftershock sequence, when enough aftershocks have occurred so that the predicted seismicity rate increases due to the importance of secondary aftershocks. The saturation of the number of aftershocks at early times in Figure 10 (for both the model and the data) is due to the increase of the threshold magnitude m c (see equation ( 12)), which recovers the usual value m c = m 0 = 2 about 10 days after Landers. Adding the corrective term ρ * (m) defined by ( 13), in order to account for the contribution of undetected early aftershocks in the rate of triggered seismicity, better predicts the rate of aftershocks just after Landers, but gives on average a smaller probability gain than without including this corrective term (see models 3 and 5 in Table 1 and Figure 10).

Figure 11 shows the predicted number of earthquakes and the probability gain (see equation ( 23)) in the time window 1992-1995, for model 3. The model underestimates the rate of aftershocks for Joshua-Tree (m = 6.1) and Landers (m = 7.3) mainshocks, slightly over-estimates for Northridge (m = 6.6) and provides a good fit (not shown) for Hector-Mine (m = 7.1) and for San-Simeon (m = 6.5). All models overestimate by a factor larger than two the aftershock productivity of the 1987 m = 6.6 Superstition Hills earthquake. This shows that there is a variability of aftershock productivity not taken into account by the model, which may partly be due to error in magnitudes. This implies that a model which estimates the parameters of each aftershock sequences (aftershocks productivity, Omori p exponent and GR bvalue), such as the STEP model [Gerstenberger et al., 2004] may perform better than the ETES model which uses the same parameters for all earthquakes (except for the increase in productivity ρ(m) with magnitude).

Figure 12 shows the predicted number of m ≥ 2 earthquakes per day for model 3 (see Table 1) as a function of the observed number. Most points in this plot are close to the diagonal, i.e., the model generally gives a good prediction of the number of events per day. A few points however have a large observed number of earthquakes but a small value of the predicted number. These points correspond to days for which a large earthquake and its first aftershocks occurred while the seismicity was close to its background level, when the predicted seismicity rate was small.

We can complexify the model to take into account fluctuations of aftershock productivity, as done in the STEP model, by using early aftershocks to estimate the productivity ρ(m) of large earthquakes. However, it needs to be investigated whether magnitude errors, biases and systematic effects significantly contribute to prediction efficiency. A method which adjust parameters to available data may seemingly perform better, especially in retrospective testing when various adjustments are possible. But if aftershock rate fluctuations are being caused by various technical factors and biases, this forecast advantage can be spurious.

Proportion of aftershocks in seismicity

The background seismicity is estimated to be µ s = 2.81 m ≥ 2 earthquakes per day for model 3, compared to the long-term average rate µ LT = 9.4 earthquakes, i.e., the proportion of triggered earthquakes is 70%. This number under-estimates the actual fraction of triggered earthquakes, because it does not count the early aftershocks that occur a few hours after a mainshock, between the present time t p and the end of the prediction window t p + T (see Figure 6). We have also removed from the catalog aftershocks smaller than the threshold magnitude m c (t, m) given by ( 12).

Scaling of aftershock productivity with mainshock magnitude

There may be several reasons for the small value α = 0.43 selected by the optimization, compared to the value α = 1 estimated by [START_REF] Felzer | A common origin for aftershocks, foreshocks, and multiplets[END_REF] and [START_REF] Helmstetter | Importance of small earthquakes for stress transfers and earthquake triggering[END_REF]. A smaller α value corresponds to a weaker influence of large earthquakes. A model with a small α has thus a shorter memory in time, and can adapt faster to fluctuations of the observed seismicity. A smaller α predicts a larger proportion of secondary aftershocks after a large mainshock. Therefore, it can better account for fluctuations of aftershock productivity. Indeed, if the rate of early aftershocks is low, a model with a small α will predict a small number of future aftershocks (less secondary aftershocks). In contrast, in case of a high rate of early aftershocks, it will predict a high aftershock activity due to the effect of secondary aftershocks. A model with a smaller α is also less sensitive to error in magnitudes. An error on the mainshock magnitude of 0.3 gives an error for the rate of direct aftershocks of a factor two for α = 0.8, and a factor 1.3 for α = 0.4. Finally, a model with a smaller α may provide a better forecast for the spatial distribution of aftershocks. Because the real spatial distribution of aftershocks is significantly different from the isotropic model (used for m ≤ 5.5 earthquakes), a model with a smaller α may perform better than the model with the true α. A small α gives more importance to secondary aftershocks, and can thus better model the heterogeneity of the spatial distribution of aftershocks. In contrast, a larger α value produces a quasi-isotropic distribution at short times, dominated by the mainshock contribution.

The corrective contribution ρ * (m) > ρ(m) (13), introduced to take into account the contribution of missing aftershocks, can also bias the value of α. Using this term ρ * (m) with a value of α smaller than the true value over-estimates the contribution of small earthquakes just after a large earthquake when m c > m 0 . For this reason we did not used this contribution (except for model 5 in Table 1). We also tried to use this method (13) to take into account the fact that the minimum magnitude of triggering earthquakes m 0 is probably below the detection threshold. We found that this parameter m 0 is not constrained by maximizing the LL. The best fitting value of m 0 is equal to the detection threshold m = 2, and the LL is almost independent of m 0 . Therefore, we took m 0 equal to the detection threshold m 0 = 2.

The main interest of short-term forecasts is to pre-dict the rate of seismicity after a large mainshock, when the best model with α = 0.4 clearly underestimates the observations. Therefore, we choose to constrain the value of α = 0.8 (models 3 and 4 in Table 1). This model gives a slightly smaller likelihood than the best model, but provides a best fit just after a large mainshock. Note that this value α = 0.8 is probably still smaller than the true value; a direct estimation of this parameter from a fit of the number of aftershocks as a function of mainshock magnitude gives indeed α = 1.05±0.1 [START_REF] Helmstetter | Importance of small earthquakes for stress transfers and earthquake triggering[END_REF]. Using for the forecasts a value of α smaller than the physical value allow us to improve the performance of the model because a smaller α best accounts for fluctuations in aftershock productivity as well as problems in the earthquake catalogs (incompleteness, location and magnitude errors).

Spatial distribution of aftershocks.

The power-law kernel ( 14) gives a slightly better LL than the Gaussian kernel (15) (see Table 1) for the unconstrained models 1 and 2 (α is adjustable parameter), but the Gaussian kernel works a little better when α is fixed to 0.8 (see models 3 and 4 in Table 1). The parameter f d defined in ( 16) is the ratio of the typical aftershock zone d(m) (16) and of the mainshock rupture length L(m) = 0.01 × 10 0.5m km. For the Gaussian kernel (15) f d ≈ 1, i.e., the average distance between a mainshock and its (direct) aftershocks is close to the mainshock rupture length.

For the power-law kernel (14), the average distance is not defined. In this case, d(m) represents the scaling distance at which f pl (r) starts to decrease with r. The inversion of f d using a power-law kernel gives an unrealistically small value of f d ≤ 0.06 for model 2 (see Table 1), so that d(m) ≈ 0.5 km (fixed minimum value of d(m) equal to the location error) independently of the magnitude of the triggering earthquake for m ≤ 5. It gives short-range interactions, with most of the predicted rate concentrated in the cell of the triggering earthquake. Using a complex spatial distribution of aftershocks for m ≥ 5.5 earthquakes (obtained by smoothing the location of early aftershocks, see section ) slightly improves the LL compared to the simple isotropic kernel (see models 3 and 6 in Table 1).

Probability gain as a function of magnitude

Table 1 shows the variation of the probability gain G as a function of the lower cut-off magnitude m min .

We used m ≥ 2 earthquakes in models 7-11 to estimate the forecasted rate of m ≥ m min earthquakes, with m min ranging between 3 and 6, and using the same parameters as in model 3 (but multiplying the background rate µ s by 10 -(mmin-2.0) to estimate the background rate for m ≥ m min earthquakes). The probability gain is slightly larger for m min = 3 than for m min = 2, but then G decreases with m min for m min ≥ 4. For m min = 6 (only 8 earthquakes), the long-term model (with a rate adjusted so that it predicts the exact number of observed events) performs even better than the ETES model (G < 1) for model #10 in Table 1.

We think that this variation with m min does not mean that our model predicts only small earthquakes ("aftershocks"), or that larger earthquakes have a different distribution in space and time than smaller ones, but that these results simply reflects the large fluctuations of the probability gain from one earthquake to another one: the difference in likelihood between ETES and the long-term model is mainly due to a few large aftershock sequences. We thus need a large number of earthquakes and aftershock sequences to compare different forecasts.

Table 2 compares the predicted seismicity rate at the time and location of each m ≥ 6 earthquake, estimated for the ETES model and for the long-term model. For each earthquake, we give two values of the predicted number of earthquakes, using the same parameters of the ETES model, but changing the time at which we update the forecasts, either midnight (universal time) for model #10 (see line 10 in Table 1) or at 13:00 for model #11. The large differences in the predicted seismicity rate between these two models show that the forecasts are very sensitive to short-term clustering, which has a large influence on the predicted seismicity rate. This suggests that the number of m ≥ 6 earthquakes in the catalog (8 earthquakes from 1985 to 2004) is too small to compare our short-term and long-term models for this magnitude range.

While some of these m ≥ 6 earthquakes are preceded by a short-term (hours) increase of seismicity (Superstition-Hill, Joshua-Tree, Landers, Big-Bear, Hector-Mine), the long-term model performs better than the ETES model if the forecasts are not updated between the foreshock activity and the mainshock (e.g., with model #10, between Elmore-Ranch and Superstition-Hill, and between Landers and Big-Bear). Landers occurred about two months after Joshua-Tree, and its hypocenter was just outside Joshua-Tree aftershock zone, so that the predicted seismicity rate at the location of Landers hypocenter, and before the precursory foreshock activity (which started 6 hours before Landers) was slightly lower than the long-term rate. Joshua-Tree had foreshocks, which started two hours before the mainshock, and thus were not included in the daily forecasted rate for both ETES models. Hector-Mine was also preceded by foreshocks, with m ≤ 3.6, which started about 20 hours before the mainshock. Therefore, the predicted seismicity rate (using ETES model #10) is 120 larger than the long-term rate for Hector-Mine. Other large m ≥ 6 earthquakes (Elmore-Ranch, Northridge, San-Simeon), were not preceded by any significant foreshock activity. Therefore the forecasted seismicity was smaller than the long-term average.

Updating the forecasts more often (each hour, or after each earthquake) would of course increase the performance of our short-term forecasts. But optimizing and testing the forecasts would then be much more difficult and time consuming if the duration of the forecasts (one day) is different from the time interval between two forecasts. Moreover, preliminary earthquake catalogs are much less accurate in the first few hours, especially after a strong earthquake.

Conclusion

We have proposed a model for daily earthquake forecasts in southern California. Our model is based on empirical laws of seismicity: the Gutenberg-Richter magnitude distribution, Omori's law, and the exponential increase of triggered seismicity with the mainshock magnitude. Our model includes only data from earthquake catalogs (time, magnitude and locations).

Including small m ≥ 2 earthquakes in long-term forecasts significantly improves the predictions. Our model performs better than a more complex one which incorporates geological data, when tested on m ≥ 5 earthquakes since 1996. Our model also forecasts well the spatial distribution of future aftershocks by smoothing the locations of early aftershocks. We can obtain a good forecast of the aftershocks within a few hours of a large m ≥ 5.5 earthquake, based on plentiful early aftershocks. Even if the spatial distribution of aftershocks is controlled by Coulomb stress changes, our empirical method may be more accurate, and faster, than direct calculations of the Coulomb stress change. Our method is accurate because the distribution of early aftershocks represents well the mainshock rupture surface, and because our method accounts for secondary aftershocks.

Retrospective tests for m ≥ 2 earthquakes in the time period 1985/1/1 to 2004/3/10 show that our short-term model realizes a probability gain of 11.5 over a long-term stationary Poissonian forecast. Several features of our model could be improved. First, geologic slip rate and geodetic strain rate data could be used to better constrain the long-term seismicity. Second, a better estimate of the magnitude distribution, resulting from statistical studies of the relationship between fault geometry and earthquakes, could improve the forecasting of large quakes. Third, other research [e.g. [START_REF] Gerstenberger | Real-time forecasts of tomorrow's earthquakes in California: a new mapping tool[END_REF]] suggests that aftershock productivity and magnitude distribution may vary considerably from one sequence to another. Comparing our model with others proposed to the Regional Earthquake Likelihood Models (RELM) working group [Kagan et al. 2003a;[START_REF] Jackson | Prospective Tests of Southern California Earthquake Forecasts[END_REF][START_REF] Schorlemmer | Earthquake likelihood model testing[END_REF]] should help to improve all available models. For example, the STEP model of [START_REF] Gerstenberger | Real-time forecasts of tomorrow's earthquakes in California: a new mapping tool[END_REF] can be compared directly with our model. c Using the same isotropic kernel for all earthquakes (including m ≥ 5.5 earthquakes).

d Updating the forecasts at 13:00 instead of 0:00 (universal time). . Plot of the magnitude versus time for a few days in the ANSS catalog, which illustrates the fact that a significant fraction of earthquakes that will occur in the next day (between the present time t p and t p + T ) may be triggered by earthquakes that will occur in the next day (t p < t < t p + T ).

  Np (i t , i x , i y , i m , n) ,(18)

  , t, m) dm dt dx dy .

Figure 1 .Figure 2 .Figure 3 .Figure 4 .

 1234 Figure 1. Declustered catalog, obtained with Reasenberg's [1985] algorithm, including 6861 m ≥ 3 earthquakes in the time window 1932 -1979, and 46,937 earthquakes with m ≥ 2 for 1980 -2003.

  Figure5. Plot of the magnitude versus time for a few days in the ANSS catalog, which illustrates the fact that a significant fraction of earthquakes that will occur in the next day (between the present time t p and t p + T ) may be triggered by earthquakes that will occur in the next day (t p < t < t p + T ).

Figure 6 .Figure 7 .Figure 8 .

 678 Figure 6. Magnitude versus time since mainshock for aftershocks (a) of Joshua-Tree m = 6.1, (b) San-Simeon m = 6.5, and (c) Landers m = 7.3 earthquakes. The continuous line represents the threshold magnitude estimated from (12) and includes the effect of all m ≥ 5 earthquakes. The vertical lines in (c) are due to the increase of m c (t) after large m ≥ 5 aftershocks. Dates of these earthquakes are shown in Table2.

Figure 9 .

 9 Figure 9. (a) Forecasted number of events with m ≥ 2 per cell for the day of October 23 rd 2004 (logarithmic scale). Black circles represent observed earthquakes with m ≥ 2 which occurred during this day. Two of these events are aftershocks of the m = 6.5 2003/12/22 San-Simeon earthquake (located at lat=35.7 • and lon=-121.1 • ), three are associated with the m = 6 2004/09/28 Parkfield mainshock (lat=35.81 • and lon=-120.37 • ), and one is an aftershock of a m = 3.7 2004/09/09 earthquake (lat=35.09 • and lon=-117.52 • ). The predicted number of events for this day was 8.39 and the observed number was 6. Most of these earthquakes are better predicted by the time-dependent ETES forecast than by the long-term model (i.e., N p ( r) > µ LT ( r) in the cells within which these earthquakes occurred). Only one event, which occurred at 11 km away from the San-Simeon earthquake (lat=35.81 • and lon=-121.02 • ), was better predicted by the long-term model, because it occurred just outside of the main aftershock zone. (b) Ratio of the forecasted number of events estimated using the time dependent and long-term model (logarithmic scale). High values of N p /µ LT (up to 800) are associated with recent large earthquakes, such as Parkfield, San-Simeon, Landers, Hector-Mine and Northridge.

Figure 10 .

 10 Figure10. Observed (crosses), and predicted number of m ≥ 2 earthquakes per day as a function of the time since Landers mainshock, for model 2 (circles), model 3 (pluses), and model 6 (diamonds). The saturation at t ≤ 10 days is due to the incompleteness of the catalog for small magnitudes.

Figure 11 .

 11 Figure 11. (a) Observed (black) and predicted (gray) number of m ≥ 2 earthquakes per day in southern California for model 3 (seeTable 1). Dashed line is the background rate µ s = 2.81/day. (b) Probability gain per earthquake defined in (23).

  Figure 12. Predicted number of m ≥ 2 earthquakes per day for model 3 (see Table 1) as a function of the observed number, for the time period 1985-2003. The dashed line represents the perfect fit. The horizontal line is the background rate µ s = 2.81/day.

Table 1 .

 1 Model parameters, log likelihood LL, and probability gain G. The spatial kernel f (r) is either a Gaussian ('gs', see Eq. (15)) or a power-law ('pl', see Eq. (14)).

	model f (r) m min	α	p	K	µ s	f d	LL	G
	1	gs	2.	0.44 1.18 1.18 2.42	1.03	-430928 11.66
	2	pl	2.	0.43 1.19 1.29 2.09	0.06	-430754 11.69
	3	gs	2.	0.80 a 1.18 0.45 2.81	0.41 -434301. 11.07
	4	pl	2.	0.80 a 1.18 0.45 2.83 0.006 -435180. 10.93
	5 b	gs	2.	0.80 a 1.20 0.38 2.98	0.44 -436103. 10.77
	6 c	gs	2.	0.80 a 1.18 a 0.45 a 2.81 a 0.41 a -436684. 10.68
	7	gs	3.	0.80 a 1.18 a 0.45 a 2.81 a 0.41 a -55077. 11.94
	8	gs	4.	0.80 a 1.18 a 0.45 a 2.81 a 0.41 a	-7463.	8.46
	9	gs	5.	0.80 a 1.18 a 0.45 a 2.81 a 0.41 a	-986.	4.41
	10	gs	6.	0.80 a 1.18 a 0.45 a 2.81 a 0.41 a	-149.6	0.75
	11 d	gs	6.	0.80 a 1.18 a 0.45 a 2.81 a 0.41 a	-135.9	4.16

a Fixed parameter b Using the corrective term ρ * (m) defined by (

13

).

Table 2 .

 2 Comparison of the predicted number of m ≥ 6 events per day, for the days when a m ≥ 6 earthquake occurred, at the location of the earthquake (i.e., within the cell of 0.05 o × 0.05 o ), for the ETES model (N ET ES , using models 10 and 11 in Table1) and for the long-term model (N LT ).

	earthquake	date	time m	N ET ES	N LT	N ET ES /N LT
	Elmore-Ranch a	1987/11/24 01:54 6.2 4.3 × 10 -8 1.7 × 10 -7	0.25
	Elmore-Ranch b			4.5 × 10 -8		0.26
	Superstition-Hill a 1987/11/24 13:15 6.6 1.0 × 10 -7 4.0 × 10 -7	0.25
	Superstition-Hill b			7.8 × 10 -4		1950.
	Joshua-Tree a	1992/04/23 04:50 6.1 3.0 × 10 -7 9.6 × 10 -7	0.31
	Joshua-Tree b			2.9 × 10 -7		0.30
	Landers a	1992/06/28 11:57 7.3 1.2 × 10 -6 1.6 × 10 -6	0.75
	Landers b			1.4 × 10 -6		0.87
	Big-Bear a	1992/06/28 15:05 6.5 1.7 × 10 -7 6.2 × 10 -7	0.27
	Big-Bear b			5.5 × 10 -3		8871.
	Northridge a	1994/01/17 12:30 6.6 7.8 × 10 -8 2.8 × 10 -7	0.28
	Northridge b			7.9 × 10 -8		0.28
	Hector-Mine a	1999/10/16 09:46 7.1 6.5 × 10 -5 5.4 × 10 -7	120.
	Hector-Mine b			3.1 × 10 -7		0.57
	San-Simeon a	2003/12/22 19:15 6.5 4.2 × 10 -8 1.6 × 10 -7	0.26
	San-Simeon b			4.2 × 10 -8		0.26

a Forecasts updated each day at 00:00 universal time (model 10 in Table

1

) b Forecasts updated each day at 13:00 universal time (model 11 in Table

1)

Acknowledgments. We acknowledge the Advanced National Seismic System for the earthquake catalog. This work is partially supported by NSF-EAR02-30429, by the Southern California Earthquake Center (SCEC) and by the James S. McDonnell Foundation 21st century scientist award/studying complex system. SCEC is funded by NSF Cooperative Agreement EAR-0106924 and USGS Cooperative Agreement 02HQAG0008. The SCEC contribution number for this paper is 895.


Main Objectives of this talk

▶ Perform a (non-exhaustive) survey of proposals and their security status.

▶ Provide the intuition of the "most useful attacks" against LW ciphers.

▶ Conclusions and remarks (link with hash functions).

5/60

Survey of Proposals 1 ▶ Feistel Networks -best external analysis DESLX -none ITUbee -self-similarity (8/20r) LBlock -imposs. diff. (24/32r) SEA -none SIMON and SPECK -imposs. diff., diff, 0-correl. XTEA -mitm (23/64r) CLEFIA -imposs. diff. (13/18r) HIGHT -0-correlation (27/32r) TWINE -mitm,imposs. diff.,0-corr (25/36r) 1 mainly from https://cryptolux.org/index.php/Lightweight Block Ciphers 6/60 Survey of Proposals ▶ Substitution-Permutation Network KLEIN -dedicated attack (full round) LED -EM generic attacks (8/12r, 128K) Zorro -diff. (full round) mCrypton -mitm (9/12r, 128K) PRESENT -mult. dim. lin. (27/31r) PRINTcipher -invariant-wk (full round) PRIDE -diff (18/20r) PRINCE -mult. diff (10/12r) Fantomas/Robin -none/invariant-wk (full round) 7/60 Survey of Proposals ▶ FSR-based KTANTAN/KATAN -mitm (153

More Proposals

For more details, primitives, classifications, see: Given an input difference between two plaintexts, some output differences occur more often than others.

State of the Art in

X ′ X ′′ Y ′ Y ′′ ∆X ∆Y E K E K A differential is a pair (∆ X , ∆ Y ).
12/60 13/60

Truncated Differential Attacks [K 94]

A truncated path predicts only parts of the differences.

Let's see a simple example:

14/60

Truncated path: example

X X X X X X X X 0 0 0 X X X ? ? 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 15/60 Truncated path: example X X X X X X X X 0 0 0 X X X 0 0 0 0 0 0 0 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 16/60
F k 5 F k 6 F k 7 F k 8 <<< 8 <<< 8 <<< 8 F k 9 <<< 8 <<< 8 F k 10 F k 11 F k 12 <<< 8 <<< 8 F k 13 <<< 8 <<< 8 F k F k F k F k <<< 8 <<< 8 <<< 8 F k <<< 8 Impossible Differential Attack ∆X ∆Y ∆ in ∆ out r in r out r ∆ (c in , k in ) (c out , k out )
For the Attacks to Work

We need, for a state size s and a key size |K|:

C data < 2 s and C data + 2 |k in ∪k out | C N + 2 |K|-|k in ∪k out | P 2 |k in ∪k out | < 2 |K|
where C data is the data needed for obtaining N pairs (∆ in , ∆ out ),

C N is the average cost of testing the pairs per candidate key (early abort technique [LKKD08]) and P is the probability of not discarding a candidate key.

23/60

First Rounds

<<< 8 <<< 8 <<< 8 <<< 8 K 1 K 2 K 3 K 4 L 4 L 5 L 3 L 2 L 1 R 1 R 2 R 3 R 4 R 5 1 cond. 1 cond. 2 cond.
3 cond.

24/60

Last Rounds

<<< 8 <<< 8 <<< 8 <<< 8 K 19 K 20 L 19 R 19 3 cond. 2 cond. 1 cond. 1 cond. L 20 R 20 K 21 K 22 L 21 R 21 L 22 R 22 L 23 R 23 25/60
❋ ✵ ✐ ✵ ✐ ✶ ✐ ✷ ✐ ✸ ❋ ✶ ✐ ✶ ✵ ✐ ✶ ✶ ✐ ✶ ✷ ✐ ✶ ✸ |❑ ✷✐ ✷ |❑ ✷✐ ✶ 28/60

Multiple Impossible Differentials

Formalize the idea of [Tsunoo et al. 08]: CLEFIA has two 9-round impossible differentials ((0, 0, 0, A) ̸ → (0, 0, 0, B)) and ((0, A, 0, 0) ̸ → (0, B, 0, 0)) when A and B verify:

A B (0, 0, 0, α) (0, 0, β, 0) or (0, β, 0, 0) or (β, 0, 0, 0) (0, 0, α, 0) (0, 0, 0, β) or (0, β, 0, 0) or (β, 0, 0, 0) (0, α, 0, 0) (0, 0, 0, β) or (0, 0, β, 0) or (β, 0, 0, 0) (α, 0, 0, 0) (0, 0, 0, β) or (0, 0, β, 0) or (0, β, 0, 0)

B ′ = ■ ⊕ S 0 (■ ⊕ ■) (with B = B ′ ⊕ ■) |k in ∪ k out | = 122 bits ⇒ |k in ∪ k out | = 122-16 + 8 B ′ bits 30/60

Applications of Improved Impossible Diff

▶ CLEFIA: best attack on CLEFIA (13 rounds). ▶ Camellia: Improved best attacks for Camellia. ▶ AES: attacks comparable with best mitm ones (7 rounds). ▶ LBlock: best attack (on 24 rounds).

31/60

Meet-in-the-middle attacks

K 1 K 2 K=K 1 U K 2 2 | K 1 | + 2 | K 2 | + 2 | K| -s '

34/60

With Bicliques [KRS'11] ▶ Very large number of analysis published (20+).

X Ci p h e r t e x t B F E K P l a i n t e x t K=K 1 U K 2 2 | k 1 | + 2 | k 2 | + 2 | K 1 | + 2 | K 2 | + 2 | K| -s ' K 2 k 2 k 1 s ' K 1 35/
▶ Best attacks so far: multiple linear attacks (27r/31r).

46/60

PRESENT

Block n = 64 bits, key 80 or 128 bits.

⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ S15 S14 S13 S12 S11 S10 S 9 S 8 S 7 S 6 S 5 S 4 S 3 S 2 S 1 S 0 ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕
31 rounds + 1 key addition.

47/60

Forward Computation

Backward Computation

Sieving through the Sboxes: 1 Sbox

x 3 x 2 x x 0 S(x) 3 S(x) 2 S(
x 2 x 1 x 0 → S y 1 y 0 000→00 000→11 001→01 001→10 010→10 010→11 011→00 011→11 100→00 100→01 101→00 101→11 110→01 110→10 111→01 111→10 16
values of x 2 , x 1 , x 0 , y 1 , y 0 , out of 32, correspond to a valid transition.

Sieving through the Sboxes

❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙ ❙
▶ Probability for 1 Sbox p = 16/32 = 1/2 ▶ Probability for the 6 Sboxes: 1 ▶ Very large number of analysis published (20+).

Importance of Dedicated Cryptanalysis

▶ Best attacks so far: multiple linear attacks (27r/31r).

49/60

PRESENT

Block n = 64 bits, key 80 or 128 bits.

⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ S15 S14 S13 S12 S11 S10 S 9 S 8 S 7 S 6 S 5 S 4 S 3 S 2 S 1 S 0 ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ ⊕

▶

  Lightweight Symmetric Cryptography, by Alex Biryukov and Leo Perrin https://eprint.iacr.org/2017Impossible differentials (Feistel) ▶ Mitm / guess and determine (SPN, FSR) ▶ Dedicated: (differential/linear...)

▶

  Introduced by Diffie and Hellman in 1977.▶ Largely applied tool.▶ Few data needed.▶ Many improvements: partial matching, bicliques, sievein-the-middle...

▶▶

  Few cases broken by well known attacks (ex. Puffin or Puffin2 -multiple differentials) ▶ Happily, this is rare. Most of the times, new families or new ideas on known attacks exploiting the new properties are needed. ▶ Lightweight: more 'risky' design, lower security margin, simpler components. ▶ Often innovative constructions: One of the most popular ciphers, proposed in 2007, and now ISO/IEC standard.

Impossible Differential on LBlock

▶ For 21 rounds a complexity of 2

69.5

in time with 2 63 data, for 22: 2 71.53 time and 2 60 data, for 23: 2 75.36 time and 2 59 data. ▶ Feistel constructions in general are good targets 26/60 Improvements [BN-PS14,BLN-PS17,B18] ▶ Multiple impossible differentials (related to [JN-PP13]) ▶ Correctly choosing ∆ in and ∆ out (related to [MRST09]) ▶ State-test technique (related to [MRST09]) ▶ More accurate estimate of the pairs [B18] 27/60 Example: CLEFIA-128 • block size: 4 × 32 = 128 bits • key size: 128 bits • # of rounds: 18

  One of the most popular ciphers, proposed in 2007, and now ISO/IEC standard.

60 Bicliques ▶ Improvement of MITM attacks, but also... ▶ It can always be applied to reduce the total number of computations (at least the precomputed part) ⇒ acceleration of exhaustive search [BKR'11] 2 ▶ Many other accelerated exhaustive search on LW block ciphers: PRESENT, LED, KLEIN, HIGHT, Piccolo, TWINE, LBlock ... (less than 2 bits of gain). ▶ Is everything broken? No. ▶

  Linear cyptanalysis: because of the Sbox, a linear approximation 1 to 1 with bias 2 -3 per round [O-09]. Early abort ← message modification techniques ▶ State-test tech. & choosing ∆ in,out ← Rebound attacks ▶ Mult. impos. diff. ← mult. limited birthday distinguishers ▶ Using freedom degrees for bounds?... be careful!! ▶ Merging lists from rebounds/sieve in the middle → many applications ▶ Other ex: AES distinguishers inspired on rebound attacks. Classical attacks, but also new dedicated ones exploiting the originality of the designs. ▶ Importance on generalizing: improvements, and dedicated might become well stablished techniques. ▶ Importance of reduced-round analysis to re-think security margin, or as first steps of further analysis. ▶ New ideas inspired by SHA-3: might help improving attacks further! ▶ Better identifying composite problems/ list merging situations might provide improved results.

	PRESENT PRINTcipher PRINTcipher The Invariant Subspace Attack [LAAZ'11] The Invariant Subspace Attack Zorro -Hash Functions links LED -Hash Functions links Hash functions links -Sum up To Sum Up To Sum Up 3
	31 rounds + 1 key addition. ▶ Multiple linear attacks: consider several possible 50/60 approxs simultaneously ⇒ up to 27 rounds out of 31 [BN-14]. 51/60 With probability 1: ▶ Mitm, bicliques/initial structures: ▶ Lightweight block cipher proposed [GGN-PS13] for easy masking. used for both scenarios ▶ Many PRESENT-like ciphers proposed, like Puffin, PRINTcipher ▶ PRINTcipher[KLPR'10]: first cryptanalysis: invariant subspace attack[LAAZ'11]. 52/60 48 rounds. 53/60 ▶ Weak key attack, but a very bad property for 2 51 keys... 54/60 ▶ More applications afterwards: iScream, Robin, Zorro, Midori. ▶ Importance of generalizing/understanding dedicated attacks: new families/techniques might appear. 55/60 Final remarks ▶ A modified AES with only four sboxes per round (SPN with partial non-linear layer). ▶ Bounds on number of active Sboxes? Computed using freedom degrees. ▶ Many analyses published. Problem: MC property ⇒ devastating attack [BDDLT13, RASA13] 56/60 ▶ Analysis provided with respect to known key distinguishers (rebound-like). Seems like a lot of SHA-3 knowledge put into this design. 57/60 Conclusion 59/60 A lot of ciphers to analyze/ a lot ▶ 58/60 of work to do!

▶ Usually, weaker than the original.

▶ Lightweight block cipher proposed in [GPPR12].

▶ AES-like with simpler key-schedule and more rounds.

Nice simple design.

▶

in total: C data = 2 113 becomes C data = 2 113 /24

B = ■ ⊕ S 0 (■ ⊕ ■) ⊕ ■

Most important application: best key-recovery on AES-128 in 2 126.1 instead of the naive 2 128 .

6▶ We only try 2 80-6 = 2 74 potential key candidates. ▶ 7 rounds (+1 bicliques).

Thank you to Christina Boura and Leo Perrin for their help with the figures and the slides. 60/60

/254r) Grain -correl./ cube attacks (some full) Trivium -cube attacks (800/1152) -Sprout -guess-and-determine (full round) Quark -condit. diff (25%) Fruit -divide and conquer (full) Lizard -guess-and-det. (full) 8/60 Survey of Proposals ▶ ARX Chaskey -diff-lin (7/8r) Hight -0-correl (27/32r) LEA -diff. (14/24r) RC5 -diff. (full round) Salsa20 -diff (8/20r) Sparx -imposs. diff. (15/24r) Speck -diff. (17/32r) 9/60

Discarding Wrong Keys

▶ Given one pair of inputs with ∆ in that produces ∆ out , ▶ all the (partial) keys that produce ∆X from ∆ in and ∆Y from ∆ out differ from the correct one.

▶ If we consider N pairs verifying (∆ in , ∆ out ) the probability of NOT discarding a candidat key is

Improved Bicliques [CN-PV 13]

Can we build bicliques with only one pair of P-C?

38/60

Sieve-in-the-Middle [CN-PV'13] ▶ Compute partial inputs and outputs of S ⇒ sieving with transitions instead of collisions.

When can we sieve?

▶ n in known bits out of m: at most 2 m-n in values for the n out output bits. ▶ A transition exists with probability p. ▶ Sieve when n in + n out > m ⇒ p < 1

40/60

How do we sieve?

▶ We obtain a list L A of partial inputs u and a list L B of partial outputs v ⇒ merge L A and L B with the condition (u, v) is a valid transition though S.

▶ Naive way costs |L

no gain with respect to exhaustive search.

▶ We need an efficient procedure.

Often S is a concatenation of S-boxes.

41/60

Merging the lists

Merging the lists with respect to R ▶ R is group-wise, i.e. for z groups

First studied for rebound attacks.

Group-wise relation

Merging Algorithms

▶ Problem also appears in divide-and-conquer attacks (and rebound attacks).