P Poli

H A Pedersen

M Campillo

Polenet / Lapnet

Emergence of body waves from cross-correlation of short period seismic noise

Keywords: Interferometry, Body waves, Wave propagation

Ambient noise correlation is now widely used in seismology to obtain the surface waves part of the Green's function. More difficult is the extraction of body waves from noise correlations. Using 42 temporary broad-band three components stations located on the northern part of fennoscandian region, we identify high frequency (0.5-2 Hz) body waves emerging from noise correlations for inter-station distances up to 550 km. The comparison of the noise correlations with earthquake data confirm that the observed waves can be interpreted as P and S waves reflected from the Moho. Because the crustal model of the area is well known, we also compared the noise correlations with synthetic seismograms, and found an excellent agreement between the travel times of all the observed phases. Polarization analysis provide a further arguments to confirm the observation of body waves.

INTRODUCTION

The possible extraction of the GreenÕs function through correlation of seismic noise has opened up for potentially new and exciting developments in seismic imaging and monitoring of the elastic properties in the Earth. The feasibility of the method is now understood through a series of theoretical developments (e.g. [START_REF] Weaver | Ultrasonics without a source: Thermal !uctuation correlations at MHz frequencies[END_REF][START_REF] Wapenaar | Retrieving the elastodynamic greenÕs function of an arbitrary inhomogeneous medium by cross correlation[END_REF]Roux et al. 2005b[START_REF] Gouždard | Cross-correlation of random "elds: mathematical approach and applications[END_REF][START_REF] De Verdižre | A semi-classical calculus of the correlations[END_REF] and through laboratory experiments [START_REF] Weaver | Ultrasonics without a source: Thermal !uctuation correlations at MHz frequencies[END_REF]. Shapiro and Campillo demonstrated the feasibility of the method by extracting intermediate and long period surface waves on field data [START_REF] Shapiro | Emergence of broadband rayleigh waves from correlations of the ambient seismic noise[END_REF], and numerous studies have now based imaging on the analysis of seismic surface waves extracted by noise correlations (e.g. Sabra et al., 2005b;[START_REF] Shapiro | High-Resolution Surface-Wave tomography from ambient seismic noise[END_REF][START_REF] Yang | Ambient noise rayleigh wave tomography across europe[END_REF][START_REF] Yao | Surface wave array tomography in SE tibet from ambient seismic noise and two-station analysis -II. crustal and upper-mantle structure[END_REF][START_REF] Ritzwoller | Ambient noise tomography with a large seismic array[END_REF].

It is generally assumed that noise is related to surface activity, ranging from human activity at high frequency to the forcing of oceans and atmosphere at low frequency. In the absence of deep sources, and with uneven distribution of surface sources, the reconstruction of body waves relies on the energy that has been scattered at depth. Although with an energy smaller than the one of the surface waves locally radiated by the sources, the scattered body waves are present in actual seismograms acquired at the surface, as a part of the almost equipartitioned diffuse field observed for long lapse time (e.g [START_REF] Hennino | Observation of equipartition of seismic waves[END_REF][START_REF] Campillo | Phase and Correlation in `Random' Seismic Fields and the Reconstruction of the Green Function[END_REF].

There should therefore be substantial hope of extracting the body wave part of the GreenÕs function, albeit with a lower signal to noise ration than the dominant surface waves.

Body waves have indeed been reported from short distance range correlations. Roux (2005a) identified direct P waves from noise correlation, using data from a small array in California.

Their noise derived P waves were linearly polarized, and with a velocity compatible with a known velocity model of the area. [START_REF] Draganov | Re ection images from ambient seismic noise[END_REF] used data from oil exploration to extract reflected P waves from shallow interfaces. Their observed body waves were in good agreement with the active source reflection response in the same area. [START_REF] Zhan | Retrieval of moho-re!ected shear wave arrivals from ambient seismic noise[END_REF] identified S reflected phases from the Moho interface at the critical distance in two shield areas. The S waves presented a striking agreement with earthquake data. On the contrary, the extraction of body waves over broader distance ranges has so far not been successful, even though such waves would be key for body wave tomography at crustal scale.

To study the possibility of recovering body waves over large distance ranges from seismic noise correlation, we processed one year of data acquired at POLENET/LAPNET seismological array [START_REF] Kozlovskaya | POLENET/LAPNET-a multi-disciplinary geophysical experiment in northern fennoscandia during IPY 2007-2008[END_REF]. This dataset is acquired by a seismic array including 42 broadband stations. The study area is part of the Precambrian northwestern segment of the East European Craton, and the crust is relatively well known from active seismic experiments (HUKKA, FIRE, FIRE4, POLAR), from which is it known that the velocity structure remains relatively simple and with limited lateral variations and a limited variation of Moho depth (Janik et al., 2007, and reference therein). The major seismic phases observed from active source experiments were PmP and SmS [START_REF] Janik | Crust-mantle boundary in the central fennoscandian shield: Constraints from wide-angle p and s wave velocity models and new results of re!ection pro"ling in Finland[END_REF], while weak amplitudes were reported for mantle phases (Pn and Sn) and inter-crustal reflection (Pg and Sg). As body waves are expected to be strong and impulsive in a crust characterised by weak scattering and attenuation, the geological and such observations of strong PmP combined with weak crustal scattering and attenuation [START_REF] Pedersen | Depth dependance of q beneath the baltic shield inferred from modeling of short period seismograms[END_REF][START_REF] Uski | A new local magnitude scale for the Finnish seismic network[END_REF], are particularly promising elements for the extraction of body waves from noise correlations.

We firstly present the dataset and the data processing, after which we discuss the noise correlations and the fast travelling waves that we interpret as bodywaves. This interpretation is supported by records of local seismic events, and by numerical simulations (arrival times, polarization) in a crustal model which is derived from the results of the active seismic studies.

DATA AND SIGNAL PROCESSING

We analyse three component seismic data continuously recorded during the POLENET/LAPNET temporary experiment. We used only stations equipped with broadband sensors, and included several permanent broadband stations in our data-set. The array configuration (fig. 1) is approximately a 2D grid with station separations that span from ~50 km to ~600 km. The array was installed between spring and autumn 2007, for a duration of two years. We used records for the calendar year 2008 during which the array was fully operational. [START_REF] Pedersen | Influence of the seismic noise characteristics on noise correlations in the Baltic shield[END_REF] reported the presence of strong directivity of the noise field, especially observed from strong asymmetric surface waves signals on the noise correlations at intermediate frequencies (0.02-0.1 Hz) while the high frequency (0.1-1 Hz) part of the noise was distributed over larger azimuth ranges. Their study was limited to winter month, and the major part of high frequency energy was related with the sea activity along the eastern Atlantic coast. In our case, where the seismic array was installed for two years, the average of the correlations over one year (2008) in a frequency range from 0.1 and 2 Hz present an high signal to noise ratio over long distances (~ 600 km) and in both causal and acausal parts of the correlations as expected in a fully diffuse wave-field or in presence of randomly distributed sources.

The standard pre-processing included removing the data mean and trend, prefiltering, resampling to identical sampling rate and deconvolution of the instrumental responses. The noise correlations were calculated for all combinations of radial, transverse and vertical components, which required rotation of the horizontal components for each station pair according to the azimuth at each station of the great-circle between the two stations. To analyse broadband signals while removing the effects of earthquakes, we applied two supplementary steps before correlation. We firstly split the continuous data into four hour windows and removed the ones where amplitudes were present which were larger than 3 times the standard deviation of the signal. This step additionally reduces the effect of instrumental problems such as spikes. Secondly, we apply a spectral whitening in a frequency band from 0.1 to 2 Hz. This second step also diminishes the relative predominant contribution of surface waves related to the secondary microseismic peak at ~ 0.14 Hz.

After this processing the seismic noise traces are correlated for all couple of stations and stacked over one year, without applying 1-bit normalization. We verified the quality of the correlations by comparing different processing procedures, including one where all major earthquakes were removed. In this test we used the correlations processed as described above.

Then, using the earthquakes and explosions database of the Finnish seismological service, we removed all the time windows where earthquakes or explosions occurred and we stacked the correlations over one year. Since we study high frequency noise data (0.1-2 Hz), local seismicity can dramatically reduce the quality of the noise correlations [START_REF] Bensen | Processing seismic ambient noise data to obtain reliable broad-band surface wave dispersion measurment[END_REF].

We compared the correlations stack with and without local seismic events and we observed a stable reconstruction of all the seismic phases. From this observation we consider our processing procedure as sound and not significantly contaminated by quarry blasts and seismic events.

The bandpass filtered (0.5-2 Hz) noise correlations that have a signal to noise ratio larger than five are show in Figures 2, 3 and 4. The correlations are organized so that positive times correspond to waves propagating from the easternmost to the westernmost of the two stations. Out of the nine components of the correlations, we show the vertical-vertical (ZZ), radial-radial (RR) and transverse-transverse (TT) components, plotted as a function of the inter-station distance.

RESULTS AND DISCUSSION

In all the analysed components of the correlations we can identify coherent surface waves propagating from one station to another. The fundamental mode Rayleigh waves (Rg) portion of the estimated Green's function (EGF) are observed on both ZZ and RR correlation (fig. 2 and 3 respectively), with a propagation velocity of approximately 3 km/s. On the TT component (fig. 4) of the EGF we observe fundamental mode Love waves (L), with a velocity of ~ 3.5 km/s. Both Love and Rayleigh waves appear symmetrically on the seismic sections, which indicates either a good diffusivity of the noise-field or well distributed noise sources.

The signal to noise ratio of these high-frequency surface waves remains high out to the full distance range covered by the array, i.e. over up to 600 km.

We here wish to draw attention to other coherent phases that are clearly present in the seismic sections. Firstly, we note the coherent phase which is present on the ZZ components of the correlations (fig. 2) with an apparent velocity of approximately 3.5 km/s, which corresponds to expected apparent velocities for SmS waves, i.e. waves reflected at the Moho.

These waves are stronger on the acausal part of the correlations, so they must originate in a different source or scatter distribution than the fundamental mode surface waves discussed above. This type of wave is not observed on the TT component, however such waves would be masked by the Love waves which also have velocities close to 3.5 km/s. Secondly, a signal with an apparent velocity of approximately 6 km/s, i.e. close to the expected apparent velocity of the PmP phase, is observed on the acausal part of the RR component. Frequency-time analysis shows that these two phases are non-dispersive over the frequency interval where they can be observed, which is 0.5-2 Hz. We can therefore, at this stage, hypothesize that these waves present in the noise correlations are body waves, and most likely SmS and PmP waves.

A first verification of whether the high velocity signals on the correlations are consistent with being body waves, we compare the noise correlations with earthquake data.

We use the acausal part of the correlation traces of which we use only the ones with a signal to noise ratio higher than ten (in the body wave arrival windows). We choose a shallow local event (M L =2.9) located beneath the northern part of the array (fig. 1) and for which clear signals are observed for the frequency band of interest (0.5-2 Hz). The earthquake data are preprocessed identically of the continuous noise recordings, and the horizontal components are rotated to obtain the radial and transverse components.

Figure 5 shows the seismic sections with the radial (5a, 5b) and vertical (5c, 5d) components of the noise correlations (blue) and earthquake records (black). The distance axis corresponds to the inter-station distance for the noise correlations and the epicentral distance for the earthquake records. The earthquake data clearly show the fundamental mode Rayleigh waves on both radial and vertical components and the faster P and S wave on the radial and vertical components respectively. For the earthquake data, both single SmS and PmP are emerging from a distance of approximately 110km, which corresponds to approximately critical distance, however their amplitude is high from approximately 200km distance. SmS 2 is clearly observed from 280 km distance.

The comparison with earthquake data seems to give further evidence that the observed waves that arrive prior to the surface waves could indeed be body waves. Because the crustal model of the area is well known, and as the crustal model only varies very moderately beneath the study area, we can also directly compare the noise correlations with the numerical GreenÕs functions calculated using a 1-D Earth model. We base our velocity model (see table 1) on the interpretation of HUKKA seismic reflection profiles as presented by [START_REF] Janik | Crust-mantle boundary in the central fennoscandian shield: Constraints from wide-angle p and s wave velocity models and new results of re!ection pro"ling in Finland[END_REF]. The velocities of the upper crust are modified using the parameters obtained by [START_REF] Pedersen | Depth dependance of q beneath the baltic shield inferred from modeling of short period seismograms[END_REF] who analysed high frequency Rayleigh waves from a quarry blast to obtain shear velocities and quality factor Q down to 3km depth. At larger depths we used Q values based on [START_REF] Uski | A new local magnitude scale for the Finnish seismic network[END_REF]. We calculate synthetic seismograms using the frequency-wavenumber method by [START_REF] Bouchon | A simple method to calculate greenÕs functions for elastic layered media[END_REF] The velocity is similar to the one of the early waves in the Z-Z noise correlations, and the similarity is striking as to the pattern where the SmS 2 phase gradually become dominant over the SmS phase from 350km and beyond. Weak P waves can also be observed on the vertical component synthetic seismograms, with a velocity of approximately 6 km/s, as the ones observed on the Z-Z correlations. The signal to noise ratio on the correlations is however insufficient to easily detect them over the whole of the section. More evident are the P waves, observed on the RR correlations, that propagate with a velocity of ~ 6 km/s which is close to the velocity of the PmP phases observed on the RR synthetics seismograms.

The very successful comparison of the noise correlation sections with the earthquake records and the synthetic seismograms are strong arguments in favor of explaining the early arrivals in the noise correlations as body waves. A final argument resides in a strong similarity in polarizations of synthetic seismograms and noise correlations. The analysis of polarized seismic waves requires phase and amplitude informations of the seismic traces.

Strong non linear pre-processing (as applied for noise correlation) can be a limitation, because their effect on the amplitude of the signals. Recent results [START_REF] Cupillard | The one-bit noise correlation: a theory based on the concepts of coherent and incoherent noise[END_REF][START_REF] Prieto | On the amplitude information carried by ambient seismic field[END_REF] demonstrate that standard pre-processing as one-bit or whitening have little effect on the amplitude informations of the noise correlation functions, so that attenuation can be obtained from the EGF [START_REF] Prieto | On the amplitude information carried by ambient seismic field[END_REF]. From the previously cited works emerge that is possible perform polarization analysis using ambient noise, also if pre-processing was applied to the raw data.

Figure 7 b-c shows the particle motion observed for a correlation chosen for its high signal to noise ratio for a station couple located sufficiently far apart (211 km, station pair KIF-LP51) to separately analyse the participle motion of the Rayleigh waves and the two early hypothesized body waves. We used the ZR and ZZ components of the noise correlations to obtain their particle motion and compare it with the one computed for a vertical force acting onto the EarthÕs surface. The particle motions are shown in three time windows, which correspond to the Rayleigh wave and the two hypothesized body waves.

The agreement between particle motion as observed on synthetic seismograms and noise correlations is striking. The Rayleigh waves have, as expected, an elliptical particle motion with very similar ratio between the ZR and ZZ axis. The PmP wave is linearly polarized, with a coefficient of rectilineartiy ~ 0.9 (~ 1 on the synthetic polarization), and polarization angle of ~ 52¡ as compared to vertical, which is very similar to the ~ 56¡ observed on the synthetic motion. The SmS polarization is more complex due to free surface conversion. A linearly polarized SV wave incident at free surface beyond the critical conversion angle create phase shifted reflected SV wave and an evanescent P wave (Aki and Richards, 1980). The result of this sum of differently polarized waves can be observed in the synthetic seismograms as an inclined, elongated elliptic like polarization. Remarkably, the polarization on the noise correlations is in very good agreement with the synthetics also for these waves. This is a strong argument in favour of our interpretation of these waves as SmS.

CONCLUSION

The noise correlations from the northernmost part of the Baltic Shield are dominated by fundamental mode Rayleigh and Love waves. In addition to these waves, we observe coherent phases up to 500km inter-station distance which have an apparent velocity higher than the one observed for the fundamental mode surface waves. The fast waves are relatively high frequency (0.5-2Hz) and non-dispersive in that frequency range.

The comparison with earthquake records from a local event also showed the presence of these waves, and synthetic seismograms were also in excellent agreement with the noise correlation. The synthetic seismograms very clearly points towards identifying the fast waves as single and multiply Moho reflected P and S waves, an interpretation which is supported by the wave polarization for different time windows. Note that the agreement between noise correlations and synthetic seismograms was dependent on the use of a crustal model based on active seismic studies (e.g. [START_REF] Janik | Crust-mantle boundary in the central fennoscandian shield: Constraints from wide-angle p and s wave velocity models and new results of re!ection pro"ling in Finland[END_REF], complemented with low S-wave and Q values in the uppermost crust as observed locally by [START_REF] Pedersen | Depth dependance of q beneath the baltic shield inferred from modeling of short period seismograms[END_REF] along a 200km long profile south of the present study area to obtain similar surface/body wave amplitude ratios. A first conclusion of this study in terms of the local crustal structure is therefore that the low S-velocity and low Q model is widespread over the whole study area.

The conditions that need to be met for a successful and systematic use of body waves for lithospheric studies are still uncertain. The first issue is the minimum amount of data needed to observe the body wave contribution to the GreenÕs function. Theoretically the correlation function converges to the complete GreenÕs function as the square root of the time over which the correlation is evaluated. Such duration dependency is also present in our correlations, when we calculate the amplitude ration of the PmP phases and the remnant fluctuations for different durations of analysis. For the data at hand, good PmP arrivals with an acceptable Signal to Noise Ratio (SNR) are observed after just one month of time averaging. This suggests that travel time measurements can be performed even with limited reording duration.

The second issue is how the noise source distribution and its distance from the array affect the high frequency noise correlations. A previous study south of our study area [START_REF] Pedersen | Influence of the seismic noise characteristics on noise correlations in the Baltic shield[END_REF] points towards the generation of the high-frequency seismic noise along the eastern Atlantic coastline during the winter season. This could potentially have a major impact on our observed noise correlations, and possibly explain the time asymmetry of the observed body waves. For such distant sources, two situations can be hypothesized.

Firstly, the presence of scattered energy from structures outside the study region can strongly contribute to the convergence of the correlation to the GreenÕs function by producing an isotropic, equipartitioned field around the stations. If the stations considered are close enough, and scattering and attenuation weak enough, balistic waves can still be observed with sufficient amplitude to emerge from the fluctuations. In a second situation, if the scattering and attenuation are strong along the path between the two stations considered, the identification of weak ballistic arrivals hidden in the correlation fluctuations will be impossible. These issues are explored in [START_REF] Larose | Fluctuations of correlations and GreenÕs function reconstruction:role of scattering[END_REF] who presented a heuristic model for the SNR in correlations of signals considering specifically the role of scattering in heterogeneous media. Note that the SNR is the ratio between actual Green function and the remnant fluctuation level of the correlation. The SNR is expected to decrease with increasing absorption and with distance. In presence of scattering, [START_REF] Larose | Fluctuations of correlations and GreenÕs function reconstruction:role of scattering[END_REF] also showed that the SNR is behaving following two regimes: SNR is increasing with scattering strength for distances smaller than the transport mean free path l*, while it is decreasing with scattering strength for distance larger than l*. This is in agreement with the fact that regional body waves have so far been detected in cratons (Zhan et al., 2010, this study), characterized by weak attenuation and large mean free path. In Finland, attenuation measurements for S wave in the crust [START_REF] Uski | A new local magnitude scale for the Finnish seismic network[END_REF] suggest that the mean free path is at least on the order of the aperture of the LAPNET network. Further work must be carried out to explore whether body waves can be extracted for all types of crustal structure,.More precisely, in strongly heterogeneous crustal structures, wave scattering could be sufficient to reduce the amplitude of the body waves to a point where they would be hidden in the fluctuations of the correlation functions. The result we report here is encouraging even though more work is required to demonstrate whether our results can be generalized to other geological contexts to open the possibility of the use of noise derived body waves for systematic imaging of the Earth's interior. particle motion analysis for synthetic seismograms calculated for the same distance between the station KIF-LP51, ZVF is the vertical synthetic seismogram generated using a vertical point source at the Earth surface, RVF is the radial synthetic seismogram generated using a vertical point source at the Earth surface. The red line on the PmP particle motion analysis is the linear regression of the ZVF and RVF motion, the coefficient of linearity is 0.99 and the polarization angle is 56¡ to vertical.

Table 1 : Crustal model used to calculate the synthetics seismograms. Vp is the P wave velocity, Vs is the S waves velocity, Qs is the P waves quality factor and Qs is the S waves quality factor.

  , using a vertical point source located at the EarthÕs surface. The vertical and radial component seismograms, calculated at a 100 points at 10km distance interval, correspond to the GreenÕs function G ZZ and G RR which we need to compare to the Z-Z and R-R correlations. The vertical and radial components of the correlations (blue) and synthetic seismograms (black) are shown in Figure 6. All the signals are filtered in the frequency range 0.5 to 1 Hz. The synthetic seismograms show dominant fundamental mode Rayleigh waves on the vertical component, as observed on the Z-Z correlations. The regularly spaced synthetic vertical component seismograms clearly show the single and multiply Moho reflected S waves beyond the critical distance of ~ 110 km of SmS and up to distances of 350km.

Figure 1 :

 1 Figure 1: Map of the geometry of the POLENET/LAPNET array. The black circles

Figure 2 :

 2 Figure 2: Cross correlations of vertical (ZZ) components stacked over 1 year (2008) plotted as

Figure 3 :

 3 Figure 3: Cross correlations of radial (RR) components stacked over 1 year (2008) plotted as

Figure 4 :

 4 Figure 4: Cross-correlations of transverse (TT) components stacked over 1 year (2008) plotted

Figure 5 :

 5 Figure 5: Comparison of the acausal part of (a) ZZ and (c) RR cross correlations plotted as

Figure 6 :

 6 Figure 6: Comparison of the acausal part of the (a) ZZ, (c) RR cross correlations plotted as

Figure 7 :

 7 Figure 7: a) particle motion analysis for the cross-correlation between the station KIF-LP51,

  Figure 1

  Figure 4

ACKNOWLEDGEMENT

We greatly acknowledge support from the QUEST Initial Training network funded within the EU Marie Curie Programme. This study received supported from the ANR BegDy project, the Institut Paul Emil Victor, and European Research Council through the advanced grant ÒWhisperÓ 227507. We thank E. Larose for useful discussions. Synthetic seismograms were calculated using Computer Program in Seismology (Herrmann, R. B., 1996). The POLENET/LAPNET project is a part of the International Polar Year 2007-2009 and a part of the POLENET consortium, and received financial support from The Academy of Finland (grant No. 122762) and University of Oulu, ILP (International Lithosphere Program) task force VIII, grant No. IAA300120709 of the Grant Agency of the Czech Academy of Sciences, and the Russian Federation : Russian Academy of Sciences (programs No 5 and No 9). The Equipment for the temporary deployment was provided by: RESIF Ð SISMOB, FOSFORE, EOST-IPG Strasbourg Equipe sismologie (France), Seismic pool (MOBNET) of the Geophysical Institute of the Czech Academy of Sciences (Czech Republic), Sodankyla Geophysical Observatory (FINLAND), Institute of Geosphere Dynamics of RAS (RUSSIA), Institute of Geophysics ETH ZŸrich (SWITZERLAND), Institute of Geodesy and Geophysics, Vienna University of Technology (AUSTRIA), University of Leeds (UK). The POLENET/LAPNET working group consists of: Elena Kozlovskaya, Teppo JŠmsen, Hanna Silvennoinen, Riitta Hurskainen, Helle Pedersen, Catherine Pequegnat, Ulrich Achauer, Jaroslava Plomerova, Eduard Kissling, Irina Sanina, Reynir Bodvarsson, Igor Aleshin, Ekaterina Bourova, Evald BrŸckl, Tuna Eken Robert Guiguet, Helmut Hausmann, Pekka Heikkinen, Gregory Houseman, Petr Jedlicka, Helge Johnsen, Elena Kremenetskaya, Kari Komminaho, Helena Munzarova, Roland Roberts , Bohuslav Ruzek, Hossein Shomali, Johannes Schweitzer, Artem Shaumyan, Ludek Vecsey, Sergei Volosov. We thank two anonymous reviewers who helped to improve the manuscript.


Matthieu Garcin

Dominique Guégan

EXTREME VALUES OF RANDOM OR CHAOTIC DISCRETIZATION STEPS *

Keywords: AMS 2000 subject classifications: Primary 60E15; secondary 60F99, 62H05 spacings, extreme values, copula, discretization, dynamical systems, invariant measure

By sorting independent random variables and considering the difference between two consecutive order statistics, we get random variables, called steps or spacings, that are neither independent nor identically distributed. We characterize the probability distribution of the maximum value of these steps, in three ways: i/ with an exact formula; ii/ with a simple and finite approximation whose error tends to be controlled; iii/ with asymptotic behavior when the number of random variables drawn (and therefore the number of steps) tends towards infinity.

The whole approach can be applied to chaotic dynamical systems by replacing the distribution of random variables by the invariant measure of the attractor when it is set.

The interest of such results is twofold. In practice, for example in the telecommunications domain, one can find a lower bound for the number of antennas needed in a phone network to cover an area. In theory, our results take place inside the extreme value theory extended to random variables that are neither independent nor identically distributed.

1. Introduction. Ask hundred people to give a number between zero and a thousand, and then sort their answers. What is the maximum step size between two successive values? This apparently simple question, depending of course on the probability distribution governing the choice of the people, has a concise solution, highlighting the probability of the maximum size of a step. The present article attempts to answer this question in a more general context.

We can introduce the problem differently in order to identify its practical interest. We consider a grid consisting of a certain number of independent, identically distributed random variables. This grid can be used for various purposes such as, for instance, the reconstruction of a signal. Indeed, imagine a mixing dynamical system and noting the state of the system at random observation times, we want to reconstruct the entire attractor if it exists. In order to solve this problem, it is useful to link the sample size to the size of the greatest step between two successive values of the ordered variables drawn from the invariant measure of the attractor. This will provide an idea of the approximation made in the reconstruction of the attractor. A practical question may be: How many observations do we need to have a 95% probability (we are particularly interested in high quantiles) of getting a maximum discretization step size below a certain threshold?

In this paper, we present three main results concerning the problem of extreme step. First, we provide an exact formula for the distribution function of the maximum step size, given in Proposition 1. This result being not easy to use in practice, we introduce some approximations: a first approximation as given in Theorem 1 and easy to calculate, and an asymptotic distribution as given in Theorem 2, when the number of random variables tends towards infinity. This discretization grid approach seems to have some similarities with Extreme Value Theory and Monte Carlo or Quasi-Monte Carlo methods, nevertheless it presents some conceptual differences. We list some below.

Monte Carlo or Quasi-Monte Carlo methods: the core of the method involves calculating expectations of a certain function using a discretization grid but it is quite different from the method presented here. Indeed, authors use several kinds of discretization grids for Monte Carlo. Some of them [START_REF] Glynn | Importance sampling for stochastic simulations[END_REF] transform the probability distribution of the signal in order to draw more variables in significant areas and use this importance sampling to the detriment of a small maximum step size in the grid. Other authors [START_REF] Sobol | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF][23] use deterministic low discrepancy sequences of variables in order to reduce the step sizes. In contrast to these works, the discretization grid in our study is simply composed of random variables drawn by a given distribution. Extreme Value Theory: the aim of such a theory is to describe the limit probability distribution of the maximum (scaled and translated) of a set of independent and identically distributed random variables [10][12][22]. Examples of non-independent variables have been proposed to generalize the concept, for example by replacing independence by a strong mixing condition [START_REF] Leadbetter | Limit theorems for strongly mixing stationary random measures[END_REF][19] [18][3]. Similarly, in our study, we seek the probability distribution of the maximum of a set of random variablesthe discretization steps -which are neither independent nor identically distributed. We provide two results, one with finite samples, the other one in an asymptotic setting. Nevertheless, the problem of maximum step size, also known as maximum spacing, is not new. Indeed, some theoretical articles, such as those of Lévy [START_REF] Lévy | Sur la division d'un segment par des points choisis au hasard[END_REF], Deheuvels [START_REF] Deheuvels | Strong limit bounds for maximal uniform spacings[END_REF][5] [START_REF] Deheuvels | On the influence of the extremes of an i.i.d. sequence on the maximal spacings[END_REF] and Devroye [START_REF] Devroye | Laws of the iterated logarithm for order statistics of uniform spacings[END_REF][8] [START_REF] Devroye | The largest exponential spacing[END_REF] address that problem of maximum spacing for some distributions. These authors impose an asymptotic framework. We use a different approach based on the countermonotonicity of the spacings which provides approximations suitable for small samples. Stevens [START_REF] Stevens | Solution to a geometrical problem in probability[END_REF] obtained the exact distribution for a finite sample, but only for a uniform distribution. Our results may be applied to more general probability distributions. Some applications using approximations of the distribution function of a maximum step size exist: for instance, [START_REF] Bettstetter | On the minimum node degree and connectivity of a wireless multihop network[END_REF][2] [START_REF] Santi | A probabilistic analysis for the radio range assignment problem in ad hoc networks[END_REF] consider an ad hoc network modelled by a set of identically distributed random variables, representing mobile antennas. Their goal is to find the optimal transmission range of these antennas, in order for the network to work with a high probability, in the sense that all the nodes would be connected. The maximum difference between two random variables is crucial in this context. These articles make some approximations and assumptions which are different from those in our approach. However, for the uniform distribution, we obtain similar results in the asymptotic case which is their framework. Our main contribution is thus the study of very general distributions in a non-asymptotic case, taking into account the dependence of the steps.

The article is divided into three parts, describing successively the framework (Section 2), the theory (Section 3) and some applications (Section 4).

2. Framework: building the discretization grid. We introduce some notations to define the framework: U 1 , ..., U T are T > 1 independent and identically distributed real random variables; U 1:T ≤ U 2:T ≤ ... ≤ U T :T are the corresponding order statistics; δ X and ∆ X are the probability density function and the cumulative distribution function of a random vector X; in particular, δ U (supposed to be continuous) and ∆ U are associated to the one-dimensional random variables U 1 , ..., U T ; Supp(f ) is the support of the function f ; V 1 , ..., V T are the discretization steps, defined by:

V t = U t:T -U (t-1):T if t > 1 U 1:T -inf(Supp(δ U )) else.
A question may arise about the expression of δ U in practice. We mentioned in the introduction the possibility that these random variables would be drawn according to the invariant measure of the attractor of a dynamical system, if such invariant measure exists. Despite the fact that such invariant measures are often difficult to identify, some attractors possess advantageous properties that allow us to find their invariant measure. The book [START_REF] Guégan | Les Chaos en finance : approche statistique[END_REF] addresses this problem and obtains the invariant measure of some chaos via the Frobenius-Perron operator. We discuss these examples later.

Our goal is to determine the cumulative distribution function of the maximum step size, ∆ max t∈{1,...,T } Vt , and more specifically the probability that the maximum step size is below a certain threshold v > 0:

∆ max t∈{1,...,T } Vt (v) = ∆ V 1 ,...,V T (v, ..., v).
Though we are able to identify δ V 1 ,...,V T , such a cumulative distribution function, ∆ V 1 ,...,V T , is difficult to calculate because of the successive integrations of the density δ U . We propose different ways of approaching this distribution function. The first one reduces the problem in determining a unique integral, and the other one adds a deterministic search for maximum, which is algorithmically much less costly than the calculation of nested integrals.

An additional difficulty concerns the dependence between all the step sizes: indeed, if the step sizes were independent, the joint probability would be equal to the product of the Documents de Travail du Centre d'Economie de la Sorbonne -2012.33 marginal probabilities, and we would easily find a bound to the step size quantiles by changing each of these marginal probabilities by the one of an ideal step corresponding to the lowest density of the invariant measure (where step sizes are bigger). This dependence forces us to use the theory of copulas. However, the expression of the copula for the vector (V 1 , ..., V T ) requires explicit knowledge of the marginal probability of each component V t . Since this is not trivial, we simply build our approximation using a classical result of the theory of copulas based on the Fréchet-Hoeffding theorem.

3. Main results. Before we turn specifically to the useful approximation of the probability distribution of the maximum step size, we present an intermediary result.

3.1. Exact probability distribution of the greatest discretization step size. The first proposition deals with the exact probability distribution function of the maximum of all the step sizes.

Proposition 1. Let U 1 , ..., U T be T > 1 independent and identically distributed real random variables, with a probability density δ U , and V 1 , ..., V T be random variables (the discretization steps) defined by:

(1)

V t = U t:T -U (t-1):T if t > 1 U 1:T -inf(Supp(δ U )) else.
Then, the probability density function of the maximum of the step sizes,

(V 1 , ..., V T ), is, for v ≥ 0: (2) ∆ max t∈{1,...,T } Vt (v) = T ! [0,v] T T k=1 δ U k m=0 v m dv 1 ...dv T ,
where v 0 denotes inf(Supp(δ U )).

Proof. First, we recall that the random variables V 1 , ..., V T are neither independent nor identically distributed.

Let φ : R T -→ R be a Borel function. Let u 1 , ..., u T ∈ Supp(δ U ) and v 1 , ..., v T ≥ 0. We characterize the density δ V 1 ,...,V T of the vector of the steps, defined in equation ( 1), by:

E [φ(V 1 , ..., V T )] = R T φ(v 1 , ..., v T )1 {v 1 ≥0,...,v T ≥0} δ V 1 ,...,V T (v 1 , ..., v T )dv 1 ...dv T = R T φ(v 1 , ..., v T )1 {v 1 ≥0,...,v T ≥0} δ V 1 +v 0 ,...,V T (v 1 + v 0 , ..., v T )dv 1 ...dv T = R T φ(u 1 -v 0 , u 2 -u 1 , ..., u T -u T -1 )1 {v 0 ≤u 1 ≤...≤u T } δ U 1:T ,...,U T :T (u 1 , ..., u T )du 1 ...du T .
The integration is done with the following substitution in order to link order statistics and steps:

       u 1 u 2 . . . u T -1 u T        =          1 0 • • • • • • 0 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 1 • • • • • • 1 1                 v 1 + v 0 v 2 . . . v T -1 v T        . Documents de Travail du Centre d'Economie de la Sorbonne -2012.33
Since the Jacobian of such a substitution is 1, we get:

E [φ(V 1 , ..., V T )] = R T φ(v 1 , ..., v T )1 {v 1 ≥0,...,v T ≥0} δ U 1:T ,...,U T :T 1 m=0 v m , ..., T m=0 v m dv 1 ...dv T ,
which leads to:

(3) δ V 1 ,...,V T (v 1 , ..., v T ) = δ U 1:T ,...,U T :T 1 m=0 v m , ..., T m=0 v m .
We also know [START_REF] Gut | An Intermediate course in probability[END_REF] that the distribution function of the vector (U 1:T , ..., U T :T ) of the order statistics is:

δ U 1:T ,...,U T :T (u 1 , ..., u T ) = T ! T t=1 δ U (u t )1 u 1 ≤...≤u T .
This relationship, used with (3), provides the joint probability density function of the vector of the steps, (V 1 , ..., V T ):

(4) δ V 1 ,...,V T (v 1 , ..., v T ) = T ! T k=1 δ U k m=0 v m .
Then we get:

(5) ∆ max t∈{1,...,T } Vt (v) = ∆ V 1 ,...,V T (v, ..., v) = [0,v] T δ V 1 ,...,V T (v 1 , ..., v T )dv 1 ...dv T .
Finally, the equations ( 4) and ( 5) provide the cumulative distribution function of the maximum of the step sizes, that is to say equation (2).

The Proposition 1 shows that the calculation of the cumulative distribution function of the maximum step size, ∆ max t∈{1,...,T } Vt (v), is difficult to compute because of the successive integrations of the density δ U . Thus, we now provide specific approximations.

3.2.

Estimation of the distribution of the greatest discretization step size. We emphasized the algorithmic difficulty of calculating the exact value of the cumulative distribution function of the maximum step size. Indeed, it amounts to the calculation of T nested integrals. Thanks to the theory of copulas, we propose a lower bound to the probability of the maximum step size. Theorem 1 provides such a lower bound, whose value is obtained by calculating a single integral.

Theorem 1. Let v ≥ 0 and U 1 , ..., U T be T > 1 independent and identically distributed real random variables, with a probability density δ U . We consider the random variables V 1 , ..., V T (the discretization steps that are neither independent nor identically distributed) defined in equation (1): Then, there exists a lower bound ∆ -, defined by:

V t = U t:T -U (t-1):T if t > 1 U 1:T -inf(Supp(δ U )) else.
(6) ∆ - max t∈{1,...,T } Vt (v) ≤ ∆ max t∈{1,...,T } Vt (v),
or, equivalently, for p ∈ [0, 1]:

∆ max t∈{1,...,T } Vt -1 (p) ≤ ∆ - max t∈{1,...,T } Vt -1 (p),
such that:

(7) ∆ - max t∈{1,...,T } Vt (v) = max 0, 1 -T sup(Supp(δ U )) inf(Supp(δ U ))+v δ U (u) [∆ U (u -v) + 1 -∆ U (u)] T -1 du .
In this theorem, we use the inequality of Fréchet-Hoeffding. We assume that the lower bound is very close to the probability that we are interested in, because of the countermonotonicity of the variables studied, especially for high quantiles for which this effect is more obvious and which corresponds to the objective of this article. Indeed, it is well known that the lower Fréchet-Hoeffding bound does not correspond to any existing copula in any dimensions higher than two. Kettler [START_REF] Kettler | Fréchet-Hoeffding lower limit copulas in higher dimensions[END_REF] proposes a realistic lower bound, corresponding to the copula of step sizes when the support of δ U is bounded. That seems to be the best example of countermonotonic variables. However, since the marginal distribution functions in our case are too complicated, we would rather use the Fréchet-Hoeffding inequality. The approach proposed by Kettler about step sizes encourages us to believe that the Fréchet-Hoeffding bound is close to the reality that we are studying.

Proof. Let v ≥ 0. We recall that, according to the Fréchet-Hoeffding theorem, we can write:

(8) max 0, -(T -1) + T t=1 P(V t ≤ v) ≤ P (V 1 ≤ v, ..., V T ≤ v) ≤ min t∈{1,...,T } P(V t ≤ v). Moreover P(V t ≤ v) = R P(V t ≤ v, U t:T = u)du = R P(V t ≤ v|U t:T = u)δ U t:T (u)
du, which, by inverting the finite sum and the integral and using the Lemma 1 postponed at the end of the article, leads to: (9)

T t=1 P(V t ≤ v) = sup(Supp(δ U )) inf(Supp(δ U ))+v T t=1 δ U t:T (u) -T t=1 ∆ U (u-v) ∆ U (u) t-1 δ U t:T (u) du + inf(Supp(δ U ))+v inf(Supp(δ U )) T t=1 δ U t:T (u) -T t=2 ∆ U (u-v) ∆ U (u) t-1 δ U t:T (u) du.
We know the expression of the probability density function of the order statistic, given in the proof of the Lemma 1 in equation (29):

δ U t:T (u) = T ! (T -t)!(t -1)! [∆ U (u)] t-1 [1 -∆ U (u)] T -t δ U (u).
Documents de Travail du Centre d'Economie de la Sorbonne -2012.33

Thus:

(10)

T t=1 δ U t:T (u) = δ U (u)T T -1 s=0 (T -1)! s!(T -1-s)! [∆ U (u)] s [1 -∆ U (u)] T -1-s = δ U (u)T [∆ U (u) + 1 -∆ U (u)] T -1 = δ U (u)T,
and, for any ρ > 0:

(11)

T t=1 ρ t-1 δ U t:T (u) = δ U (u)T [ρ∆ U (u) + 1 -∆ U (u)] T -1 .
On the other hand, for u < v + inf(Supp(δ U )), ∆ U (u -v) = 0 and therefore: [START_REF] Gnedenko | Sur la distribution limite du terme maximum d'une série aléatoire[END_REF] T t=2

∆ U (u-v) ∆ U (u) t-1 δ U t:T (u) = 0.
Thus, using equations ( 9) to [START_REF] Gnedenko | Sur la distribution limite du terme maximum d'une série aléatoire[END_REF], the sum of the cumulative distributions of the step sizes becomes:

(13)

T t=1 P(V t ≤ v) = T -T sup(Supp(δ U )) inf(Supp(δ U ))+v δ U (u) [∆ U (u -v) + 1 -∆ U (u)] T -1 du.
Finally, we use the lower Fréchet-Hoeffding bound [START_REF] Devroye | A log log law for maximal uniform spacings[END_REF] with the equation ( 13) and we get ( 6) and [START_REF] Devroye | Laws of the iterated logarithm for order statistics of uniform spacings[END_REF].

Even though the exact probability is much closer to the lower Fréchet-Hoeffding bound than to the upper one, we can provide an upper bound to the probability of the maximum step size, for comparison, thanks to the upper Fréchet-Hoeffding bound:

∆ max t∈{1,...,T } Vt (v) ≤ ∆ + max t∈{1,...,T } Vt (v), with: (14) ∆ + max t∈{1,...,T } Vt (v) = min t∈{1,...,T } T ! (T -t)!(t -1)! R ∆ U (u) t-1 -∆ U (u -v) t-1 (1 -∆ U (u)) T -t δ U (u)du .
Theorem 1 provides an interesting expression, but in certain cases it will be difficult to calculate it analytically. A good solution would be to discretize the integral which appears in the theorem, in a way that preserves the fact that our estimator is a lower bound. Thus instead of an integral, we can consider a sum together with a deterministic search for maximum. This is the subject of the following corollary:

Corollary 1. Let v ≥ 0 and U 1 , ..., U T be T > 1 independent and identically distributed real random variables, with a probability density δ U . We consider the random variables V 1 , ..., V T defined in equation [START_REF] Bettstetter | On the minimum node degree and connectivity of a wireless multihop network[END_REF]. Assume that the support of δ U is bounded. Let Θ ∈ N and (S 1 , ..., S Θ ) be a set of intervals forming a partition of [inf(Supp(δ U )) + v, sup(Supp(δ U ))].

We define ∆ -,(S 1 ,...,S Θ ) max where:

ρ S θ (v) = max u∈S θ δ U (u) [∆ U (u -v) + 1 -∆ U (u)] T -1 .
Then:

(16) ∆ -,(S 1 ,...,S Θ ) max t∈{1,...,T } Vt (v) ≤ ∆ - max t∈{1,...,T } Vt (v)
and:

lim max θ∈[1,Θ] [sup(S θ )-inf(S θ )]→0 ∆ -,(S 1 ,...,S Θ ) max t∈{1,...,T } Vt (v) = ∆ - max t∈{1,...,T } Vt (v).
Proof. All the results are a direct consequence of the discretization of the integral appearing in Theorem 1, more precisely in equation ( 7):

∆ - max t∈{1,...,T } Vt (v) = max 0, 1 -T sup(Supp(δ U )) inf(Supp(δ U ))+v δ U (u) [∆ U (u -v) + 1 -∆ U (u)] T -1 du = max 0, 1 -T Θ θ=1 S θ δ U (u) [∆ U (u -v) + 1 -∆ U (u)] T -1 du ≥ max 0, 1 -T Θ θ=1 [sup(S θ ) -inf(S θ )] max u∈S θ δ U (u) [∆ U (u -v) + 1 -∆ U (u)] T -1 , because, for u ∈ R, δ U (u) [∆ U (u -v) + 1 -∆ U (u)] T -1 ≥ 0.
Then, with the definition of ∆ -,(S 1 ,...,S Θ ) max t∈{1,...,T }

Vt in ( 15), we get (16).

3.3. Towards a limit distribution. One objective of the extreme value theory is to find the limit distribution function of the maximum of a set of random variables. In the theory we developed above, we did not consider limit distributions, but we tried to describe precisely what happens for a given number of steps. Nevertheless, we can definitely present a limit theorem describing the shape of the lower bound of the probability when drawing an infinite number of variables.

Theorem 2. Let v ≥ 0 and U 1 , ..., U T be T > 1 independent and identically distributed real random variables, with a probability density δ U supposed to be differentiable. We consider the random variables V 1 , ..., V T (the discretization steps that are neither independent nor identically distributed) defined in equation (1):

V t = U t:T -U (t-1):T if t > 1 U 1:T -inf(Supp(δ U )) else.
When v is greater than a certain threshold v(T ) such that lim

T →∞ T g (v(T )) = 1, where g : v -→ R δ U (u)e -δ U (u)v du, then: (17) ∆ - max t∈{1,...,T } Vt v T T →∞ ∼ max(0; 1 -T g(v)).
Documents de Travail du Centre d'Economie de la Sorbonne -2012.33

In the equation ( 17), we used the notation a(T )

T →∞ ∼ b(T ), for functions a and b. This means that, when lim

T →∞ b(T ) = 0, then lim T →∞ a(T ) b(T ) = 1.
Proof. According to Theorem 1: (18)

∆ - max t∈{1,...,T } Vt v T = max 0, 1 -T sup(Supp(δ U )) inf(Supp(δ U ))+v/T δ U (u) [∆ U (u -v/T ) + 1 -∆ U (u)] T -1 du .
Moreover, for u ∈ R, we get the Taylor expansions:

(19) [∆ U (u -v/T ) + 1 -∆ U (u)] T →∞ = 1 -δ U (u) v T + o v T , because δ U is differentiable. Thus: (20) 1 -δ U (u) v T T -1 T →∞ = exp (-vδ U (u)) + o v T .
Putting ( 19) and ( 20) together leads to:

(21) [∆ U (u -v/T ) + 1 -∆ U (u)] T -1 T →∞ = exp (-vδ U (u)) + o v T .
We can now integrate the Taylor expansion [START_REF] Lévy | Sur la division d'un segment par des points choisis au hasard[END_REF] and use it in ( 18): ( 22)

∆ - max t∈{1,...,T } Vt v T T →∞ = max 0, 1 -T sup(Supp(δ U )) inf(Supp(δ U )) δ U (u) exp (-vδ U (u)) du + o (1) .
When v is greater than a certain threshold v(T ) such that lim

T →∞ T g (v(T )) = 1, which
implies that the right part of the equation ( 22) is greater than zero, we get (17).

Examples for diverse invariant measures.

We can apply the previous results to any type of distribution known analytically. We give examples for the attractors characterized by maps defined on [0, 1]. See Guégan [START_REF] Guégan | Les Chaos en finance : approche statistique[END_REF] and Lasota and Mackey [START_REF] Lasota | Probabilistic properties of deterministic systems[END_REF][17] for interesting developments on the subject. We restrict to the logistic attractor and the 2adic attractor and we provide the distribution of their maximum step size. We also study the case of a density with unbounded support: the exponential distribution. Finally, we discuss the building of the discretization grid.

In addition, all our theoretical work can be applied without any modification to more general probability density functions and thus to the search for the maximum step size for any random variable.

4.1. Logistic attractor. The logistic attractor of parameter α is defined by the recurrence relation:

U n+1 = αU n (1 -U n ).
Documents de Travail du Centre d'Economie de la Sorbonne -2012.33

In a very general case, no theoretical result seems to reveal an invariant measure for the logistic attractor. However, if its parameter is 4, which is then literally a chaotic system, we know the existence and shape of the invariant measure, µ 4 :

µ 4 : x ∈ [0, 1] -→ 1 π x(1 -x) ,
in which we recognize a beta density function of parameters 1/2 and 1/2 -see Ulam and von Neumann [START_REF] Ulam | On combination of stochastic and deterministic processes[END_REF]. We represent the density in Figure 1. We note M 4 as the associated cumulative distribution function: We apply Corollary 1 to this attractor and represent the result in Figure 2. We note that the upper bound of quantiles (the corresponding curve is below the empirical distribution) approaches reality for high probabilities, especially when increasing the number Θ of intervals in the partition (S 1 , ..., S Θ ). That proximity between the empirical quantiles and our approximation corroborates the intuition that the step sizes are countermonotonic and thus that the lower inequality of Fréchet-Hoeffding applied to the distribution of the step sizes V t is close to reality. Indeed, when a step V t is much larger than others, there is less room for other steps, therefore other steps will potentially be smaller.

M 4 : x ∈ [0, 1] -→ x 0 1 π y(1 -y) dy = 1 2 + 1 π arcsin(2x -1).
For probabilities higher than 85%, the empirical quantile and its upper bound, calculated for a partition of the support in 100 intervals, are nearly equal, as we can see in Figure 3. Its invariant measure, when r = 2, is the constant function [START_REF] Guégan | Les Chaos en finance : approche statistique[END_REF]:

µ 2 : x ∈ [0, 1] -→ 1.
In fact, that is also the probability density function of a uniform random variable. We note M 2 as the corresponding cumulative distribution function, which is the identity function.

In this simple case, we can apply directly Theorem 1 without using a partition of the support. We get the following expression for the lower bound of the distribution function:

(23) ∆ - max t∈{1,...,T } Vt (v) = max 0, 1 -T [v,1] [1 -v] T -1 du = max 0, 1 -T (1 -v) T .
The limit distribution function is obtained thanks to Theorem 2:

∆ - max t∈{1,...,T } Vt v T T →∞ ∼ max(0; 1 -T e -v ).
For this attractor, the upper bound is obtained from equation ( 14) by successive integrations by parts or by recognizing Beta functions:

(24) ∆ + max t∈{1,...,T } Vt (v) = min t∈{1,...,T } T ! (T -t)!(t-1)! R M 2 (u) t-1 -M 2 (u -v) t-1 (1 -M 2 (u)) T -t µ 2 (u)du = min t∈{1,...,T } T ! (T -t)!(t-1)! [0,1] u t-1 (1 -u) T -t du -[v,1] (u -v) t-1 (1 -u) T -t du = min t∈{1,...,T } 1 -(1 -v) T = 1 -(1 -v) T .
We represent both the upper and lower bounds in Figure 4. We note that the upper bound of quantiles, corresponding to the lower Fréchet-Hoeffding inequality, converges on the empirical quantiles for high probabilities. On the other hand, the lower bound of quantiles, corresponding to the upper Fréchet-Hoeffding inequality, is very far from the empirical quantiles. These two facts are consistent with our intuition that the step sizes are countermonotonic.

For probabilities higher than 85%, the empirical quantile and its upper bound are nearly equal, as we can see in Figure 5.

Exponential distribution.

In order to give an example of unbounded support of probability density, we study the case of exponential random variables of parameter λ, defined by the probability density function:

µ λ : x ≥ 0 -→ λe -λx ,
and by the cumulative distribution function:

M λ : x ≥ 0 -→ 1 -e -λx .
Documents de Travail du Centre d'Economie de la Sorbonne -2012.33 reach a target of proximity for all the nodes. As far as our approximation is an overestimation of the quantile, then the minimum number of observations claimed by the same approximation is also slightly higher than in reality. Thus, our approximation allows us to build a discretization grid which is at least as accurate as we want.

For example, for a draw of random variables in the invariant measure of the 2-adic attractor, we recall the estimation obtained in equation ( 23):

∆ - max t∈{1,...,T } Vt (v) = max 0, 1 -T (1 -v) T .
For a number of observations T and a given probability p, this leads to an upper bound of the maximum step size:

(28) v ≤ 1 - 1 -p T 1 T
.

We present this result in Figure 7, for various values of p.

As an illustration, with a 95% probability, a maximum step size of 0.02 is guaranteed for T > 450. If the probability is 99.9%, then the condition becomes T > 663. Each curve is for one probability, from the top right to the bottom left, 99.9%, 99%, 95% et 90%.

Ask hundred people to give a number between zero and a thousand, and then sort their answers. What is the maximum step size between two successive values?

The answer can be read, scaled, in Figure 7 if the distribution of such numbers is uniform: the maximum step size is less than 67 with a probability greater than 90%, less than 74 with a probability greater than 95%, less than 88 with a probability greater than 99% and less than 109 with a probability greater than 99.9%.

We have then obtained three main innovative theoretical results:

in the Proposition 1, we give an exact formula for the distribution function of the maximum step size; in Theorem 1, we give an approximation based on the countermonotonicity of the step sizes and which is easy to calculate; in Theorem 2, we give a limit distribution, when the number of random drafts tends towards infinity.

The contribution is interesting for the extreme value theory, since we consider the maximum of non-independent and non-identically distributed random variables. We also proposed numerical applications mainly for invariant measures of chaos. More general probability distributions can be used too. We presented a practical application to quantify the accuracy of a random or chaotic discretization grid. In further work, we will invoke the results of this paper to estimate the maximum error made when reconstructing attractors of chaos by the technique of wavelets, applied to some observations forming a discretization grid.

Documents de Travail du Centre d'Economie de la Sorbonne -2012.33

  Documents de Travail du Centre d'Economie de la Sorbonne -2012.33

  t∈{1,...,T } Vt by: (15) ∆ -,(S 1 ,...,S Θ ) max t∈{1,...,T } Vt (v) = max 0, 1 -T Θ θ=1 [sup(S θ ) -inf(S θ )] ρ S θ (v) , Documents de Travail du Centre d'Economie de la Sorbonne -2012.33

Figure 1 .

 1 Figure 1. Probability density function µ4: invariant measure of the logistic attractor of parameter 4.

4. 2 . 2 -

 22 adic attractor. The r-adic attractor is defined by the recurrence relation:U n+1 = rU n mod 1.Documents de Travail du Centre d'Economie de la Sorbonne -2012.33

Figure 2 .

 2 Figure 2. Simulated cumulative distribution function (in black) of the maximum step size for the logistic attractor of parameter 4, with T = 100. The simulation has been carried out with 10, 000 sets of T variables drawn in the invariant measure of the attractor. In red, we represent the lower bound of the cumulative distribution function (the upper bound of quantiles), specified in Corollary 1, for partitions of 10 (solid line) or 100 (dotted line) intervals of equal size.

Figure 3 .

 3 Figure 3. Empirical and approximated quantile of the maximum step size for the logistic attractor of parameter 4, with T = 100. The simulation has been carried out with 10, 000 sets of T variables drawn in the invariant measure of the attractor. The approximated quantile is calculated with Corollary 1 and a partition of 100 intervals of equal size.

Figure 4 .

 4 Figure 4. Simulated cumulative distribution function (in black) of the maximum step size for the 2-adic attractor, with T = 100. The simulation has been carried out with 10, 000 sets of T variables drawn in the invariant measure of the attractor. In red, we represent the lower bound of the cumulative distribution function (the upper bound of quantiles, equation (23)). In green, we represent the upper bound of the cumulative distribution function (the lower bound of quantiles, equation (24)).

Figure 5 .

 5 Figure 5. Empirical and approximated quantile of the maximum step size for the 2-adic attractor, with T = 100. The simulation has been carried out with 10, 000 sets of T variables drawn in the invariant measure of the attractor.

Figure 6 .

 6 Figure 6. Simulated cumulative distribution function (in black) of the maximum step size for exponential random variables, with intensity λ = 0.1 and T = 100. The simulation has been carried out with 10, 000 sets of T variables drawn in the exponential probability distribution. In red, we represent the lower bound of the cumulative distribution function (the upper bound of quantiles, equation (26)). Matching the black line, we represent the approximation by the truncated product of the analytic asymptotic cumulative distribution function, equation (27)), which coincides with the simulated distribution.

5 .

 5 Conclusion. We introduced the problem of extreme values of random or chaotic discretization steps with the facetious question: Documents de Travail du Centre d'Economie de la Sorbonne -2012.33

Figure 7 .

 7 Figure 7. Upper bounds of the maximum step size (for the 2-adic attractor) as a function of the number of observations (T ), as in equation (28). Each curve is for one probability, from the top right to the bottom left, 99.9%, 99%, 95% et 90%.

We can apply directly Theorem 1, without using a partition of the support. We obtain the following expression for the lower bound of the distribution function: Thanks to the substitution x = e -λu in [START_REF] Sobol | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF], we get:

∆ - max t∈{1,...,T } Vt (v) = max 0, 1 + T [e -λv ,0) 1 + x 1 -e λv T -1 dx = max 0, 1 -1-e -λvT 1-e -λv e -λv .

In [START_REF] Devroye | The largest exponential spacing[END_REF] we have an exact asymptotic formula for the cumulative distribution function of the maximum step of an infinite set of exponential random variables which is:

We note that ∆ - max t∈{1,...,T }

Vt may be written as:

which is the first-order expansion of:

T k=1

1 -e -λvk when v is big. We find empirically that this product, inspired by the asymptotic formula of [START_REF] Devroye | The largest exponential spacing[END_REF], is an excellent approximation of the true distribution function. For this reason, we can assert, for the example of exponential random variables, that high quantiles of discretization steps are well approximated by our formula. We can check that in Figure 6.

The upper bound of quantiles has the advantage of being very fast to compute, and although it is an approximation, we know where it stands in relation to reality: we chose to overestimate the quantiles. Moreover, for probabilities higher than 90%, the empirical quantile matches its upper bound.

4.4.

Back to the building of the discretization grid. One of the motivations of this article was: How many observations do we need to have a 95% probability of getting a maximum discretization step size below a certain threshold?

As we are able to provide an upper bound of the quantiles (even more accurate for large probabilities), we can get a minimum number of observations (T ) that allows the grid to Documents de Travail du Centre d'Economie de la Sorbonne -2012.33

APPENDIX A: LEMMA

In Lemma 1, used in the demonstration of Theorem 1, we introduce the conditional probability density of a step, given the position on its right.

Lemma 1. Let U 1 , ..., U T be T > 1 independent and identically distributed real random variables, with a probability density δ U . We consider the random variables V 1 , ..., V T (the discretization steps that are neither independent nor identically distributed) defined in equation [START_REF] Bettstetter | On the minimum node degree and connectivity of a wireless multihop network[END_REF]. Let t ∈ {1, ..., T }, v ≥ 0 and u ∈ R. The probability density of V t conditional to U t:T is:

[∆ U (u)] t-1 and the corresponding cumulative distribution function, defined by the conditional probability:

can be expressed by:

.

Moreover, when t = 1, then V t = U t:T -inf(Supp(δ U )) and the previous result can also be written:

Proof. We prove the lemma in four steps:

1. We first need to have the probability density function of any order statistic. This is a well-known result:

(29)

For the proof, see [START_REF] Gut | An Intermediate course in probability[END_REF]. 2. We also need to obtain the joint density function of two following order statistics, for t ≥ 2. We first write the joint density function of all the order statistics:

For the proof, see [START_REF] Gut | An Intermediate course in probability[END_REF]. Then we integrate this equation: (30)

3. We now want to obtain the joint density function of a step and the corresponding order statistic on the right of the step. Let t ≥ 2. Let φ : R 2 -→ R be a Borel function. We characterize the density δ U t:T ,Vt by: We implement the following substitution:

Since the Jacobian of such a substitution is 1, we get:

which leads to:

This relationship, used with (30), provides the joint probability of (U t:T , V t ): (31)

4. According to Bayes' theorem, the conditional density function δ Vt|U t:T is:

With equations ( 29) and (31), we conclude for t ≥ 2 that:

The cumulative distribution function is then obtained by integration:

.

The case t = 1 is trivial anyway. Indeed, U 1:T = V 1 + inf(Supp(δ U )), and we get:


Vincent Luboz
email: vluboz@imag.fr

Emmanuel Promayon
email: epromayon@imag.fr

Grégory Chagnon
email: gregory.chagnon@grenoble-inp.fr

Thierry Alonso
email: thierry.alonso@ujf-grenoble.fr

Denis Favier
email: denis.favier@grenoble-inp.fr

Christine Barthod
email: christine.barthod@univ-savoie.fr

Yohan Payan
email: ypayan@imag.fr

Validation of a Light Aspiration device for in vivo Soft TIssue Characterization (LASTIC)

Keywords: soft tissue characterization, aspiration, mechanical behavior, inverse analysis, planning and assistance

In-vivo characterization of soft tissue is a key step towards accurate biomechanical simulation enabling planning and intra-operative assisted surgery. This chapter presents the new version of LASTIC, a device measuring soft tissue deformations using a negative pressure. Its capabilities are compared with standard tensile tests on five samples with different elastic properties, i.e Young modulus, from 10 kPa to 1 MPa, in order to estimate its accuracy and define the functional measurement range. Results show that LASTIC overestimates Young modulus by an average of 24 % compared to the tensile devices. This error, although rather large, allows a first estimation of the elastic modulus of different materials, especially living tissues, even during surgery. Directions for improvements are given that will allow for better patient-specific biomechanical simulations of soft tissues.

Introduction

Living soft tissues are known to exhibit non-linear, inhomogeneous, anisotropic, patient-specific, time and rate-dependent behaviors. Taking into account such specificities is highly challenging, especially for intra-operative computer aided devices that require an on-line estimation of the constitutive behavior of the tissues.

This chapter introduces the Light Aspiration device for in vivo Soft TIssue Characterization (LASTIC) which principles were proposed by our group and used, in its first version, to evaluate the constitutive behavior of skin [START_REF] Schiavone | In-vivo measurement of the human soft tissues constitutive laws[END_REF], tongue [START_REF] Schiavone | A light sterilizable pipette device for the in vivo estimation of human soft tissues constitutive laws[END_REF] and brain tissues [START_REF] Schiavone | In Vivo Measurement of Human Brain Elasticity Using a Light Aspiration Device[END_REF]. In its first version the device was based on the pipette aspiration principle which consists in measuring the tissues deformations induced by a negative pressure. The surgeon held the instrument and established contact with the tissues surface while the device measured the negative pressures and displacement responses. This implied that (1) the device had to undergo a full sterilization and (2) the data processing to be sufficiently fast to provide an interactive estimation of the tissues constitutive equation. These two constraints have encouraged us to propose a new version of LASTIC [START_REF] Schiavone | LASTIC: A Light Aspiration Device for in vivo Soft TIssue Characterization[END_REF]. This chapter aims at validating this new version by comparing, for synthetic materials, the constitutive parameters it estimates with the parameters provided by usual tensile test machines. This validation is indeed a prerequisite before any intra-operative use of the new version of LASTIC.

Next section introduces the samples and the devices that were used to estimate the materials' mechanical parameters. The remaining provides the corresponding results (section 3) before a discussion (section 4) and a conclusion (section 5).

Material and methods

Elastic samples

Although human soft tissues are considered to be nonlinear viscoelastic, most of the work done to simulate them analyzes their behavior on the linear elastic side. To provide a better coverage of human soft tissues' elasticity and behavior, experiments were conducted on five different samples. Their elastic moduli ranged from few kilopascals to several hundred kilopascals (kPa thereafter). According to the literature, several materials correspond to this elasticity. We therefore chose to use the following ones:

• RTV#1: a RTV-EC00 silicone (artificina.com), made from a mix of 50 % of base and 50 % of catalyst (also called curing agent). It has a linear behavior up to 25 % of engineering strain (see Figure 6). • RTV#2: the same RTV-EC00 silicone, but made from a mix of 40 % of base and 60 % of catalyst, to create a softer silicone than RTV#1. It has a linear behavior up to 25 % of engineering strain (see Figure 6).

• RTV#3: a RTV 141 silicone (artificina.com), consisting of 90 % of base and 10 % of catalyst. It has a linear behavior up to 120 % of engineering strain, where rupture of the material occurs. • Ecoflex: an Ecoflex® 00-30 silicone (Smooth-on.com) constituted by two bases mixed in equal proportion. It has a linear behavior up to 15 % of engineering strain (see Figure 6). • Candle gel: a gel (glorex.com) that is fluid at temperatures over 95 °C. It is fragile at room temperature and tends to tear. It has a linear behavior up to 10 % of engineering strain (see Figure 6). For each material, two samples were created by two molds in order to perform two different tests: one to be used by LASTIC and one to be used for the tensile tests, see subsections 2.2 and 2.3 for details on the mold shapes. This process ensured that the same material mix is poured in the two molds so that the elasticity measurements are done on the same material in each testing machine. Each test is done at least three times on each sample to take the variability of the measurements of the elastic modulus into account, except for the RTV#3 which was extensively characterized in previous studies [START_REF] Meunier | Mechanical experimental characterisation and numerical modelling of an unfilled silicone rubber[END_REF].

RTV#1, RTV#2 and RTV#3 were all created following the same process: first the carefully weighted amount of base and catalyst were mixed in a container according to the proportion described previously, then the mix was exposed to vacuum for several minutes (about five minutes for RTV#1 and RTV#2, and around 30 minutes for RTV#3) in order to evacuate any air bubbles. The mixes were then poured in the two different molds. The silicones RTV#1 and RTV#2 cured in two hours at room temperature while the RTV#3 cured in 24 hours, including three initial hours at 150° C.

The Ecoflex silicone is created from its two bases by mixing them in a container and allowing it to cure for four hours at room temperature. The mix is also submitted to a vacuum chamber. Five minutes seemed enough to evacuate any air bubbles. Its properties are stationary 24 hours later.

The candle gel is heated to 95° C in a boiling water bath. Once liquefied, it is poured in the two molds and becomes an elastic solid in few minutes. Because of its fragility, precautions must be taken when extracting it out of the mold.

LASTIC

Based on a prototype aspiration device quantitatively evaluated during surgery on the brain [START_REF] Schiavone | In Vivo Measurement of Human Brain Elasticity Using a Light Aspiration Device[END_REF]) and during more classical experiments on other organs such as tongue, cheeks and forearm skin, we designed a more elaborate device, called LASTIC (Light Aspiration device for in vivo Soft TIssue Characterization) [START_REF] Schiavone | LASTIC: A Light Aspiration Device for in vivo Soft TIssue Characterization[END_REF]. LASTIC is intended to be used in operating rooms and therefore to meet the very rigorous sterilization and handling process imposed during surgery. The basic design is not very far from the one proposed by [START_REF] Vuskovic | Device for in-vivo measurement of mechanical properties of internal human soft tissues[END_REF], except that it is more compact, see Figure 1.

LASTIC is built in a very compact metallic cylinder of 33 mm in height and 34 mm in diameter. This cylindrical case is divided in two compartments. The lower compartment is an airtight cylindrical chamber, open at the bottom by a circular aperture and closed at the top by a glass window. The upper compartment holds the electronic part consisting of a miniature digital camera and a LED used as a light source. The camera is a 9.5 mm × 9.5 mm 2 megapixels mobile phone camera sensor (model VS6750, STMicroelectronics, Geneva, Switzerland) with a resolution of 1600 × 1200 pixels (UXGA).

The complete setup is composed of a programmable syringe pump (Aladdin AL-1000, World Precision Instruments, Inc., FL, USA), a USB digital manometer (model 8215, AZ Instrument Corp, Taiwan) and LASTIC, all linked by Luer-Lok connected flexible medical tube, see Figure 2.

Figure 1 -Cross section of LASTIC: the upper part contains the camera while the lower part is the aspiration chamber with the mirror (left). View of the two chambers when separated (right).

While the lower compartment is in contact with the tissue, a negative pressure P is generated by a software controlled syringe pump programmed to withdraw. The soft-tissue is consequently aspired inside the chamber, where the deformation, measured as the aspiration height, is imaged by the digital camera, by means of a 45° inclined mirror (see Figure 3). The pressure from the manometer is acquired synchronously. Once the micro-leaks, occurring at the interface between the decompressing chamber and LASTIC, are sealed off by a sufficient minimal negative pressure and LASTIC weight, the measurements can be carried out.

The height of the aspired tissue is then segmented on the image and a basic camera calibration is performed to determine the pixel size. The experiment is performed for increasing level of negative pressure. A displacement height is therefore obtained for each applied negative pressure.

We only considered the apex displacements higher than 0.5 mm, which corresponds to the position of the bottom of the mirror relatively to the initial tissue surface. If this level is not reached, there is a risk of measuring a point lower than the deformation dome but appearing higher on the mirror, due to an incorrect inclination of the setup. For the present experiments, the mold produces a material sample of cylindrical shape of approximately 20 mm of height, 60 mm of diameter.

A Finite Element Analysis of the aspiration experiment using a Neo-Hookean constitutive law [START_REF] Treloar | The elasticity of a network of long chains molecules I and II[END_REF]) is used in order to build a library of displacement heights. The Neo-Hookean equation is written as:

W= C 1 (I 1 -3) (1) where I 1 is the first invariant of the left Cauchy-Green deformation tensor, C 1 is a material parameter, and W is the strain energy.

The library of displacements is basically a double entry table depending on C 1 and the applied negative pressure. From this library, a least-square minimization method is used to find the value C 1 in the library that best matches the measurements. This is performed in interactive time (less than 1 s).

Tensile devices

Because it is a well-established characterization technique, tensile test measurements are used as a reference to compare with LASTIC measurements and estimate the resulting precision on the material parameter. Two extension devices were used for the tensile measurements in two different laboratories.

The first tensile device is an Instron® machine (Instron, Norwood, MA, USA), used in the SYMME Laboratory, see right panel of Figure 4. It applies a given displacement to the material and records the force needed to reach it. The material sample had a rectangular shape of 60 mm of length, 15 mm of width and 7 mm of thickness. It was fixed by tensile jaws on each side of its length.

The second extension device is an Eplexor 500N testing device (Gabo Qualiter Testanlagen GmbH, Ahlden/Aller, Germany), used in the 3SR laboratory with a 25 N sensor. It works similarly to the Instron machine: applying a given displacement and recording the relative force. The material sample was also fixed between two tensile jaws, at each extremity, see left panel of Figure 4. It had a rectangular shape of 25 mm of length, 7 mm of width and 2 mm of thickness. To evaluate the elastic modulus from the measurements of both devices, the tensile test stress/strain values are studied by fitting a Neo-Hookean constitutive equation:

π 11 = 2C 1 (λ -1/λ 2 )
(2) where π 11 is the First Piola-Kirchhoff stress, C 1 is a material parameter, and λ is the principal stretch. The strain energy density function is the same as Equation (1).

During infinitesimal deformations, one can assume that the deformed and initial configurations coincide. The Cauchy stress can consequently be written as: σ11 = π 11 λ (3) As the elongation λ is equal to 1+ ε where ε is the engineering strain (measured by the testing device),

σ 11 = 2C 1 (λ 2 -1/λ)=2C 1 ((1+ε) 2 -1/ (1+ ε))=6C 1 ε +o(ε) (3)
We deduce from this equation the link between the C 1 parameter and the Young modulus E = 6C 1 .

Results

The results of LASTIC measurements on RTV#1, RTV#2, Ecoflex and candle gel are presented in Figure 5 a), b), c) and d) respectively.

For the RTV#1 sample, the minimal apex displacement, or aspiration height, was reached with a negative pressure of 50 mbar. In total, a volume of 10 ml of air was aspirated with LASTIC, resulting in a negative pressure of nearly 200 mbar and an aspiration height of slightly more than 1.3 mm. A Young modulus of 100 kPa was estimated for this material.

For the RTV#2 sample, a negative pressure of 30 mbar was needed to reach enough apex displacement. In total, a volume of 5ml was pulled with LASTIC, resulting in a negative pressure of about 95 mbar and an aspiration height around 2 mm. A Young modulus of 32 kPa was estimated for this material.

For the RTV#3 sample, the minimum negative pressure could not be reached even though a volume of 60 ml was aspired with the pump, resulting in a negative pressure of about 745 mbar and an apex displacement of approximately 0.4 mm. A higher negative pressure could not be obtained. Because of the stiffness of this material and the high negative pressure needed to deform it, the aspiration tests were unsuccessful. No Young modulus was therefore estimated for this material.

For the Ecoflex sample, a negative pressure of 60 mbar was needed to reach the minimal aspiration height. In total, a volume of 5 ml was aspirated, resulting in a negative pressure of about 140 mbar and an apex displacement close to 1.5 mm. A Young modulus of 67.5 kPa was estimated for this material.

For the candle gel sample, a minimum negative pressure of about 5 mbar was needed. In total, a volume of 5 ml of air was aspirated, resulting in a negative pressure of nearly 50 mbar and an aspiration height of slightly less than 2.2 mm. A Young modulus of 14.5 kPa was estimated for this material. Figure 6 shows the measurements made with the two tensile devices. RTV#1, RTV#2 and RTV#3 were tested with the Eplexor 500N Gabo device while the Ecoflex and the candle gel were tested with the Instron machine. Each sample was placed between the tensile jaws before being deformed in extension. A Neo Hookean curve was then fitted on the measurements to estimate the best Young modulus as explained previously. The results of the tensile device measurements on RTV#1, RTV#2, RTV3, Ecoflex and candle gel are presented in Figure 6 a), b), c), d) and e) respectively.

The RTV#1 sample was pulled at two different strain rates to estimate the influence of viscosity on this silicone. The "slow" tensile test was performed at approximately 3 mm.min -1 while the "fast" tensile test was performed at 40 mm.min -1 . No significant difference was recorded between the two measurements: this silicone had no strain rate dependence in this range of strain rates. A Young modulus of 90 kPa was estimated for this sample, using the Neo Hookean fitting. At the first stage of these tensile tests, some high stresses were locally recorded (visible as two spikes on Figure 6 a), they were probably due to some initial sliding of the sample in the tensile jaws.

The RTV#2 measurements were similarly conducted. Again, the measurements at different strain rates showed no viscosity influence for this material. The Neo Hookean law fitting estimated the Young modulus to be 25 kPa.

The measurements on the RTV#3 were highly reproducible and only one test is presented in Figure 6 c). A Young modulus of 950 kPa was estimated with the Neo Hookean law.

The measurements for Ecoflex estimated by the Neo Hookean constitutive law provided a Young modulus of 55 kPa. Finally, the tensile tests on the candle gel gave a Young modulus estimation of 10.5 kPa.

Figure 6 -Tensile measurements and Neo Hookean fitting for (a) RTV#1, (b) RTV#2, (c) RTV#3, (d) Ecoflex, and (e) candle gel.

Figure 7 summarizes the Young modulus estimated by the tensile test devices and by LASTIC. For the candle gel, a Young modulus of 10.5 kPa and 14 kPa were found for the tensile devices and LASTIC respectively, resulting in 33 % difference. Similarly, for the other materials, differences of 28 %, % and 11 % were found for the RTV#2, Ecoflex, and RTV#1 respectively. It leads to an average 23.8 % ± 9.5 % SD overestimation of the Young modulus as given by LASTIC compared to the tensile tests.

Discussion

The average discrepancy of 23.8 % in the value of the Young modulus estimated between LASTIC and the tensile devices is quite high. However, this overestimation seems to be highly linearly dependent of the value of E. If this relatively constant inaccuracy was confirmed by further study, this could be taken into account during soft tissue characterization of living tissues. Nevertheless, for a precise estimation of the Young modulus, we should improve LASTIC by considering at least three main sources of error. Two sources of error are linked with the measurement devices, namely the manometer precision and the camera calibration, and the last one is due to the inverse problem solving.

The first area of possible improvement concerns the manometer. The digital manometer used in our set up has a precision of ± 6 mbar in measuring the negative pressure created with the pump (± 0.3 % of the full scale of the manometer range). Although rather small, this still represents 3 % of the maximum applied negative pressure for the RTV#1, 5 % for the Ecoflex, 7 % for the RTV#2, and 15 % for the candle gel. These variations could therefore be responsible for part of the total discrepancy. The fact that these variations are higher for the more elastic materials seems to corroborate our assumption that the manometer is responsible for part of the measurement deviation: a precision of ±6 mbar has more negative impact on a softer material. A more accurate manometer could help to reduce this error.

The second area for possible improvement is relative to the camera calibration. If the camera is exactly positioned on top of the mirror, with a point of view at 45 degrees, then the aspired distance measured as seen on the mirror exactly reflects the vertical aspired height. The incorrect orientation or angulation generates an additional perspective distortion that is not taken into account by our very simple manual calibration method which only considers the camera as an ideal pinhole. It should be possible to use a more sophisticated calibration procedure in order to better estimate the intrinsic as well as extrinsic camera parameters.

Because our simple calibration of the extrinsic parameter is mainly based on a pixel/mm ratio that depends on manually picked pixels corresponding to fiducial points inside the vacuum chamber, this ratio can be misestimated. Likewise, manual segmentation of the dome position can introduce errors in the Young modulus estimation. We are currently working on the development of a specific semiautomatic image segmentation algorithm that should decrease the dependency to manual picking. A second camera could be another solution since recording the aspiration from two viewpoints would provide a way to more accurately find the top of the deformation dome.

Finally, the third possible improvement could come from a better fitting of the negative pressure/apex displacement curves resulting from LASTIC by improving the FEM used to build the library of displacements. As shown on Figure 8, the observed strain in the central part of the model is around 20 %, which is in the same order as the strain applied with the tensile devices. Nevertheless, the highest strains appear in the elements that are at the interface with the vacuum chamber edges. Even if this is somehow expected, it can also indicate that the contact resolution is a crucial point and should be carefully monitored. Adding more elements or fine tuning the friction for these elements could also improve the FEA results and consequently the inverse problem resolution. Using a different nonlinear constitutive equation could also help in estimating the full scale of strains of the materials, whereas we are currently only focusing on the linear part of the strains. The FEM can also be improved: it is currently considering a pure symmetrical problem while the LASTIC's vacuum chamber hole is not centered, which might result in non-symmetrical load conditions. Avoiding the errors coming from the FEA could also be achieved by using an analytical solution such as the one presented in [START_REF] Zhao | Comparison of analytical and inverse finite element approaches to estimate cell viscoelastic properties by micropipette aspiration[END_REF]. This kind of approach could give a more precise solution to estimate the Young modulus, avoiding approximation at the interface between the tissue and LASTIC device. Nevertheless, finding the proper analytical solution or an approximation could be a difficult task.

Another source of error in the estimation of the elasticity, not depending on LASTIC, could be coming from the tensile tests and their evaluation. These tests, while being the closest approximation available of the actual characteristics of the different materials, cannot be considered as gold standards. Some errors could be introduced while estimating the stress/strain curve of each material with a Neo Hookean constitutive equation. Furthermore, as shown in Figure 6, these measurements are not always reproducible and could also be the cause for some errors in the Young modulus estimation. This non reproducibility could be coming from the position of the samples between the tensile jaws which can sometimes be nonsymmetrical, see Figure 4, and produce some slight errors in the measurements.

The impossibility to estimate the Young modulus of the RTV#3 by LASTIC measurements is mainly due to the high rigidity of this material: E = 950 kPa as estimated by the tensile test. Even with an aspired volume of 60 ml resulting in a negative pressure of 745 mbar, only 0.4 mm of apex displacement were measured with LASTIC, which is below the fixed threshold for an accurate estimation of the elasticity. Higher values of negative pressure are therefore needed to develop a more consequent aspiration height. Nevertheless, this was not possible with the type of pump used with LASTIC. Considering the range of elasticity found in living tissues, a limitation to linear elasticity lower than 1 MPa seems however to be acceptable.

This chapter presents a device, called LASTIC, to estimate the Young modulus of elastic materials by aspirating an area of their surface and recording the corresponding displacement. With Young modulus overestimated on average by 24 % compared to the value given by standard tensile devices, LASTIC needs to be improved before being used for accurate measurements. Nevertheless, because of the consistency in the overestimation and because of its small size, small cost and capability to be sterilized, LASTIC remains a good tool to estimate the elasticity of living tissues, even during surgery, in order to create a biomechanical simulation of these tissues. The previous version of LASTIC has already been used to characterize tissues of the skin, tongue, and brain, giving interesting insights and a good starting point for simulations.

Figure 2 -

 2 Figure 2 -Complete set up for a measurement: LASTIC is on the silicone sample (here the RTV#2), the manometer is on the left, and the pump is at the top. Luer-Lok connected flexible medical tubes link the three components. They are also all connected to the PC (not shown) that commands the synchronized acquisition.

Figure 3 -

 3 Figure 3 -Images acquired during LASTIC measurements on RTV#2 at different negative pressure values from 0 mbar to -102 mbar. The deformation undergone by the sample can be observed as reflected on the mirror and is measured as the height of the top of the dome.

Figure 4 -

 4 Figure 4 -The Eplexor 500N tensile device applying the displacement load to the RTV#2 silicone (left). Instron® tensile device applying the displacement load to the Ecoflex silicone (right).

Figure 5 -

 5 Figure 5 -LASTIC measurements and the corresponding Young modulus estimated by the FE analysis for (a) RTV#1, (b) RTV#2, (c) Ecoflex, and (d) candle gel.

Figure 7 -

 7 Figure 7 -Summary of the Young modulus found by LASTIC and the tensile devices for each material.

Figure 8

 8 Figure8shows an example of the deformation of a sample using the finite element software ANSYS. ANSYS post processor was used to compute the Von Mises strains resulting from the negative pressure on the FE model. The higher values of the Von Mises strain naturally occur at the interface of the material and the inner edge of the vacuum chamber, where the sample is squeezed while the deformation dome is created. Observed strain in the central part of the model is around 20 % as shown on Figure8for a FEA simulating the RTV#2 silicone (For practical reasons, we used a Mooney Rivlin material for the simulation with E = 32 kPa and υ = 0.49).

Figure 8 -

 8 Figure 8 -Example of aspiration (Von Mises strain) of the RTV#2, using AN-SYS FEA, with E = 32 kPa and υ = 0.49, with a negative pressure of -102 mbar. The aspiration appears at the center of the figure, between the two bars, representing LASTIC vacuum chamber.

Acknowledgments ANR TecSan IDS, Région Rhône-Alpes (projet SIMED, Cluster ISLE), Amandine Dufaug and Luc Maréchal (SYMME), Jérôme Giraud (LSP), Jacques Ohayon (TIMC).