Gwladys Ravon

Yves Coudière

Angelo Iollo

Oliver Bernus

Richard D Walton

Issues in Modeling Cardiac Optical Mapping Measurements

Keywords: mathematical modeling, optical mapping, inverse problem

Optical mapping allows to visualize cardiac action potentials (AP) on cardiac tissue surfaces by uorescence using voltage-sensitive dyes. So far, the surface measurements are directly related to surface AP. In a previous study was developed a method to reconstruct threedimensional depolarization front: the main idea was to solve an inverse problem using the experimental measures on the surfaces. Although the method was very accurate on in silico data, it showed diculties to recover real optical mapping measurements. Here we describe the dierent directions we followed to improve the results.

Introduction

Optical mapping is an important tool for the understanding of cardiac arrhythmias [START_REF] Rosenbaum | Optical mapping of cardiac excitation and arrhythmias[END_REF]. It provides surface optical recordings that are linked to surface AP [START_REF] Bishop | Synthesis of voltage-sensitive optical signals: application to panoramic optical mapping[END_REF]. Although the photons are known to interact with the tissue up to a few mm in depth, it remains very challenging to actually retrieve 3D information from optical recordings. In [START_REF] Khait | Method for 3-dimensional localization of intramyocardial excitation centers using optical imaging[END_REF], Khait et al derive a formula to determine the depth of some xed electrical sources in a phantom. Instead, we propose to solve an inverse problem, so as to recover more accurate and complete information on the electrical sources. The approach is also expected to apply to more general experimental conditions. In a previous study we presented a rst attempt to obtain a complete 3D reconstruction following this approach. Although we obtained excellent agreement between the actual location of the source and the location found by the inverse method, on in silico data, the results on the experimental data provided by the authors of [START_REF] Khait | Method for 3-dimensional localization of intramyocardial excitation centers using optical imaging[END_REF] were disappointing. Actually, we observed a large mismatch between the experimental measures and the measures computed from the optical mapping model and the known experimental locations. Here we explore the model and its optical parameters as a cause of this mismatch.

Optical mapping

For optical mapping of AP, uorescent voltage-sensitive dyes that attach to the cells' membrane are injected into a slab of tissue. The tissue is put in a bath, and cameras and lights are placed on both sides of the preparation (epicardium and endocardium, g. 1). A lter is placed in front of each camera, that allows to choose the wavelength to be recorded. The epi-and endocardium are alternatively illuminated. The dyes then emit a uorescent light assumed to be proportional to both the incident light and the transmembrane potential (TMP). Optical images are the surface recordings of this light by the cameras [START_REF] Walton | Experimental validation of alternating transillumination for imaging intramural wave propagation[END_REF]. Images are recorded on the illuminated surface (reexion mode) and on the opposite one (transillumination mode). For each time step, optical mapping hence produces four images.

The medium has a natural uorescence (called background F 0 ) recorded when the tissue is at rest. Fluxes captured during an AP are denoted by F . The signal due only to the AP itself is F -F 0 . We shall rather use the usual renormalization:

g = F 0 max 0, F -F 0 F 0 .
(

) 1
Indeed the max(. . .) amounts to ignore negative, physically irrelevant, optical signals (due to noise). The multiplication by the average F 0 of the background signal is a way to retrieve the correct amplitude of the signal. Our main goal is to reconstruct the 3D front of the AP from these 2D optical data. 3 Model

Forward problem

In order to write the mathematical model of these observations, we assume the following: the cameras record photon uxes through the surfaces, the light interacts with the tissue material in the diusive regime [START_REF] Khait | Method for 3-dimensional localization of intramyocardial excitation centers using optical imaging[END_REF], and a Robin boundary condition can be used to model the interaction between the tissue and its environment. Hence the illumination light is described by its photon density φ 0 that solves the diusion equation

   -D 0 ∆φ 0 + µ 0 φ 0 = 0 in Ω, φ 0 + d 0 ∂φ 0 ∂n = 0 on ∂Ω \ Γ, and φ 0 = I 0 δ 0 D 0 on Γ, (2)
where Ω ⊂ R 3 represents the slab of tissue, Γ is the illuminated surface, and n is the unit normal to ∂Ω, outward of Ω. The uorescent light is assumed to be proportional to the TMP and the illumination light (multiplicative factor β > 0). Its photon density solves:

   -D∆φ + µφ = β(V m -V 0 )φ 0 in Ω, φ + d ∂φ ∂n = 0 on ∂Ω. (3)
In both equations, the optical parameters D, µ, d stand respectively for diusion coecient, absorption coecient and extrapolation distance. The attenuation length is the parameter δ 0 dened by δ 0 = D0 µ0 . The intensity of the illumination, assumed uniform, is the parameter I 0 . The multiplicative factor β > 0 is known for the dye used during the experiments. The dyes are assumed to be uniformly distributed in the tissue. Finally the uxes measured through the surfaces are given by Fick's law:

g = -D ∂φ ∂n
on the epi or endocardium.

Remark that the experimental ux g given by ( 1) does not satisfy equation ( 3), because of the renormalization. The quantity F -F 0 does. However we shall consider g as a good approximation of g, following the recommendations of the experimenters.

Since we consider a rectangular slab of tissue, we may have used structured meshes. We choose to work with unstructured meshes in order to allow more general geometries. This is necessary to study data from heart tissues. The diffusion equations are solved with with P1-Lagrange nite elements method using the solver FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF].

Inverse problem

The problem of retrieving the 3D spatial distribution of the TMP, denoted by V m (t, x) from the 2D optical signals at time t > 0 is under-determined. Hereafter, x = (x, y, z) denotes a point in Ω with Cartesian coordinates (x, y, z). Instead of nding the complete distribution V m (t, x), we look for a depolarization front at each time. Specically, we assume that a surface S(t) = {x ∈ Ω : f (t, x) = 0} dened as the level 0 of the function f splits the domain Ω into the region Ω r = {x ∈ Ω : f (t, x) > 0} of tissue at rest, and the region

Ω p = {x ∈ Ω : f (t, x) < 0} of excited tissue. It follows that V m (t, x) = V p if x ∈ Ω p , and V m (t, x) = V 0 if x ∈ Ω r .
We consider simple depolarization fronts S modeled either by the sphere centered in x 0 ∈ Ω and expanding with the velocity c > 0 after the given time t = t 0 ≥ 0, dened by the level-set function

f (t, x) = |x -x 0 | -c(t -t 0 ),
or by the xed ellipsoid centered in x 0 ∈ Ω and with radiuses r x , r y , r z > 0, dened by the level-set function

f (t, x) = (x-x0) 2 r 2 x + (y-y0) 2 r 2 y + (z-z0) 2 r 2 z
This level-set approach generalizes to more complex AP, once these simple cases are completely understood. In both cases, the inverse problem reduces to the identication of a small parameters set P = (x, c, t 0 ) ⊂ R 5 (sphere) or P = (x, r x , r y , r z ) ⊂ R 6 (ellipsoid). In order to identify these parameters, we minimize the least squares dierence e(P) between the actual measurements and the measurements computed from equations ( 2), (3), and (4) with a TMP as above:

e(P) = 1 2 4 i=1 g i P -g ,i 2 L 2 (Si) , (5)
where the functions g ,i are the data. Here i refers to one of the four images (i ∈ {1, 2, 3, 4}), and the surface S i is either the epicardium or the endocardium, as detailed in Table 1. Although this is the natural way to dene the cost function,

# surface measured surface 1 epicardium S1 = epicardium 2 epicardium S2 = endocardium 3 endocardium S3 = endocardium 4 endocardium S4 = epicardium Table 1.
References of the measures the value I 0 of the illumination in equation ( 2) is unknown, while the optical parameters are. Consequently, and since equations ( 2), (3), and (4) are linear, the density φ 0 , or φ, can only be computed up to a multiplicative constant. The mapping I 0 → g i is also linear, the measurement in-silico g i is consequently proportional to I 0 , and we can change the cost function to account for this unknown value. A rst idea is to identify the intensity I 0 , and consider the following modied cost function:

e(P) = 1 2 2 i=1 2 j=1 I i 0 g ij P -g ,ij 2 L 2 (Sj ) , (6)
where i, j ∈ {1, 2}, i stands for the illuminated surface while j stands for the measured surface (i = j for the reexion mode, and i = j for the transillumination mode). Since the intensities are dierent, we have two additional parameters to retrieve, I 1 0 and I 2 0 . A second idea consists in comparing the normalized uxes:

e(P) = 1 2 4 i=1 g i g i - g ,i g ,i 2
L 2 (Si) , (7)
where the inner norms • are also L 2 norms on the surface S i . In this case there is no additional parameter to be identied, but the problem becomes nonlinear with respect to V m . In all cases, a xed-step gradient method followed by the BFGS algorithm is used to solve the inverse problem. We need the gradient of all the cost functions e with respect to the unknown parameters P, which is computed with the adjoint method. [START_REF] Hecht | New development in FreeFem++[END_REF] Summary of the results and issues

Summary of previous results

In this part we quickly recall the rst results we obtained. We rst compared the formula derived in [START_REF] Khait | Method for 3-dimensional localization of intramyocardial excitation centers using optical imaging[END_REF] and our method for the expanding sphere and the cost function [START_REF] Bishop | Synthesis of voltage-sensitive optical signals: application to panoramic optical mapping[END_REF] on data in silico.

Figure 2 shows the depth of the source computed from Khait's formula (diamonds) and from our method (squares) as a function of time. Results are presented for inclusions at four dierent locations, excitation time t 0 = 0 and velocity c = 0.5 m.s -1 . The last example was carried out on a cylinder, in order to illustrate the case of a complex geometry. For all cases we retrieved the complete Excit at ion at (10,10,2.54) (r ounded slab)

depth (m m ) Fig. 2. Comparison between with Khait's formula (diamonds) and our method (squares). The vertical lines mark the breakthrough of the wave on the observed surface.

unknown location x 0 of the source, at any time, and even after the breakthrough, with an accuracy up to machine precision. Additionally, the velocity c could be recovered from any time-sequence of data. These were very good results, based on data in-silico.

The method was then tested on the experimental data from the optical phantom experiments set up by Khait and coworkers [START_REF] Khait | Method for 3-dimensional localization of intramyocardial excitation centers using optical imaging[END_REF]. In this case, xed ellipsoidal sources were considered and the intensity I 0 of the illumination was unknown. We chose to work with the normalized cost function (7), so that we had to identify six parameters.

Figure 3 shows the results obtained for a phantom located at a depth z 0 = 13 mm. Although the reconstruction of the photon uxes looks qualitatively correct (rst and second rows of images), the reconstructed depth is z = 16.7 mm. In order to understand this large error (3 mm over a total depth of 20 mm), we computed the theoretical observations associated to the exact experimental location of phantom source: from P known, we compute the TMP distribution V m (x), solve equations ( 3) and ( 4), and nally compute the observation g with equation (4). These normalized uxes (third row of images on gure 3) clearly have dierent amplitudes than the experimental signal. These uxes are also more diuse than the experimental ones. We observe the same behavior with seven other experimental phantoms. We tried to replace the measure g (eq. ( 1)) by the simpler dierence max (F -F 0 , 0) and obtained similar errors on the location. We deduce from these results that the optical model ( 2) and ( 3) is questionable. In the next section, we identify and study in depth several possible sources of error.

Model improvements to x the mismatch

Uniform Illumination. First, we addressed the approximation of uniform distribution of the illumination I 0 . We tested the eect of several spatially distributed illumination intensities I 0 (x, y): constant, narrow Gaussian, diuse Gaussian, supposed to mimic the experimental lights. They read

I 0 (x, y) = A exp - (x -x 0 ) 2 σ 2 x + (y -y 0 ) 2 σ 2 y + (x -x 0 )(y -y 0 ) 2σ xy . (8)
We observed the same uxes g for the three cases. Indeed, the source term (V m (x) -V 0 )φ 0 (x) in the uorescence equation is the same in all cases because the photon density φ 0 becomes quickly constant inside the tissue (see Figure 4). Finally it is not a limitation to consider an uniform illumination. Optical parameters. Another possibility is that the values of the optical parameters are not correct. In order to understand the role of these parameters, we derive a dimensionless version of the equations. Consider a space scale L > 0 and let rescale x as : x = x L . Equations ( 2) and ( 3) rewritten

             - D 0 L 2 ∆φ 0 + µ 0 φ 0 = 0 in Ω, φ 0 + d 0 L ∂φ 0 ∂n = 0 on ∂Ω \ Γ, φ 0 = I 0 δ 0 D 0 on Γ,      - D L 2 ∆φ + µφ = β(V m -V 0 )φ 0 in Ω, φ + d L ∂φ ∂n = 0 on ∂Ω.
We dene γ = I0δ0 D0 and rescale the density as φ 0 = φ0 γ and φ = φ γ , in such a way that:

           - D 0 L 2 ∆φ 0 + µ 0 φ 0 = 0 in Ω, φ 0 + d 0 L ∂φ 0 ∂n = 0 on ∂Ω \ Γ, φ 0 = 1 on Γ,      - D L 2 ∆φ + µφ = β(V m -V 0 )φ 0 in Ω, φ + d L ∂φ ∂n = 0 on ∂Ω.
Dividing the diusion equation by µ 0 (resp. µ) the dimensionless system reads:

         -δ 0 2 ∆φ 0 + φ 0 = 0 in Ω, φ 0 + d 0 ∂φ 0 ∂n = 0 on ∂Ω \ Γ, φ 0 = 1 on Γ,    -δ 2 ∆φ + φ = (V m -V 0 )φ 0 in Ω, φ + d ∂φ ∂n = 0 on ∂Ω.
where δ 0

2 = D0 µ0L 2 , d 0 = d0 L , δ 2 = D µL 2 and d = d
L are dimensionless optical parameters. Finally the uxes are given by g = φ on the epi or endocardium. The other parameters, I 0 and β are hidden in a dimensionless number γ and we compare the experimental uxes g to γg . The optical system is characterized by the four optical parameters. Adding the two terms γ (one for each illumination), we end up with six parameters. We tried to solve a second inverse problem: knowing the characteristics P of the inclusion, identify the six new parameters by minimizing the cost function

e 1 (δ 0 2 , d 0 , δ 2 , d, γ 1 , γ 2 ) = 1 2 2 i=1 2 j=1 γ i g ij -g ,ij 2 L 2 (Sj ) . (9)
To date, the numerical solutions to this optimization problem are still being computed.

Other possibilities. In our model we do not consider the distance between the preparation and the camera. We could consider the diusion of the photon density in the air by ensuring the continuity of the uxes at the border medium/air. Instead, we impose a Robin condition on all the surfaces. We consider that the uxes are recorded directly through the tissue surfaces, and not through the Plexiglas, because the Plexiglas has a negligible absorption coecient.

Discussion

Modeling the optical measurements by the diusion equations ( 2) and ( 3) is widely used in cardiac optical mapping. But when we confront the measurements obtained with this model to the experimental ones in a well controlled setup, we observe an important mismatch (g. 3). We studied several ideas that might explain the dierences, always with a negative result. To our opinion, this is likely to suggest that the diusive regime is a too coarse approximation of the interaction between light and matter in the cardiac context. If this is conrmed by further experiments, the complete radiative transfer equation (RTE) might be used to model the measurements.

Conclusion

The aim of this study was to detail our investigations concerning the current model of cardiac optical mapping measurements. Having recalled previous results we described our attempts to improve our mathematical model. The most likely assumption was that the illumination on the tissue was not constant. Few tests showed that it was not the key. We eliminated several other sources of error, but there remains some more. We keep on working on the dimensionless problem in order to identify its parameters. Otherwise the RTE might be used.

Fig. 1 .

 1 Fig. 1. The optical imaging setup: (1) CCD camera, (2) emission lter, (3) LED illumination, (4) tissue sample, (5) ECG electrode, (6) bipolar stimulating electrode.

  at ion at (x,y,z)= (10,10,0.5)

Fig. 3 .

 3 Fig. 3. Results for one set of experimental data. Domain size: 40x40x20mm. Columns 1 and 2: records on the epicardium. Columns 3 and 4: records on the endocardium.

Fig. 4 .

 4 Fig. 4. Narrow Gaussian example for an epi-illumination. Domain size: 40x40x20mm. On the left: spatial distribution of I0, on the right: photon density φ0, cut plane in z-direction.

Acknowledgments

This work is partially supported by the grant number ANR-10-IAHU-04 from the French government, and by the Conseil Régional Aquitaine.


Marie-Madeleine Mialon

Karen Helle

Karen Helle Sloth

Isabelle Veissier
email: isabelle.veissier@clermont.inra.fr

Real time positioning to detect early

come

Introduction

Biological rhythms are essential for life. While diurnal animals are more active during the day, nocturnal ones are more active at night. This pattern is species specific. There is evidence that the circadian rhythm of activity is disrupted under chronic stress. For instance, in rodents (nocturnal species), more activity is observed during the dark period and less during daylight when animals are submitted to a repeated stressors or inescapable electrical shocks [START_REF] Stewart | Circadian rhythmicity and behavioral depression: I. Effects of stress[END_REF][START_REF] Solberg | Circadian rhythms and depression: effects of exercise in an animal model[END_REF]. We also noticed in cattle that the circadian variations of activity were less marked two days before the occurrence of symptoms of pneumonia [START_REF] Veissier | The use of circadian behaviour to measure adaptation of calves to changes in their environment[END_REF]. These latter results were obtained by serendipity, during an experiment on the stabilisation of activity rhythms after the animals were turned from pasture to indoor conditions and during which an outbreak of pneumonia occurred in the course of the observations. These results had never been confirmed, since they require heavy observations. We investigated if we could detect a circadian rhythm of activity of dairy cows in commercial settings using a real-time positioning system and if this rhythm depends on the state of the animal.

Material and methods

During the EU-PLF project (http://www.eu-plf.eu/), INRA and GEA Farm Technologies collaborate to run field tests using the GEA CowView system (http://www.geafarmtechnologies.com/hq/en/mediacenter/news/2012/dlgcowview.aspx). The CowView is based on real time positioning of animals, with a precision of 50 cm. The position of the cow is allocated to the different zone structures in the barn (cubicles, feeders, corridors…). The activity of a cow is inferred from its allocation and its movements to:

-the cow is considered to be resting if it is found in a cubicle, -the cow is supposed to be feeding if it is classified within the feed bunk zone, -the cow is supposed to be standing in the alley way if it does not move in the alley zone, -the cow is assumed to be walking in the alley way if it moves in the alley zone.

The data from the 350 cows from a commercial farm based in Denmark were collected for five months. They consisted of hourly accumulated activities as determined by the CowView system and all events recorded by the farmer (interventions on the herd, oestrus, lameness, mastitis, accident, calving, respiratory problem, diarrhoea,…).

To determine activity levels and their variations during the day, we needed to attribute a weight to each of the four activities distinguished. Those weights were calculated according to the method proposed by Veissier et al. [START_REF] Veissier | Calves' responses to repeated social regrouping and relocation[END_REF]. The weights were obtained from a factorial analysis of correspondences. On the first axis of the analysis, the four activities obtained the following weights: resting, -.19; standing, +0.12; walking, +.26; and feeding, +.29. This axis was considered to express the level of activity. For each cow, we then calculated its level of activity per day and per hour within each day and the size of circadian variations (the sum of the squares of the deviation of each hour from the average of the day).

Finally, we ran a variance analysis to relate the average activity level of a cow on a given day and its circadian variations to the events recorded by the farmer (Proc Mixed procedure of SAS for repeated measurements).

Results

The cows were active from 08:00 to 21:00, with two peaks of activity, one in the morning and one in the afternoon. The activity was minimal at night (Figure 1).

We divided events recorded on the farm in five categories: control days, oestrus, lameness, infectious disease, and interventions that may affect cows' behaviour such as vaccination, claw trimming… (the data from the days with an intervention were discarded from the analyses). We observed a higher activity in case of oestrus and infectious diseases. The circadian variations were less marked in case of infectious diseases. They were slightly less marked in case of oestrus and lameness (Table 1).

Discussion

From the results obtained in this study it seems possible to distinguish the state of a cow: a diseased cow seems to be more active (at least at the very beginning of a disease) but its activity is less variable during the day. The same is observed for a cow in oestrus but to a lesser extent. Finally a lame cow seems to have a less variable activity during the day. If these results are confirmed in other settings, they offer the possibility to design models to predict physiological or pathological states of cows from real time positioning.

Figure 1 :

 1 Figure 1: Diurnal variations of activity-level in 350 cows in a Danish dairy herd observed for five months.

Table 1 :

 1 Average activity level and circadian variations in activity level in the 350 cows from a Danish dairy herd observed for five months.

		Control	Oestrus Lameness Infectious	SEM	F	P
		days			diseases
	Activity	-0.285 a*	0.90 b	-0.17 a	1.73 c	0.107	249	<0.0001
	level
	Circadian	2.83 c	2.67 b	2.63 b	2.32 a	0.05	54	<0.0001
	variations
	*on each row values with no common letter differ significantly

Acknowledgements

This study was part of the European project no: 311825 EU-PLF (Animal and farm-centric approach to precision livestock farming in Europe), co-financed by the European Commission.


O Lengliné
email: lengline@unistra.fr

J.-P P Ampuero
email: ampuero@gps.caltech.edu

Insights on Earthquake Triggering Processes from 1 Early Aftershocks of Repeating Micro-Earthquakes 2

Characterizing the evolution of seismicity rate of early aftershocks can yield important information about earthquake nucleation and triggering. However, this task is challenging because early aftershock seismic signals are obscured by those of the mainshock. Previous studies of early aftershocks employed high-pass filtering and template matching, but had limited performance and completeness at very short times. Here we take advantage of repeating events previously identified on the San Andreas fault at Parkfield and apply empirical Green's function deconvolution techniques. Both Landweber and sparse deconvolution methods reveal the occurrence of aftershocks as early as few tenths of a second after the mainshock. These events occur close to their mainshock, within one to two rupture lengths away. The aftershock rate derived from this enhanced catalog is consistent with Omori's law, with no flattening of the aftershock rate down to the shortest resolvable time scale ∼ 0.3 s. The early aftershock rate decay determined here matches seamlessly the decay at later times derived from the original earthquake catalog, yielding a continuous aftershock decay over time scales spanning nearly 8 orders of magnitude. Aftershocks of repeating micro-earthquakes may hence be governed by the same mechanisms from the earliest time resolved here, up to the end of the aftershock sequence. Our results suggest that these early aftershocks are triggered by relatively large stress perturbations, possibly induced by aseismic afterslip with very short characteristic time. Consistent with previous observations on bimaterial faults, the relative location of early aftershocks shows asymmetry along-strike, persistent over long periods.

Introduction

Earthquake aftershock sequences are one of the most abundant manifestations of seismic activity and earthquake interactions. A robust characteristic is that their seismicity rate, λ(t), decays as a function of time t after the mainshock as a power-law well described by the Omori-Utsu law [START_REF] Utsu | The centenary of the Omori formula for a decay law of aftershock activity[END_REF],

λ(t) = K/(t + c) p .
(

) 1
where the exponent p is usually ∼ 1, the aftershock productivity K generally depends on mainshock magnitude, and the time scale c marks the onset of the power-law regime.

Despite the robustness of this empirical observation, the detailed mechanism responsible for aftershock sequences is still elusive and represents a major challenge for understanding the physics of earthquake nucleation and triggering. Several seismicity models, invoking very different physical mechanisms, have been equally successful at explaining the powerlaw decay of aftershock rates. These include rate-and-state nucleation models driven by static stress transfer from the mainshock [START_REF] Dieterich | A constitutive law for rate of earthquake production and its application to earthquake clustering[END_REF], Coulomb earthquake triggering models driven by postseismic slip [e.g. [START_REF] Schaff | Postseismic response of repeating aftershocks[END_REF] or by pore fluid diffusion and poroelastic stress transfer [e.g. [START_REF] Bosl | Aftershocks and pore fluid diffusion following the 1992 Landers earthquake[END_REF], and mechanisms triggered by dynamic stresses carried by the mainshock wavefield such as modifications of the fault zone properties [e.g. [START_REF] Parsons | A hypothesis for delayed dynamic earthquake triggering[END_REF] and aseismic slip transients [e.g. [START_REF] Shelly | Triggered creep as a possible mechanism for delayed dynamic triggering of tremor and earthquakes[END_REF]. Hence the power-law regime of aftershock decay rates contains limited information to discriminate among earthquake triggering models. However, at early times after the mainshock the predictions of aftershock models diverge [e.g. [START_REF] Helmstetter | Afterslip and aftershocks in the rate-and-state friction law[END_REF], suggesting that quantification of early aftershock rates can provide constraints on the physical D R A F T August 31, 2015, 4:43am D R A F T in Eq. 1 has a physical origin [e.g. [START_REF] Dieterich | A constitutive law for rate of earthquake production and its application to earthquake clustering[END_REF][START_REF] Ampuero ; Aftershocks X -27 Narteau | Common dependence on stress for the two fundamental laws of statistical seismology[END_REF] or mainly results from observational bias [START_REF] Kagan | Short-term properties of earthquake catalogs and models of earthquake source[END_REF][START_REF] Kagan | Relation between mainshock rupture process and omori's law for aftershock moment release rate[END_REF][START_REF] Helmstetter | Comparison of short-term and time-independent earthquake forecast models for southern california[END_REF]. Indeed, incompleteness of earthquake catalogs at early times produces an apparent saturation of the earthquake rate right after the mainshock.

Previous efforts to detect early aftershocks involved either amplitude-threshold detectors on high-pass filtered waveform envelopes [e.g. [START_REF] Peng | Anomalous early aftershock decay rate of the 2004 Mw6. 0 Parkfield, California, earthquake[END_REF][START_REF] Peng | Seismicity rate immediately before and after main shock rupture from high-frequency waveforms in Japan[END_REF][START_REF] Enescu | Quantifying early aftershock activity of the 2004 mid-Niigata Prefecture earthquake (Mw6. 6)[END_REF][START_REF] Enescu | Omori-Utsu law c-values associated with recent moderate earthquakes in Japan[END_REF] or matched filter detectors [e.g. [START_REF] Peng | Migration of early aftershocks following the 2004 Parkfield earthquake[END_REF][START_REF] Lengliné | Decay and expansion of the early aftershock activity following the 2011, Mw9. 0 Tohoku earthquake[END_REF].

Most of these studies recover a large number of early aftershocks and show an apparent decrease of aftershock rate at early times. However, the completeness achieved by these methods soon after the mainshock is much poorer than at later times, hampering the quantitative analysis of aftershock decay rates over a broad time scale range. [START_REF] Enescu | Omori-Utsu law c-values associated with recent moderate earthquakes in Japan[END_REF] inferred, from the study of early aftershocks of relatively large earthquakes, c-values on the order of their lowest resolvable time scale ∼ 1 min. Hence no flattening of the Omori decay was observed and the c-value could be even smaller. [START_REF] Sawazaki | Imaging the high-frequency energy radiation process of a main shock and its early aftershock sequence: The case of the 2008 iwate-miyagi nairiku earthquake, japan[END_REF] found a transition of aftershock rate behavior between 10 and 40 s after a Mw6.9 earthquake and attributed it to the transition between the dynamic and static triggering regimes.

Here we propose a novel strategy to quantify early aftershock activity. We consider a composite aftershock sequence obtained by stacking multiple sequences relative to mainshock time and location. We focus on earthquakes with very similar waveforms, which facilitates the uniform detection of aftershocks within their mainshock's coda. We present

D R A F T August 31, 2015, 4:43am D R A F T LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS X -5
two deconvolution approaches which capture aftershocks in the first 20 s following a mainshock. We show that we can detect events as early as 0.3 s after their mainshock, that is, a time scale of about ten times the mainshock rupture duration. Comparing aftershock rates in our newly identified events with those at later times, we find a constant power-law decay of the aftershock rate from the earliest resolvable time (0.3 s) up to about 100 days.

Our results suggest that the mechanism driving aftershock activity remains similar over a broad range of time scales spanning 8 orders of magnitude.

Data

We used the repeating earthquake catalog of [START_REF] Lengliné | Inferring the coseismic and postseismic stress changes caused by the 2004 M w = 6 Parkfield earthquake from variations of recurrence times of microearthquakes[END_REF], which comprises events with magnitudes ranging from M L = 0.9 to M L = 3.2 that occurred between 1984 and June 2007. This dataset contains 2414 events distributed in 334 repeating earthquake sequences (RES) formed by linking earthquake pairs with similar waveforms, overlapping source areas and similar magnitudes (see details in [START_REF] Lengliné | Inferring the coseismic and postseismic stress changes caused by the 2004 M w = 6 Parkfield earthquake from variations of recurrence times of microearthquakes[END_REF]). The dataset is provided in the supplementary material. The analysis is based on waveforms recorded by the short-period vertical sensors of the Northern California Seismic Network (NCSN) (Figure 1). We selected 20.48 s long signal segments, sampled at 100 Hz (2048 samples), starting 1 s before the P-wave arrival. We used all NCSN stations with at least one phase pick for the earthquakes considered here. We used all available records at these NSCN stations to form the dataset. For reference, 12 out of these 396 stations recorded half of the events; they are located at less than 20 km from the event epicenters.

D R A F T August 31, 2015, 4:43am D R A F T X -6
LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS

Methods

Deconvolution of Repeating Earthquake Signals

Identifying the seismic signature of very early aftershocks is not straightforward, because it overlaps with the coda waves of the mainshock (Figure 3). For instance, a template matching approach has limited detection capacity on such early signals (see Appendix A).

Here we propose a deconvolution method for reliable early aftershock detection. Previous related work includes [START_REF] Fischer | Modelling of multiple events using empirical greens functions: method, application to swarm earthquakes and implications for their rupture propagation[END_REF] and [START_REF] Wang | Compound earthquakes on a bimaterial interface and implications for rupture mechanics[END_REF]. Our analysis focuses on nearby events and exploits waveform similarity to facilitate the search. We note u k i (t) the waveform of the i th event of a given sequence recorded at station k. In the Fraunhofer approximation, this waveform can be represented as the temporal convolution ( * ) between the apparent source time function (ASTF), f k i (t), which depends on the source-receiver configuration, and a Green's function, G k i (t), that incorporates all effects related to seismic wave propagation and instrument response:

u k i (t) = [G k i * f k i ](t) . = ∫ +∞ -∞ G k i (t -t ′ )f k i (t ′ )dt ′ . ( 2
)
Our first objective is to recover f k i (t) by deconvolution. The source time functions will then be analyzed to identify if any aftershock occurs during the 20.48 s long time window studied. The deconvolution process is based on empirical Green's functions (EGF), i.e.

the waveform of an earthquake located close to earthquake i and recorded at the same station k, such that the nearly common ray paths for the two earthquakes implies a similar G k i . The ASTFs obtained by EGF deconvolution are not absolute, but relative to the ASTFs of the EGF event. We denote the resulting relative source time function as RSTF. In studies addressing in detail the rupture process of a target earthquake i, the EGF is chosen as a significantly smaller event to allow a point-source approximation. In

D R A F T August 31, 2015, 4:43am D R A F T LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS X -7
contrast, our goal here is to detect secondary events in the wake of earthquake i which is part of a RES, without attempting to resolve their individual rupture process. We hence choose as EGF event another earthquake of the RES. As all events in a RES have very similar magnitudes, their individual source time functions are very similar. Hence, if no aftershock occurs during the investigated time window the RSTF should be a single Dirac delta function. If an aftershock sufficiently similar (located closely) to the other events of the RES occurs, another peak should emerge in the RSTF (Figure 4).

To retrieve the RSTFs we apply the projected Landweber deconvolution algorithm [e.g. [START_REF] Piana | Projected Landweber method and preconditioning[END_REF][START_REF] Bertero | Application of the projected Landweber method to the estimation of the source time function in seismology[END_REF][START_REF] Vallée | Stabilizing the Empirical Green Function Analysis: Development of the Projected Landweber Method[END_REF]. This deconvolution procedure is an iterative process and we set f = 0 as a starting guess. The target waveform, u k i (t), and the EGF waveform, u k j (t), are aligned to the nearest sample based on the time shift of the maximum of their cross-correlation function. If the normalized cross-correlation of the entire 20.48 s of signals is higher than 0.7, we proceed with the deconvolution process.

This criterion is fulfilled by most of the earthquake pairs we considered, because events in a RES have very similar waveforms by definition. The purpose of this step is to remove stations with very noisy signals that produce poor reconstruction of the RSTF.

We assess through simulations whether this selection criterion excludes waveforms with shortly separated earthquake doublets that may have low correlation coefficient due to destructive interference. We build a synthetic doublet waveform as the sum of the original waveform and a shifted version of it. We then compute the correlation coefficient between the original and the synthetic doublet waveform. We test each possible time shift and repeat this operation for a random set of 100 events in our catalogue at all possible stations. We find that the correlation coefficient is higher than 0.7 if the second event occurs

D R A F T August 31, 2015, 4:43am D R A F T X -8
LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS after more than 0.3 s (Figure 2). At shorter inter-event times the correlation coefficient falls from 1.0 at zero inter-event time to around 0.65 at an inter-event time of 0.1 s. In reality, the correlation coefficient is higher than in these synthetics because aftershocks are usually smaller than their mainshock, hence produce less interference. Our selection criterion thus may exclude aftershocks at times shorter than 0.3 s. However, as shown later, the deconvolution results are quite noisy in the first 0.3 s anyway. Hence, this selection step does not limit the aftershock detection capability, it simply removes noisy waveforms.

The waveform u is filtered with a 4th-order Butterworth low-pass filter with cut-off frequency of 20 Hz. This stabilizes the deconvolution procedure, yielding less noisy RSTFs.

At each iteration we estimate

f n+1 = P [ f n + τ G T * (u -G * f n ) ] (3)
where f n is the RSTF estimated at iteration n, P a projection operator, G the EGF, G T its transpose, and τ is a relaxation parameter. Following [START_REF] Vallée | Stabilizing the Empirical Green Function Analysis: Development of the Projected Landweber Method[END_REF], we set

τ = 1/(max ω | Ĝ(ω)| 2 )
, where Ĝ is the Fourier transform of G and ω is frequency. The projection operator P is defined as follows:

Pf (t) = { f (t) if f (t) > 0 and 0 < t < T 0 elsewhere (4)
This imposes two constraints on f (t), positivity and an upper bound on its duration (T =20 s). The latter is mostly to prevent spurious values at the end of the time interval where data constraints are weak. The iteration process stops either when we reach a maximum number of allowed iterations or when the residual, ε =∥ u -G * f ∥ 2 does not decrease significantly any more. Once convergence is reached, we normalize f by its maximum

D R A F T August 31, 2015, 4:43am D R A F T LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS X -9
amplitude. We then stack RSTFs computed from all EGFs, at all stations (Figure 5).

By stacking RSTFs obtained at different stations we focus on aftershocks located close to their mainshock, otherwise stacking would not result in constructive interference. This is also consistent with the fact that our analysis relies on waveform similarity. The final product is, for each event in the original catalog, a source time function (STF) that can reveal the occurrence of similar events in the first ∼ 20 s.

Identifying Early Aftershocks

Our next objective is to identify if aftershocks occurred during the investigated time period, i.e. in the first 19.48 s after the mainshock. We scan through the stacked STF for prominent peaks. We first remove from the STF any long period trend. As the noise level might vary across the whole duration of the STF, we proceed to the detection of peaks in consecutive windows of 1 s. For each window we compute the mean µ and standard deviation σ of the STF function. If a peak (local maximum) is higher than µ + 5σ we keep it as a possible detection. We declare a detection if the peak is the largest one (excluding the first peak corresponding to the mainshock). We only found 1 instance of multiple aftershocks in our dataset. Following this procedure, we identified 68 early aftershocks out of 2414 events in the dataset. Repeating the same operation with a different threshold setting, m + 9M AD where m is the median and M AD the median absolute deviation of the STF within the 1 s window, we identify 78 early aftershocks. A total of 64 events with early aftershocks are identified in common by both criteria. The remaining events have large noise, so we exclude them from our further analysis.

D R A F T August 31, 2015, 4:43am D R A F T X -10
LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS

Verifying the Identification of Early Aftershocks

Although the waveforms considered here are very similar and tend to produce stable deconvolution results, in some instances the deconvolution process leads to noisy RSTFs, making it difficult to distinguish if a peak is significant. In order to assess the robustness of our RSTFs, we test a second deconvolution technique on the 64 events identified above.

As we expect only few aftershocks in the first 20 s, most of the RSTF values should be null. We hence employ a sparse deconvolution procedure that favors RSTFs with a low number of non-zero values, following a formulation similar to [START_REF] Rodriguez | Simultaneous recovery of origin time, hypocentre location and seismic moment tensor using sparse representation theory[END_REF]. We assume that the target waveform u can be written as a linear combination of a subset of the waveforms ϕ l , l =1,2, . . . , L, derived from the EGF waveform G by

ϕ l (t) = G(t -l.δt), (5)
i.e. the same Green's functions time-shifted to represent all the possible occurrence times of the aftershock. Here, L = 2048 is the number of samples of u and δt = 0.01 s is the sampling time interval. We build a matrix Φ with columns formed by the basis functions ϕ 1 , . . . , ϕ L . The target waveform is related to the RSTF s ∈ R L by

u = Φs. (6)
We hypothesize that the signal can be reconstructed by a small number m of basis functions, m ≪ L, that is, the RSTF has a maximum of m non-zero values. We express our deconvolution problem as

minimize ||r|| 2 = ||Φs -u|| 2 , ( 7
)
subject to:

Card{s ̸ = 0} ≤ m, (8) s > 0. (9) D R A F T August 31, 2015, 4:43am D R A F T X -11
The two constraints enforce sparsity and positivity, respectively. We solve the problem in Eqs. (7-9) by the Orthogonal Matching-Pursuit algorithm [START_REF] Tropp | Signal recovery from random measurements via orthogonal matching pursuit, Information Theory[END_REF].

The approach is similar to that of [START_REF] Kikuchi | Inversion of complex body wavesiii[END_REF]. We start from an initial null guess s = 0, residual r 0 = u and set the iteration counter to h = 1. We then identify the basis function that best matches the current residual,

n l = arg max l=1,L ϕ T l r h-1 . ( 10
)
We add the identified basis function to the set already obtained at the previous iterations,

Γ t = Γ h-1 ∪ n l , (Γ 0 = ∅).
We then solve a positive least squares problem, minimizing

||Φ(Γ h )s -u|| 2 , ( 11
)
under positivity constraint, s > 0. This yields a new estimate of the RSTF, s h , and the 6). We finally repeat the iteration procedure until h = m or until convergence (no significant update) is reached. Here we set m = 10.

updated residual r h = u h -Φ(Γ h )s h (Figure
Like for the Landweber deconvolution technique, we stack the RSTF obtained at different stations to obtain the STF of each event.

For most of the 64 analyzed events, the sparse deconvolution technique validates the STF obtained previously, i.e. it yields a very similar, but obviously more sparse STF (Figure 6). Ten STFs obtained by sparse deconvolution do not present any clear peak, and also correspond to weak peaks in the STF obtained by Landweber deconvolution.

The remaining 54 STFs show a clear secondary peak in both deconvolution methods and constitute our final set of early aftershocks. We finally extract from the identified peaks, the occurrence time ∆t of the aftershocks relative to their mainshock.

D R A F T August 31, 2015, 4:43am D R A F T X -12
LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS

Aftershock Locations and Magnitudes

Computing Time Delays

We are interested in locating the detected aftershocks relative to their mainshock. We already have some information. First, the high similarity between mainshock and aftershock waveforms suggests that the two events are very close to each other. Furthermore, the emergence of a peak in the STF by stacking the RSTFs of individual stations is only possible if the two events are close. The peaks in the STFs are clear and their width does not exceed 0.03 s (3 samples). To locate the aftershocks we first recover their waveforms. This is a difficult task as aftershock signals are strongly contaminated by the mainshock coda. We set the first 10 samples (0.1 s) of the RSTF to zero, then convolve it with the EGF. This effectively removes the contribution of the mainshock. This technique performs better than a simple subtraction of the scaled EGF waveform. We then isolate a 2.56 long window around the P-wave arrival of the aftershock and mainshock waveforms.

At each station, we determine the arrival time difference between mainshock and aftershock as the time lag that maximizes the correlation function of their waveforms. The search is restricted to time lags within 0.05 s of the relative interval ∆t of the stacked STF. Sub-sample precision differential times, dt, are estimated by quadratic interpolation of the correlation function. For each mainshock-aftershock pair we finally obtain a set of differential arrival times dt at each station, computed from all available EGFs.

Aftershock Relative Locations

Relative locations are here obtained by a cascaded Metropolis algorithm [START_REF] Mosegaard | Probabilistic Approach to Inverse Problems[END_REF]. We perform a random walk, that samples the posterior probability

dt calc (m) = 1 c p [δx sin(θ) sin(ϕ) + δy cos(θ) sin(ϕ) + δz cos(ϕ)] + δt 0 ; ( 12
)
where θ is the azimuth between the earthquake and the station, ϕ is the take off angle and c p is the P-wave velocity in the source area. The negative log-likelihood, l(m) for this new set of parameters is obtained as

l(m) = |dt calc -dt| σ dt ( 13
)
where σ dt is the uncertainty associated to a measurement of dt and is based on the correlation coefficient of the measured time delays. We adopt the L 1 norm to reduce the influence of outliers during the inversion procedure. The new set of parameters is accepted as a sample of the distribution ρ if l(m) < l 0 or if ln (X) < (l 0 -l(m)) with X a random number taken from a uniform distribution in [0-1] and l 0 the negative log-likelihood of the latest obtained sample of ρ. We then iterate the process and repeatedly draw another set of random parameters until we have at least 1000 samples of ρ.

We find that the mean horizontal distance between mainshocks and aftershocks is 56 remaining three events also have the least well resolved relative locations, with the highest uncertainties. For a robust analysis of the event locations we only keep those events for which the location uncertainties in both δx and δy (measured as the standard deviation of the distributions) are smaller than 20 m (Figure 7). This selected subset consists of 38 events out of the initial 54. The aftershocks that have poorer location accuracy are mostly associated with the lowest peaks in the STF, corresponding to aftershocks of relatively small magnitude compared to their mainshock. The relative locations of the majority of the well-located aftershocks appear to be in agreement with the inferred mean orientation of the San-Andreas fault near Parkfield, N150E [START_REF] Thurber | Three-Dimensional Compressional Wavespeed Model, Earthquake Relocations, and Focal Mechanisms for the Parkfield, California, Region[END_REF] (Figure 7). All selected aftershocks occur close to their mainshock, confirming that the identified peaks in the STF are truly associated to aftershocks. We estimated mainshock rupture dimensions from catalog magnitude converted to moment using the relation of [START_REF] Bakun | Seismic moments, local magnitudes, and coda-duration magnitudes for earthquakes in central california[END_REF] and assuming a circular rupture with constant stress drop of 1 MPa. The well-located aftershocks are located within 1 to 2 rupture lengths of their mainshock, but it is difficult to assess if they occur within the mainshock rupture area or at its edge as observed by [START_REF] Rubin | Aftershock asymmetry/rupture directivity among central San Andreas fault microearthquakes[END_REF].

We also observe an asymmetry of aftershock locations, with 26 out of 38 (2/3) aftershocks occurring to the NW of their mainshock. We quantify the significance of the observed asymmetry by testing against the null hypothesis that the aftershock relative locations result from a Bernoulli trial with two outcomes of equal probability: NW or SE location. Under the null hypothesis, the probability, P , of observing at least 26 out of 38 events in a preferred direction is given by where p = q = 0.5 and n = 38. This gives P = 3.5 %. It is thus quite unlikely that the observed asymmetry is the result of random drawing from an equal outcome binomial distribution. This supports the observed asymmetry as a real feature of the early aftershock seismicity.

P = 2 38 ∑ k=26 ( n k ) p k q n-k (14) D R A F T August

Aftershock Magnitudes

We estimate the magnitude of the newly detected aftershocks. We proceed by computing the amplitude of the aftershock peak in the STF relative to the mainshock peak, A AM . To account for the width of STF peaks, we define A AM based on the sum the STF amplitude over 3 time samples (0.03 s) around a peak. The relative STF amplitudes are equivalent to a relative moment. We convert them to a magnitude difference following [START_REF] Bakun | Seismic moments, local magnitudes, and coda-duration magnitudes for earthquakes in central california[END_REF]'s relation for earthquakes in Central California. The aftershock magnitude M A is

M A = M M + 1 1.2 log 10 (A AM ) (15)
where M M is the mainshock magnitude. We repeat the operation at all stations and for all possible EGFs, then average to obtain our final estimate of M A for each aftershock. We obtain an average value of ⟨M A ⟩ = 1.05 for all 54 identified aftershocks. In comparison, the mean magnitude for the 2414 events in the RES catalog is 1.52. Hence our method detects events of lower magnitude than in the original catalog.

Temporal Distribution of Sub-Events

We analyze the temporal organization of the detected early aftershocks and compare it to that of later aftershocks listed in the original catalog (> 20 s). This comparison accounts for differences in detection threshold between the two catalogs. On the one hand, early aftershocks are more difficult to detect due to interference by the mainshock coda. On the other hand, our method detects smaller earthquakes than in the original catalog.

To account for changes of detection capability across time scales we first compute the magnitude distribution of events listed in the original catalog. We observe that the magnitude distribution is well fitted by a model that combines a probability function of detecting earthquakes and the Gutenberg-Richter distribution [START_REF] Ogata | Analysis of temporal and spatial heterogeneity of magnitude frequency distribution inferred from earthquake catalogues[END_REF]. The best fit to this model gives b = 1.0 for the Gutenberg Richter law. The estimated magnitude of completeness is m c = 1.2 (Figure 8). We assume that the magnitude distribution of the 54 newly detected events can be described by a similar model with b = 1.0 but with a different value of m c . Due to the low number of events for this new population of events, we impose the value of b to reduce the uncertainty on the estimation of m c . We found by fitting the distribution to the model of [START_REF] Ogata | Analysis of temporal and spatial heterogeneity of magnitude frequency distribution inferred from earthquake catalogues[END_REF] that m c = 0.8. The difference of the magnitude of completeness of the two catalogs is ∆m c = 1.2 -0.8 = 0.4.

In order to compare the rates of aftershocks in the two catalogs, we correct the aftershock rate derived from the original catalog by the factor α = 10 ∆mc = 2.5. We then build a composite aftershock sequence. We treat each event in the repeating earthquake catalog as a potential mainshock and all subsequent events of the same sequence as aftershocks.

Based on the time delay between events in a given sequence we then compute the rate of earthquakes following a each mainshock. We then stack results obtained for all sequences to obtain a composite aftershock sequence. In a second step we add to the already existing repeating sequences the new events detected with our deconvolution technique. We only factor α to account for the difference of completeness of the two catalogs. The resulting evolutions of aftershock rate are shown in figure 9. The early aftershock rate implied by the new events we detected agrees with the extrapolation to early times of the Omoritype behavior deduced from the original catalog. Notably, we resolve the emergence of the power-law decay starting at ∼ 0.3 s after the mainshocks, that is, a time scale one order-of-magnitude longer that the mainshock rupture durations. The earthquake rate several years after a repeating earthquake deviates from Omori's law because repeating micro-earthquakes at Parkfield occur quasi-periodically with a recurrence time of about a year.

We tested, through synthetic tests, if the detection threshold of our method depends on the time since mainshock and on the mainshock magnitude. The tests are presented in Appendix B. We find that the detection threshold does not vary with time since mainshock. We also find that the changes of detectability related to the mainshock magnitude are not very large and when taken into account give similar results as the one presented in figure 9.

Discussion

Our detection method recovers early aftershocks of micro-earthquakes (M L = 0.9 to Most previous studies of early aftershock rates focused on a single aftershock sequence, notably following large magnitude events. To compare our study with these others studies, we normalize the earliest resolvable time of the aftershock rate by the mainshock rupture duration. Considering a typical dimension of 100 m for the events considered in our study and a rupture speed of 3 km/s, the typical rupture duration is 33 ms. This is 10 times shorter than the time scale of emergence of power-law decay found here, ∼ 0.3s. This factor of 10 separation is similar, for example, to the one resolved for the 2004 M w 6

Parkfield earthquake: [START_REF] Peng | Anomalous early aftershock decay rate of the 2004 Mw6. 0 Parkfield, California, earthquake[END_REF] resolved aftershock rates starting 130 s after the mainshock onset and the estimated mainshock rupture duration is ∼ 6 s Fletcher et al.

[ 2006]. Despite this similarity of relative time scales, most studies of early aftershocks of large events infer a departure from Omori's power-law decay at early times, i.e. a slower decay rate just after the mainshock than predicted by extrapolating the later Omoritype decay [e.g. [START_REF] Peng | Anomalous early aftershock decay rate of the 2004 Mw6. 0 Parkfield, California, earthquake[END_REF]. In contrast, our results show that the aftershock rate shortly after the end of the mainshock rupture decays according to Omori's power-law, i.e. no flattening is observed down to the lowest resolvable time scale (figure 9). Multiple explanations of such a difference can be proposed: i) repeating earthquake sequences, which are the specific focus of our study, might have different aftershock sequences than the rest of the seismicity; ii) aftershock sequences may behave differently for large and small magnitude mainshocks; iii) the magnitude correction applied in large-mainshock studies is different from the one used in this study because of the severe effect of mainshock coda, and may result in underestimation of the early aftershock rates in previous studies.

The constant aftershock decay exponent we resolve over the entire time span suggests that micro-earthquake triggering results from the same physical processes from scales of

D R A F T August 31, 2015, 4:43am D R A F T
LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS X -19 0.3 s up to 100 days. It is naturally conceivable that at even earlier times aftershock activity is controlled by different triggering processes related to dynamic stresses carried by mainshock waves. This is supported by the analysis of [START_REF] Wang | Compound earthquakes on a bimaterial interface and implications for rupture mechanics[END_REF] to resolve aftershocks occurring at times shorter than 0.2 s and located close to their mainshocks.

They detected events at these short time scales through a parametric inversion procedure applied to sub-sampled waveforms.

The relative timing and location of these very early aftershocks are consistent with the passage of shear waves radiated by the mainshock.

Wang et al.

[2014] detected overall 153 aftershocks within a distance of 2 mainshock rupture radii, out of a total of 20990 possible mainshocks. They also found that 100 out of these 153 events occurred near the mainshock S-wave front and their relative locations showed a pronounced asymmetry, favoring aftershock triggering to the SE of their mainshocks. This corresponds to the preferred rupture direction predicted by theory and simulation of dynamic rupture on bimaterial faults [START_REF] Rubin | Aftershock asymmetry on a bimaterial interface[END_REF][START_REF] Ampuero | Cracks, pulses and macroscopic asymmetry of dynamic rupture on a bimaterial interface with velocity-weakening friction[END_REF], consistent with rupture directivity observations [START_REF] Wang | Rupture directivity of microearthquakes on the San Andreas Fault from spectral ratio inversion[END_REF], and the bimaterial structure in the Parkfield area [START_REF] Mcguire | High-resolution imaging of the bear valley section of the san andreas fault at seismogenic depths with fault-zone head waves and relocated seismicity[END_REF][START_REF] Lengliné | Rupture directivity of microearthquake sequences near Parkfield, California[END_REF][START_REF] Kane | Rupture directivity of small earthquakes at Parkfield[END_REF]. Here we found more aftershocks located in the opposite direction (NW). We investigate if this asymmetry changes with time since the mainshock. We compute the ratio between the number of NW events and the number of SE events that occurred since the mainshock time for various time intervals. We find that this asymmetry is largest at the earliest times in our study, from 0.3 s to about 1 min, and then decreases progressively (figure 9-bottom). Such aftershock asymmetry, opposite to the preferred rupture direction, has been previously observed over long time

D R A F T August 31, 2015, 4:43am D R A F T X -20
LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS scales of 10 s to 9 hrs [START_REF] Rubin | Aftershock asymmetry/rupture directivity among central San Andreas fault microearthquakes[END_REF][START_REF] Rubin | Aftershocks of microearthquakes as probes of the mechanics of rupture[END_REF], and up to two days [START_REF] Zaliapin | Asymmetric distribution of aftershocks on large faults in california[END_REF]. The earthquakes resolved here hence belong to this population of longer-term aftershocks, rather than to the very early events triggered dynamically. [START_REF] Wang | Compound earthquakes on a bimaterial interface and implications for rupture mechanics[END_REF] proposed that the excess of aftershocks to the SE at very early times could account for the events missing in the SE in the longer term (10 s-9 hours) aftershock population, if sites are allowed to break only once during the whole aftershock sequence.

Our results imply that the shift of the direction of aftershock asymmetry occurs less than 0.3 s after the mainshock. Aftershocks in the preferred rupture direction on bimaterial faults (here, SE) were proposed to result from a tensile stress pulse propagating with the rupture front of the mainshock in the preferred rupture direction [START_REF] Rubin | Aftershock asymmetry on a bimaterial interface[END_REF][START_REF] Wang | Compound earthquakes on a bimaterial interface and implications for rupture mechanics[END_REF]. Aftershocks at longer time scales (here, 0.3 s 100 days) result from a different mechanism.

The early aftershocks identified in our study may result from static stress transfer directly from the mainshock rupture or indirectly from its afterslip. Afterslip can reproduce the 1/t aftershock rate decay, and is constrained by our results to start earlier than 10 times the mainshock rupture duration [START_REF] Helmstetter | Afterslip and aftershocks in the rate-and-state friction law[END_REF]. Aftershock migration is predicted by models where aftershocks are driven by an expanding afterslip front [e.g. [START_REF] Kato | Expansion of aftershock areas caused by propagating post-seismic sliding[END_REF]. The persistence of aftershock asymmetry on bimaterial faults up to long times remains unexplained by current numerical models [START_REF] Rubin | Aftershock asymmetry on a bimaterial interface[END_REF],

but the possibility of asymmetric afterslip has not been explored. While we do not observe aftershock migration in our results, resolving it would require the detection of more distant aftershocks, but this is limited by our selection criteria and our focus on similar In the context of rate-and-state friction models of aftershocks triggered by coseismic static stress steps, the onset time c of power-law aftershock rate decay is related to frictional parameters [START_REF] Dieterich | A constitutive law for rate of earthquake production and its application to earthquake clustering[END_REF] by

c = t a exp ( -∆τ aσ ) , ( 16
)
where t a is the aftershock duration, in the range 0.5 to 5 years in this region [START_REF] Toda | Response of the San Andreas fault to the 1983 Coalinga-Nuñez earthquakes: An application of interaction-based probabilities for Parkfield[END_REF], ∆τ is the average shear stress change caused by the mainshock rupture, a is a rate-and-state friction parameter quantifying the importance of the immediate logarithmic velocity-strengthening effect, and σ is the effective normal stress. Considering c < 0.3 s and the lower bound t a = 0.5 years, we estimate ∆τ /aσ > 20. This lower bound on the stress stimulus is consistent with the idea that the stress transferred to the immediate vicinity of a rupture is higher than the stress drop within the rupture, which scales in the rate-and-state friction model as (b -a)/σ log(V dyn /V load ) ∼ 20(b -a)σ [e.g. [START_REF] Perfettini | Dynamics of a velocity strengthening fault region: Implications for slow earthquakes and postseismic slip[END_REF]. Most of the aftershocks we resolve are located near the edge of their mainshock rupture but we acknowledge important relative location uncertainties, notably because we estimate mainshock sizes assuming a constant stress drop for all events. High stress concentrations triggering aftershocks within the nominal mainshock rupture area are also possible if its slip distribution is heterogeneous [e.g. [START_REF] Herrero | A kinematic self-similar rupture process for earthquakes[END_REF][START_REF] Mai | A spatial random field model to characterize complexity in earthquake slip[END_REF]. Considering the regional upper bound of earthquake duration, that high are not expected to prevail over distances larger than one mainshock radius, but may be accounted by stress transfer via afterslip.

Conclusion

By analysing precisely early aftershocks of similar micro-earthquakes in Parkfield, California, we extend the resolution of aftershock rates down to 0.3 s after the mainshock origin time, that is, ∼ 10 times the mainshock rupture duration. Over a time scale span of nearly 8 orders of magnitude, from 0.3 s up to more than 100 days, the aftershock rate decay is well described by a single Omori power-law with no flattening at early times.

If a characteristic time for the onset of the power-law regime exists, it is necessarily shorter than 0.3 s. Our results suggest that aftershocks occurring beyond the time scales of dynamic triggering arise from relatively large stress perturbations, possibly caused by aseismic afterslip with very short characteristic time. We also observe an asymmetry of aftershock relative locations along-strike, persistent over long periods and consistent with previous observations on bimaterial faults.

Appendix A: Comparison with the Template-Matching Approach

We consider the example shown in Figure 4. We compute the cross-correlation function between the signal with sub-event and the signal without sub-event (from another earthquake of the same repeating sequence). We also compute for reference the auto-correlation function of both signals. A second peak of the correlation function appears at the time of the sub-event (Figure 10) and is about twice as high as other nearby secondary peaks.

In contrast, the peak associated with the sub-event in the RSTF obtained by Landweber deconvolution (also shown in Figure 10) is about one order of magnitude higher than all D R A F T August 31, 2015, 4:43am D R A F T

  density distributions of the model parameters ρ(m), where m = {δx, δy, δz, δt 0 } is the D R A F T August 31, 2015, 4:43am D R A F T LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS X -13 vector containing the model parameters, i.e. the relative 3D position and origin time between mainshock and aftershock. The data are given by the computed time delays, d = dt. The initial guess is δx = 0, δy = 0, δz = 0, δt 0 = ⟨dt⟩, i.e. co-located with the mainshock. We then draw a set of random parameters around the initial guess. The new parameters are drawn from a Gaussian distribution centered on the parameters estimated at the previous step and with a standard deviation of 100 m for relative distances and 0.1 s for relative times. Given this new set of parameters we compute the expected time delays, dt calc ,

  m. Out of the 54 events, 51 have a horizontal separation shorter than 100 m. The D R A F T August 31, 2015, 4:43am D R A F T X -14 LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS

  31, 2015, 4:43am D R A F T LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS X -15

  consider earthquakes occurring before the 2004 Parkfield M w 6 earthquake, to avoid the influence of the major stress perturbations it caused. This selection reduces the number of new events to 44. The aftershock rate computed from catalog data is multiplied by the D R A F T August 31, 2015, 4:43am D R A F T LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS X -17

  3.2), previously missed because hidden in the mainshock signals, down to ∼ 0.3 s after the mainshock initiation. Our method is currently limited at short time-scales by noise and frequency content of the waveform data. The resolvable mainshock-aftershock time separation is limited by the width of the source time functions retrieved by deconvolution, which depends on noise level and dominant waveform period, in turn controlled by attenuation or rupture duration and filter high-frequency cutoff. D R A F T August 31, 2015, 4:43am D R A F T X -18 LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS

  events. Thus we cannot rule out the afterslip-driven model on the basis of lack of observed D R A F T August 31, 2015, 4:43am D R A F T LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS X -21 aftershock migration. The asymmetric residual stresses left by the mainshock rupture on a bimaterial fault Rubin and Ampuero [2007] could contribute to asymmetry of afterslip, even without migration.

t

  a = 5 years, would lead to even larger estimates of stress transfer. Coseismic stresses D R A F T August 31, 2015, 4:43am D R A F T X -22 LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS

Figure 1 .

 1 Figure 1.Map of the area, seismic stations (triangles) and earthquake epicenters

Figure 2 .

 2 Figure 2. Top : Waveform correlation coefficient between a single event and a synthetic

Figure 3 . 31 Figure 4 .

 3314 Figure 3. Waveforms of a repeating earthquake sequence comprising 9 events, ordered

Figure 5 .Figure 6 .

 56 Figure 5.Result of EGF deconvolution for an earthquake of a repeating sequence.

Figure 7 .

 7 Figure 7. Map view of the locations of aftershock centroids relative to their mainshock.

Figure 8 . 33 Figure 9 .

 8339 Figure 8. Magnitude distribution for all events included in the repeating earthquake

Figure 10 .Figure 11 .

 1011 Figure 10. Red and blue curves represent the autocorrelation function of the signal

Figure 12 .

 12 Figure 12. Same as figure 9 but with rate corrections applied to the newly detected

  This task is non-trivial because analysis of early aftershock signals is severely obstructed by the coda of the mainshock. It is not yet clear if the characteristic time c

	X -4	LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS
	mechanism.

Acknowledgments. Data used in this study are from the Northern California Seismic Network, U.S. Geological Survey, Menlo Park (http://www.ncedc.org/ncsn/). We thank A. Rubin for numerous comments and Z. Duputel for discussion on the Metropolis Algorithm. We thank the associate editor and reviewers Y. Kagan and B. Enescu for suggestions. We acknowledge funding from NSF (grant EAR-1015698).

X -23

other RSTF peaks. This example illustrates that, when the aftershock signal is dominated by the mainshock signal, the template matching approach has a poorer sub-event detection capability than our deconvolution method.

Appendix B: Synthetic Tests on Detection Thresholds

We perform synthetic tests to estimate the detection capability of our deconvolution method. We select randomly 50 events in the repeating earthquakes catalog. We ensure that none of these events contains a sub-event in the first 20 s. We add a duplicate of the mainshock signals with various time delays and amplitudes relative to the mainshock.

Relative time delays are set at 0.05, 0.1, 0.2, 0.3, 0.5, 1, 2, 3, 5, 10, and 15 s. Relative maximum amplitudes are set to 0.005, 0.01, 0.03, 0.1, 0.3, and 0.5. We then run our deconvolution sub-event detection process. We first find that, for all events, the interevent time has no influence on the detection capability as long as the sub-event occurs later than 0.2 s after the mainshock (see Figure 11).

The minimum relative amplitude of the sub-event that is detected varies with the mainshock magnitude. For a magnitude 1.2 event we can detect sub-events down to a relative amplitude of 0.1, which corresponds to a minimum magnitude of 1.2 -1 1.2 log 10 (0.1) = 0.4.

For a larger mainshock (m=3.0) detection is possible down relative amplitude of 0.01, or a sub-event magnitude as low as 3.0 + 1 1.2 log 10 (0.01) = 1.3. The minimum resolvable magnitude evolves between these two values as a function of mainshock magnitudes, although the scatter is important. We define the magnitude detection threshold for a mainshock of magnitude m based on a linear fit of the detection thresholds obtained for the magnitude 1.2 and 3.0 events. We then corrected our aftershock rate, for the newly detected events, taking into account this change of detection threshold. We give a weight to each detected

LENGLIN É AND AMPUERO: REPEATING AFTERSHOCKS earthquake that depends on its probability of being detected (based on the mainshock magnitude). This results in the new aftershock rate shown in Figure 12. We note that this correction leaves unchanged the observation of a continuous decay of the aftershock rate over the whole considered period.