Mathieu Chekroun

Loïc Le Marrec

Bruno Lombard
email: lombard@lma.cnrs-mrs.fr

Joël Piraux

Time-domain numerical simulations of multiple scattering to extract elastic effective wavenumbers

Keywords: ultrasounds, multiple scattering, effective medium, homogenization, numerical methods, finite-difference time-domain schemes, scientific computing, signal processing

Elastic wave propagation is studied in a heterogeneous 2-D medium consisting of an elastic matrix containing randomly distributed circular elastic inclusions. The aim of this study is to determine the effective wavenumbers when the incident wavelength is similar to the radius of the inclusions. A purely numerical methodology is presented, with which the limitations usually associated with low scatterer concentrations can be avoided. The elastodynamic equations are integrated by a fourth-order time-domain numerical scheme. An immersed interface method is used to accurately discretize the interfaces on a Cartesian grid. The effective field is extracted from the simulated data, and signal-processing tools are used to obtain the complex effective wavenumbers. The numerical reference solution thus-obtained can be used to check the validity of multiple scattering analytical models. The method is applied to the case of concrete. A parametric study is performed on longitudinal and transverse incident plane waves at various scatterers concentrations. The phase velocities and attenuations determined numerically are compared with predictions obtained with multiple scattering models, such as the Independent Scattering Approximation model, the Waterman-Truell model, and the more recent Conoir-Norris model.

Introduction

We consider the propagation of elastic waves across a medium containing randomly distributed circular inclusions, the size of which is similar to that of the wavelength. The effective field, which is obtained by averaging the fields in all the possible disordered configurations, corresponds to that of waves propagating in an effective homogeneous medium.

There exists three possible approaches for obtaining the effective wavenumbers (and equivalently, the effective phase velocity and attenuation):

• theoretical approach, based on multiple-scattering models such as the Foldy [START_REF] Foldy | The multiple scattering of waves. I. General theory of isotropic scattering by randomly distributed scatterers[END_REF],

Waterman-Truell [START_REF] Waterman | Multiple scattering of waves[END_REF], and Fikioris-Waterman [START_REF] Fikioris | Multiple scattering of waves. II. Hole corrections in the scalar case[END_REF] models. It provides closed-form expressions useful in practical applications. The main assumption is that the scatterer concentration is low, i.e. typically less than 10 % [START_REF] Martin | Multiple Scattering: Interaction of Harmonic Waves with N Obstacles[END_REF]. At higher concentrations, more sophisticated models developed in acoustics by Linton and Martin [START_REF] Linton | Multiple scattering by random configurations of circular cylinders: second order corrections for the effective wavenumber[END_REF][START_REF] Martin | Multiple scattering by random configurations of circular cylinders: weak scattering without closure assumptions[END_REF], and extended to elastodynamics by Conoir and Norris [START_REF] Linton | Multiple scattering by random configurations of circular cylinders: second order corrections for the effective wavenumber[END_REF][START_REF] Martin | Multiple scattering by random configurations of circular cylinders: weak scattering without closure assumptions[END_REF][START_REF] Conoir | Effective wavenumbers and reflection coefficients for an elastic medium containing random configurations of cylindrical scatterers[END_REF][START_REF] Conoir | Generalization of the Waterman and Truell formula for an elastic medium containing random configurations of cylindrical scatterers[END_REF], are required. But, to our knowledge, no rigorous error estimate is available, and the limits of validity are not accurately known; • experimental approach [START_REF] Derode | Influence of correlations between scatterers on the attenuation of coherent wave in random medium[END_REF]. It introduces no limitation about the concentration of scatterers, but it is very difficult to control accurately the various parameters involved: positions and geometries of scatterers, values of the physical parameters; • numerical approach [START_REF] Lunéville | Mercier Finite element simulations of multiple scattering in acoustic waveguides[END_REF][START_REF] Schubert | Numerical time-domain simulation of diffusive ultrasound in concrete[END_REF]. It allows a simpler control of the parameters and is fast; in practice, however, specific tools are required to perform efficient numerical computations and to render the numerical artifacts much smaller than the quantities of physical interest.

The first aim of the present paper is to describe a numerical approach of this kind.

The second aim is to highlight the efficiency of this methodology to explore the validity domain and limitations of analytical multiple scattering models [START_REF] Kawahara | Attenuation and dispersion of antiplane shear waves due to scattering by many two-dimensional cavities[END_REF].

For this purpose, we will proceed as follows. In section 2, the problem of obtaining random configurations is discussed; naive algorithms converge slowly when the scatterer concentrations are greater than 40 %. The statistical behavior of the configurations is determined by performing a detailed analysis of the radial distribution function. In section 3, time-domain numerical methods are introduced. Elastodynamic equations are integrated using a high-order finite-difference timedomain scheme whose numerical artifacts are known in the case of a homogeneous medium. The discretization of the interfaces between the host matrix and the scatterers is a key issue: special care has to be taken here to prevent the interfaces to introduce large numerical artifacts for physical, geometrical and numerical reasons [START_REF] Lombard | Modélisation Numérique de la Propagation et de la Diffraction d'Ondes Mécaniques[END_REF]. In section 4, signal-processing tools are applied to the simulated data, yielding the effective wavenumbers.

This numerical method is applied to a simple model of concrete, consisting of mortar containing composite inclusions [START_REF] Schubert | Numerical time-domain simulation of diffusive ultrasound in concrete[END_REF]. In section 5, numerical experiments are performed at various inclusion concentrations (ranging from 3 % to 60 %), with longitudinal and transverse incident plane waves. Studies are performed to ensure that the averaged field obtained from a finite number of disordered configurations is representative of the theoretical effective field. In section 6, wavenumbers are extracted from the simulated data. Comparisons are made with multiple-scattering models in terms of the concentration and the adimensional frequency. In particular, the advantages of recent developments over the traditional Waterman-Truell model are confirmed. In section 7, conclusions are drawn and some future lines of research are described. Technical details about the computation of theoretical wavenumbers are given in the appendix A.

Random configurations

Algorithm

Wave propagation is investigated in an infinite medium consisting of a matrix containing circular inclusions with a constant radius a, in the xy plane. Matrix and inclusions are in perfect bonded contact; they both consist of linear elastic isotropic homogeneous media.

In practice, the time-domain numerical simulations are performed in a bounded computational domain 1). A minimum exclusion distance 2 a + ξ between the centers of the scatterers is required by the numerical methods; ξ increases with the mesh size (section 3.1). The x and y coordinates of the centers of circles

[X 1 , X 2 ] × [Y 1 , Y 2 ]. For this purpose, N scatterers are introduced into the rectangular subdomain D = [X 1 , X 2 ] × [Y inf , Y sup ],
Y 1 < Y inf < Y sup < Y 2 (figure
C i (i = 1, • • • , N ) are uniformly distributed in [X 1 , X 2 ] and [Y inf + a, Y sup -a].
Lastly, periodicity of the configuration is imposed along the x-axis, where the period is

X 2 -X 1 .
Various methods to simulate the C i may be found in the literature [START_REF] Siqueira | Numerical simulation of scatterer positions in a very dense medium with an application to the two-dimensional Born approximation[END_REF]. Here we propose two algorithms:

Algorithm 1. ⊲ choose C 1 randomly in D; ⊲ for i = 2 to N do -choose randomly C i in D; -if C i C j ≥ 2 a + ξ (j = 1, • • • , i -1) then C i is kept; -otherwise choose another C i .
Algorithm 1 is very simple and gives quasi-uniform distributions. However, poor convergence is obtained at surface concentration φ greater than 30 %, especially with large ξ. Surface concentrations φ greater than 50 % are beyond the reach of this algorithm ( [START_REF] Torquato | Random Heterogeneous Materials, Microstructure and Macroscopic Properties[END_REF], p.67), even with ξ = 0.

Algorithm 2. ⊲ a compact hexagonal packing pattern consisting of N circles is initially introduced into D. The side length of each hexagon is 2 a + ξ; ⊲ repeat until sufficiently uniform distributions are obtained:

-for i = 1 to N do • perturb the position of C i ; • if C i C j ≥ 2 a + ξ (j = 1, • • • , N , j = i) then C i is kept.
The mean values and variances of the C i coordinates are measured at each iteration. When their third decimal value no long varies, then the perturbation process is stopped. Configurations thus-obtained closely resemble to uniform distributions [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]. Algorithm 2 can be used to reach surface concentrations up to roughly 66 %.

In practice, we recommend the use of algorithm 2 whatever the φ. In order to obtain a sufficiently large number of disordered patterns, the algorithm selected is applied N times, which gives N independent configurations for each incident wave at each scatterer concentration (section 2.2).

Radial distribution function

In this section, we examine numerically whether the above algorithms give uniform distributions. For this purpose, let us take the normalized radial distribution function (RDF)

g(r) = p(r) N n 0 ,
where p is the conditional probability (appendix A), and n 0 is the number of scatterers per unit area. In an infinite statistically homogeneous domain, g → 1 as r → ∞. The RDF is calculated numerically by counting the number n(r) of inclusion centers present in a circular ring with radius r and thickness ∆r:

g(r) = n(r) n 0 1 2 π r ∆r , 2 a + ξ ≤ r ≤ r max , (1)
where r max , the distance to the nearest boundary of D, is used in order to prevent bounding effects. In practice, we take ∆r = a / 20. The RDF is calculated for each inclusion in the simulation domain in order to obtain a representative value.

In the case of dilute media, the number N of inclusions is too low to obtain a smooth curve, and it is not possible to determine the typical behavior of the RDF. The number of configurations N is increased until a standard deviation on g(r) of around 5% is obtained when 7 ≤ r/a ≤ 10: in this range, the RDF is stabilized at 1. The parameters of these calculations are given in table 1. As can be seen in figure 2, the RDF depends greatly on the concentration of the inclusions. At low concentration (φ 10 %), the RDF can be satisfactorily approximated by a Heaviside function if N is sufficiently large. In other words, the conditional probability is uniform in this case. At higher concentrations, the local density of the neighbors in the vicinity of a given inclusion is increased. If φ 30 %, this local increase is proportional to φ, as predicted by the virial expansion [START_REF] Varadan | Scattering and attenuation of elastic waves in random media[END_REF], whereas the conditional probability is uniform if r 4 a. With more densely packed media (φ 40 %), attenuated oscillations occur periodically. As the concentration increases, the oscillations occur farther away from the inclusion, and the period decreases. Our computation of theoretical effective wavenumbers will take into account the distance of exclusion between scatterers (appendix A). But incorporation of non-uniform RDF is beyond the scope of this paper; references on that topic may be found in [START_REF] Siqueira | Numerical simulation of scatterer positions in a very dense medium with an application to the two-dimensional Born approximation[END_REF][START_REF] Caleap | Coherent acoustic wave propagation in media with pair-correlated spheres[END_REF].

φ
The random arrangement of circular inclusions has been studied in detail in many contexts. The behavior obtained here is in agreement with the results of statistical analyses. For further information about these distributions, see [START_REF] Torquato | Multifunctional composites: optimizing microstructures for simultaneous transport of heat and electricity[END_REF]. The most compact 2D-arrangement corresponds to a hexagonal lattice. The RDF of this crystal is composed of Dirac distributions. As shown in figure 2 (right), the positions of the peaks in the simulated RDF do not correspond to the position of the Dirac distributions in the lattice. Even in the case of densely packed media, the heterogeneous structure cannot be approximated by a pseudo-periodic medium.

Time-domain simulations

Integration of elastodynamic equations

A velocity-stress formulation of elastodynamics is followed. The physical parameters are the density ρ, the speeds of longitudinal waves c L and of transverse waves c T . The unknown are the horizontal and vertical velocity (v x , v y ) and the independent components of the stress tensor (σ xx , σ xy , σ yy ). One has to solve the first-order linear hyperbolic system

∂ ∂ t U + A ∂ ∂ x U + B ∂ ∂ y U = 0, (2)
where U = (v x , v y , σ xx , σ xy , σ yy ) T , and A and B are 5 × 5 matrices depending on the physical parameters. The system (2) is solved on a uniform Cartesian grid of N x × N y nodes, with mesh sizes ∆x = (X 2 -X 1 ) / N x and ∆y = (Y 2 -Y 1 ) / N y , and a time step ∆t. In practice, ∆x = ∆y. An explicit two-step finite-difference ADER (Arbitrary DERivatives) scheme is used, giving fourth-order accuracy in both space and time [START_REF] Schwartzkopf | Fast high order ADER schemes for linear hyperbolic equations[END_REF][START_REF] Käser | An arbitrary high order discontinuous Galerkin method for elastic waves on unstructured meshes I: the two-dimensional isotropic case with external source terms[END_REF]. With this scheme, the minimal extra distance between two scatterers is ξ = 3 ∆x (section 2.1). The CFL limit of stability is

θ = c max ∆t ∆x ≤ 1, (3)
where c max is the maximum speed of the waves in the domain. A plane wave analysis of this scheme has been performed in the case of a homogeneous medium [START_REF] Lombard | Modélisation Numérique de la Propagation et de la Diffraction d'Ondes Mécaniques[END_REF], in terms of θ and G = ∆x/λ, G ∈]0, 0.5], where λ is the wavelength. The maximum artifacts are obtained when the direction of the propagation coincides with the grid axes, that is in the case of 1-D configurations. In this case, one has

q(θ, G) = 1 - 2 π 4 15 (θ 2 -1) (θ 2 -4) G 4 + O(G 6 ), α(θ, G) = 4 π 6 9 θ (θ 2 -1) (θ 2 -4) G 6 + O(G 8 ), (4)
where q is the ratio between the exact and discrete phase velocities, and α is the discrete attenuation [START_REF] Strickwerda | Finite Difference Schemes and Partial Differential Equations[END_REF]. The relations (4) have crucial effects on the accuracy of the simulations, because they bound the numerical artifacts in homogeneous media.

Discretization of interfaces

Three classes of drawbacks are classically associated with interfaces in finitedifference schemes. First, since the geometrical description of arbitrarily-shaped interfaces is poor, spurious diffractions are generated. Secondly, since the jump conditions are not enforced numerically, convergence may occur towards a nonphysical solution. Lastly, non-smoothness of the solution across interfaces decreases the accuracy of the scheme, leading to spurious oscillations and even to instabilities. These three drawbacks increase with the scatterer concentration and preclude the use of simulations as metrological tools in highly heterogeneous media. An alternative strategy consists in using numerical methods with unstructured meshes, such as finite-element methods, Galerkin discontinuous methods, and spectral element methods [START_REF] Lunéville | Mercier Finite element simulations of multiple scattering in acoustic waveguides[END_REF][START_REF] Komatitsch | Introduction to the spectral-element method for 3-D seismic wave propagation[END_REF]. But the computational cost of these methods would be much higher due to the meshing, and the stability condition ( 3) is penalized.

To overcome these drawbacks, we use a r-th order immersed interface method [START_REF] Lombard | Modélisation Numérique de la Propagation et de la Diffraction d'Ondes Mécaniques[END_REF][START_REF] Lombard | Numerical treatment of two-dimensional interfaces for acoustic and elastic waves[END_REF]. This numerical method modifies the ADER scheme at grid points close to the interfaces, based on the jump conditions up to the r-th order, the elastodynamic equations, and the Beltrami equations. This procedure associates the efficiency of Cartesian grid methods and the accuracy of an interface meshing. The work is mainly carried out during a preprocessing step, before the numerical integration. At each time step, O(L / ∆x) matrix-vector products are done, where L is the total perimeter of the interfaces, and the matrices are small, typically 5 × 100. The results are then injected into the scheme. After optimizing the codes, the additional CPU time required by the immersed interface method can be made negligible in comparison with the CPU time required by the scheme (less than 1%).

Intensive computing

To obtain reliable effective wavenumbers values, the numerical methods used must meet the following specifications:

• large computational domains, such as grids consisting of 10 4 × 10 4 nodes, involving 10 GB of data; • long integration times, consisting for example of 10 4 time steps;

• several simulations, so as to increase the number of independent disordered configurations, to N = 3, for example; • processing a large number of scatterers with the immersed interface method: 1500 interfaces when φ = 48 %, for example; • performing many simulations in the parametric studies, e.g. in terms of the scatterer concentration or the incident wave polarizations.

To meet these specifications, the computer codes are parallelized. Domain decomposition is performed in the x-direction, associating each slide with one computational process. All the slices have the same size and contain approximately the same number of inclusions, so that the computational cost of each process is roughly the same. After each time step, data are exchanged between neighboring processes.

In practice, the simulations presented in sections 5 and 6 were performed on a cluster of 4 PC bi-processor quadricores, amounting to 32 processes. The optimum speed-up 32 was reached. The communication time between processes was negligible in comparison with the computational cost of each process. After parallelizing, the configurations investigated in section 5 required 1 hour of preprocessing (due to the use of the immersed interface method) and 24 hours of integration (due to the use of the ADER scheme).

Data processing 4.1. Numerical coherent field

At each time step, the components of U have to be stored inside the subdomain containing the inclusions. For this purpose, a uniform network consisting of N l lines and N c columns of receivers is placed in the subdomain D. The position of the receivers is given by (

x i = X 1 + i ∆ c , y j = Y inf + j ∆ l )
, where i = 0, . . . , N c -1 and j = 0, . . . , N l -1. The bottom and top receivers are sufficiently far from the sides of the computational domain to prevent spurious effects from being recorded. These columns of receivers are visible in figure 5-(a).

The acquisition setup has to meet some specifications in order to prevent the occurrence of aliasing and low resolution problems [START_REF] Forbriger | Inversion of shallow-seismic wavefield: I. Wavefield transformation[END_REF]. Aliasing occurs when the distance ∆ l is larger than the shortest wavelength under consideration, while the resolution is limited by the total length N l ∆ l of the acquisition setup (i.e. the distance between the first and last receivers).

The field recorded on each array (each column of receivers) corresponds to a field propagating in a given disordered configuration. Summing the time histories of these N c arrays in all the N simulations gives a coherent field propagating in the y direction. The relationship between the coherent field and the effective field will be discussed in section 5.3.

The polarization of the coherent field is the same as the polarization of the incident plane wave. In the following section, we will therefore deal with either a L (longitudinal) or a T (transverse) incident wave.

Extraction of the coherent wavenumbers

The coherent phase velocity c(ω) is computed by applying a pω transform to the space-time data on the coherent field, where p is the slowness of the waves (p = 1 / c) and ω is the angular frequency [START_REF] Mcmechan | Analysis of dispersive waves by wave field transformation[END_REF][START_REF] Mokhtar | Seismic velocity and Q model for the shallow structure of the Arabian shield from short-period Rayleigh waves[END_REF]. The time Fourier transform of the coherent field s(y j , ω) is denoted by

s(y j , ω) = A(y j , ω) e -i ω p0(ω) yj , (5)
where A(y j , ω) is the amplitude spectrum at y j , and p 0 (ω) needs to be determined.

A p-ω stack quantity ŝ(p, ω) is then defined as

ŝ(p, ω) = Nl-1 j=0
A(y j , ω) e i ω(p-p0(ω)) yj .

ŝ(p, ω) is computed at several p values. Given ω, the maximum value of the modulus |ŝ(p, ω)| is reached at p = p 0 (ω) = 1 / c(ω). The phase velocity dispersion curve is then obtained by taking the maximum locus on the 2-D map |ŝ(p, ω)|. An error estimate is also deduced [START_REF] Herrmann | Computer programs in seismology[END_REF].

Let us now examine the attenuation of the coherent field. This quantity is estimated from the decrease in the amplitude spectrum of the coherent field during the propagation of the waves. In the frequency domain, the amplitude A(y j , ω) in ( 5)-( 6) is assumed to satisfy an exponential decay with distance

A(y j , ω) = A 0 (ω) e -α(ω) (yj -y0) , (7)
where A 0 (ω) is the amplitude of s at the first receiver located at the offset y 0 , and α(ω) is the attenuation factor. Although two points suffice to be able to calculate α, the slope of a least-squares fit of ln(A(y j , ω)) over the whole range of reception gives more accurate results. An error estimate is also deduced. In the first test, transverse wave propagation was simulated in a homogeneous 1-D cement matrix. The dispersion and attenuation measured were due only to numerical artifacts occurring in the ADER scheme. Comparisons between the theoretical (4) and measured dispersion and attenuation values is made in figure 3. The error between the theoretical and measured curves is less than 10 -3 % in the frequency range of interest. The signal processing tools used and the acquisition setup chosen are therefore suitable for accurately assessing the dispersion and the attenuation, and the risk of adding significant signal processing artifacts is thus avoided. In the second test, the wave propagation was simulated in a 2-D medium with a single scatterer centered at (0, 0). The source was the plane compressional wave described in section 5. The diffracted fields were stored in the receiver R at (-0.25, -0.25). Figure 4 gives a snapshot of the stress σ xx after 600 time steps (a) and compares the numerical and analytical values of σ xx during 1400 time steps (b). The exact solution was computed by performing standard Fourier-Bessel decompositions. The excellent agreement observed confirms the validity of both the ADER scheme and the immersed interface method.

Numerical setup

The numerical method presented in sections 2 to 4 is now applied to some physically relevant configurations. A simple model of concrete is studied, where circular aggregates with a radius a = 6 mm are embedded in a homogeneous cement matrix. The physical parameters are (ρ 0 , c 0,L , c 0,T ) = (2050 kg.m -3 , 3950 m.s -1 , 2250 m.s -1 ) in the cement matrix, (ρ 1 , c 1,L , c 1,T ) = (2610 kg.m -3 , 4300 m.s -1 , 2470 m.s -1 ) in the aggregates.

A parametric study is performed in terms of the concentration, from φ = 3% to 60%. The domain of investigation presented in table 2 is discretized on N x × N y = 7200 × 7200 nodes, hence ∆x = ∆y = 10 -4 m. The CFL number (3) is θ = 0.95, giving ∆t = 2.21 10 -8 s. A third-order immersed interface method is implemented (r = 3 in section 3.2). The source is a plane longitudinal (L) or transverse (T) wave propagating along the y-axis, initially outside the domain D (figures 1 and 5-a). The time evolution of the source is a Ricker with a central frequency 250 kHz. The frequencies of interest range between 50 kHz and 600 kHz.

X 1 X 2 Y 1 Y 2 Y inf Y sup -0.36 +0.36 -0.02 +0.7 +0.2 +0.6
Table 2. Coordinates of the physical domain and those of the subdomain D (section 2), in meters.

The acquisition network contains N c = 100 columns and N l = 400 lines, with the spacing ∆ l = 0.001 m and ∆ c = 0.0072 m, respectively. Each column is an array of receivers that follows the wave propagation in a particular realization of disorder. Performing N = 3 simulations yields 3 × 100 independent disordered configurations. This acquisition setup gives the following bounds on the standard errors: from 2 m.s -1 at f = 50 kHz to 0.2 m.s -1 at f = 600 kHz in the case of the phase velocity, and from 0.05 Np.m -1 at 50 kHz to 0.1 Np.m -1 at 600 kHz in that of the attenuation.

Signals recorded along 2 different arrays logically show different behaviors. Figure 6-a gives the time histories at various receivers along one particular array. A main wave train is clearly visible in each of the time histories, followed by a coda.

As recalled in section 4.1, a coherent signal can be obtained by averaging the signals recorded on the various arrays [START_REF] Derode | Random multiple scattering of ultrasound. I. Coherent and ballistic waves[END_REF]. The v y component is used in the case of an incident L-wave, whereas the v x component is used in that of an incident T -wave (stress components provide the same results). An example of coherent seismogram is presented in figure 6-b. The coda has disappeared, and the main wave train behaves like a plane wave propagating in a homogeneous (but dispersive and attenuating) medium.

Convergence of the coherent field to the effective field

The coherent field is obtained by averaging the signals recorded along the 300 arrays of receivers. If this number is too low, the estimated properties will still be dependent on the configuration of the scatterers encountered. Theoretically, the effective wavenumber can be defined by taking an infinite number of configurations, which means that all the possible configurations of scatterers will be taken into account; but this approach is obviously impracticable. The aim of this paragraph is to show that 300 arrays suffice for estimating the effective field and hence, the effective wavenumber.

In each case (in terms of the density φ of the scatterer and the incident wave β = L, T ), a coherent signal is computed with an increasing number N a of arrays ranging from 1 to 300. The N a arrays are chosen randomly among the 300 available ones, to avoid taking consecutive arrays which are located too near each other in the medium. The properties c β and α β are then evaluated from the averaged signal, and their evolution with N a is then studied at a given frequency. Figure 7 shows how α T evolves with N a , at φ = 36% and f = 300 kHz. As this evolution depends on the N a arrays selected in the averaging procedure, the study is repeated 10 times. The value obtained by summing the fields over the 300 arrays is taken as a reference value for the effective medium: in this case, it was α e,T =3.06 Np.m -1 . In figure 7, the red curve gives the envelope of all 10 curves in gray, corresponding to 2 standard deviations of α T (N a ). As was to be expected, as N a increases, the value of α T tends towards the reference value α e,T . This figure shows that when there are too few configurations, α e,T cannot be accurately assessed; for instance, taking only 50 configurations results in an uncertainty greater than 1 Np.m -1 . Taking ±0.2 Np.m -1 to be an acceptable level of uncertainty for α e,T , the optimum number N opt of configurations requested must be greater than N opt ≃ 210 in the case of the present example.

Figure 8 summarizes the values of N opt obtained in all the cases studied, at several frequencies, using the procedure above. The criterion used to obtain an accuracy of about c e,β was ±5 m.s -1 ; in the case of α e,β , we took ±0.2 Np.m -1 . Except for the higher frequencies (> 400 kHz), only about 20 configurations are required for assessing the effective velocity phase c e,β , whereas greater values of N opt are required for assessing the effective attenuation. If N opt is obtained at a given scatterer concentration, its value will increase with the frequency at a given polarity, and it will be almost twice as high with T -waves as with L-waves at a given frequency. Convergence therefore depends mainly on the size of the wavelength, as the effects of multiple scattering are greater at shorter wavelengths.

The results obtained on α T at high frequencies (500 and 600 kHz) show that the optimum number of configurations N opt was almost 300, which was the maximum number of configurations available, whatever the density of the scatterers. In the present 2 cases, it was not possible to say whether the optimum number of configurations was actually reached, so as to be able to assess the attenuation with a sufficiently high level of certainty. In all the other cases, the c e,β and α e,β values obtained based on 300 scatterer configurations were fully representative of the effective wavenumber. 6. Numerical results

Simulated effective wavenumbers

Results obtained in the numerical simulations with the two polarizations of the incident wave (β = L, T ) and at various scatterer concentrations φ are presented in figure 9. In a first approximation, the effective phase velocity was found to be proportional to the density of the inclusions, increasing monotonically with φ. The effective phase velocity also showed a dispersive behavior, which became more conspicuous as φ increased. This effect was stronger at the lower frequencies (at k 0,L a 1 and k 0,T a 2), where the value of the overall phase velocity was lower than at high frequencies. A maximum value of c e was reached at k 0,L a ≈ 1 and k 0,T a ≈ 2. At higher frequencies, c e,β remained almost constant but small fluctuations are visible, up to 10 m.s -1 at 60% which amount to less than 0.5% of the phase velocity. However, the positions of these local extrema are almost φinvariant. At high frequencies, the mean phase velocity c = (1φ) c 0,β + φ c 1,β is a good approximation in the case of dilute media and an upper limit in that of more densely packed media.

The attenuation is more difficult to explain: contrary to what occurs with the phase velocity, the attenuation does not depend monotonically on the concentration. The frequency dependence shows φ-invariant peaks corresponding to the local maxima of the phase velocity, mainly at k 0,L a ≈ 1 and k 0,L a ≈ 3 in the case of Lwaves, and k 0,T a ≈ 2 and k 0,T a ≈ 5 in that of T waves. In dilute media (φ 20%), the attenuation is proportional to φ. In denser heterogeneous media, the behavior effective phase velocity (L wave) effective phase velocity (T wave) Figure 9. Effective properties obtained with numerical simulations and signal processing tools, with longitudinal (left) and transverse (right) incident waves at various inclusion concentrations φ. Top: phase velocity c e,β ; the horizontal dashed lines give the mean phase velocity c. Bottom: attenuation α e,β where φ = 30% (that is, around the critical threshold mentioned in the discussion) is given by a dashed line; the insert is a zoom in the low-frequency range.

depends on both the inclusion concentration and the frequency range. Three types of overall behaviors were observed:

• around the previously mentioned peaks, α e,β increases with φ;

• at high frequencies between these peaks, α e,β remains at an almost constant value when 30% φ 60%; • at low frequencies, α e,β reaches a maximum at φ ≈ 30%, and then it decreases.

With an incident T -wave and k 0,T a > 7, the behavior is less clear-cut: the attenuation reaches a peak at around φ ≃ 30%. However, the accuracy of these findings was not confirmed in the section 5.3.

Comparison between various theoretical models

In this section, the effective wavenumbers predicted by the Independent Scattering Approximation (ISA), Waterman-Truell (WT) [START_REF] Waterman | Multiple scattering of waves[END_REF], and recent Conoir-Norris (CN) [START_REF] Conoir | Effective wavenumbers and reflection coefficients for an elastic medium containing random configurations of cylindrical scatterers[END_REF] multiple-scattering models are presented. Technical details can be found in the appendix A. The results obtained with the explicit (or analytical) formulation of these models are compared in figure 10 in the case of a longitudinal incident wave, and the same comments apply in the case of transverse waves.

ISA and WT models give similar results. The phase velocities differ only at low frequencies in the case of densely packed media, where both models are inaccurate; we will therefore focus on the WT predictions. As regards the phase velocity, all phase velocity: ISA, WT, CN these models predicted the same overall behavior as the simulations: a monotonic increase with the concentration, the same φ-invariant position of the local extrema, and an increase in the dispersion with the concentration. However, the attenuation given by the ISA and WT models are mainly linear functions of the concentration, contrary to what observed with the simulations: the attenuation occurring in densely packed media was clearly over-estimated. The attenuation predicted by the CN model differed considerably from that obtained with the previous models. First, the attenuation showed a linear dependence on the concentration only at low densities φ 30%, where it reached a maximum (except for the peak at k 0,L a ≃ 1). Secondly, the attenuation was negative in some (low or high) frequency ranges, which is unphysical, but this occurs only at very high concentrations. Similar behavior has already be noticed [START_REF] Caleap | Coherent acoustic wave propagation in media with pair-correlated spheres[END_REF][START_REF] Varadan | Multiple scattering of elastic waves by cylinders of arbitrary cross section. II. Pair-correlated cylinders[END_REF]. These comparisons between models and simulations made it possible to determine the range of validity of each model (see figure 11 for the phase velocity and figure 12 for the attenuation). In view of the signal processing limitations, satisfactorily results were defined as those with an error of less than 5 m.s -1 in the case of the phase velocity, and less than 0.2 Np.m -1 in that of the attenuation. On this basis, the WT model is suitable for dealing only with very dilute media φ 12% for c e , and φ ≤ 6% for α e , whereas the CN model gives accurate results up to inclusion concentrations of 24% in terms of both the phase velocity and the attenuation. This considerable difference is probably attributable to the hole correction occurring in the CN model. As mentioned above, the simulation at φ = 6% was performed using only three realizations of the simulation domain, which resulted in the non-smooth RDF shown in figure 2; however, the CN and WT models both gave excellent results at φ = 6%. This indicates that the hole correction has the most significant effects at φ 10%. At φ ≥ 24%, there is a marked discrepancy between the RDF and the Heaviside step function. A more realistic form of p(r 2 |r 1 ) presented in (A1) might extend the range of validity of the CN model to include higher densities.

The WT and CN models are often referred to as second-order models, because of the second-order Taylor expansion (A5) used to estimate k e,β . However, to our knowledge, this definition has never been justified either theoretically or experi-attenuation: L wave attenuation: T wave

where ke,β is the effective wavenumber obtained in the numerical simulations. The difference between the numerical and the theoretical results are defined by

ε β = ( ke,β a) 2 -(k e,β a) 2 , = a 2 n 0 (δ 1,β -d 1,β ) + (n 0 ) 2 (δ 2,β -d 2,β ) + O(n 3 0 ), (9)
where k e,β is the wavenumber deduced from the WT or CN model. This difference is plotted at various dimensionless frequencies versus the concentration φ = n 0 πa 2 in figure 13. With the WT model, ε β is governed throughout the whole frequency and concentration range by a slope of 2 decades per decade on a log-log scale: ε β therefore depends mainly on δ 2,β -d 2,β . In other words, the WT model predicts d 2,β inaccurately, even at low densities. In fact, there is no reason why the WT model should be preferred to the ISA model, which is only a first-order model. With the CN model, ε β is independent of the inclusion concentration up to φ ≈ 18%, and has a +3 slope on log-log scale. This confirms the accuracy of the second-order coefficient d 2,β obtained with the hole correction. When ε β was calculated using an implicit formulation to obtain the wavenumber (A3), ε β was found to have the same concentration-dependent properties with both Waterman-Truell (L wave)

Waterman-Truell (T wave) the WT and the CN models (×-curves in figure 13). The implicit formulation did not improve the quality of the results, and CN is therefore an intrinsically-secondorder model.

Conclusion

The effective properties of random elastic media were calculated here using purely numerical methods. Combining sophisticated methods of simulation (the fourthorder ADER scheme and the immersed interface method) and signal processing tools yielded reference solutions for both the real and imaginary parts of the effective wavenumbers. With this approach, the accuracy of the simulations does not depend on the scatterer concentration. Maximum computational efficiency is obtained by performing domain decomposition and parallelizing the algorithms. This numerical method was applied in the present study to a 2-D model of concrete. In this case, the numerical simulations confirmed that traditional models (such as the Waterman-Truell model) are valid roughly up to inclusion densities of 10%, whereas the recent Conoir-Norris (extension of Linton-Martin to elastodynamics) is valid up to densities of 25 %. In particular, the present simulations confirmed the validity of the second-order term in the Conoir-Norris model, as previously done in a theoretical study [START_REF] Conoir | Effective wavenumbers and reflection coefficients for an elastic medium containing random configurations of cylindrical scatterers[END_REF]. We hope that a similar approach can be used by researchers to test the validity of their favorite multiple-scattering model [START_REF] Kanaun | Effective medium method in the problem of axial elastic shear wave propagation through fiber composites[END_REF].

The numerical method presented here can be used to handle more complex configurations, such as a granulometry or composite containing scatterers of various shapes. Other constitutive laws could also be introduced, such as viscoelastic laws accounting for dissipative effects [START_REF] Lombard | Numerical modeling of transient two-dimensional viscoelastic waves[END_REF]. Lastly, the possibility of extending the present approach to 3-D configurations is a great computational challenge. Preliminary tests have already been conducted on the numerical methods with fluid scatterers included in a fluid matrix.

Figure 1 .

 1 Figure 1. Computational domain. Subdomain D containing the scatterers. A minimum exclusion distance 2 a + ξ between the centers of the scatterers is assumed, where ξ depends on the mesh size.

Figure 2 .

 2 Figure 2. Radial distribution function g(r) at various concentrations φ, calculated with ∆r = a/20 in (1). Left: RDF at φ = 6%, with N = 3 and N = 100. Right: RDF at various surface concentrations φ, calculated with the parameters given in table 1; the thin vertical lines indicate the positions of the Dirac distributions in the case of a hexagonal lattice.

Figure 3 .

 3 Figure 3. Validation test 1/2. Numerical dispersion (a) and numerical attenuation (b) in a 1D homogeneous medium: analytical values (4) in red lines, and measured values in blue circles.

Figure 4 .

 4 Figure 4. Validation test 2/2. Snapshot and time-history of σxx: comparison between the numerical and the analytical solutions at the receiver R. In (a), the green-red palette and yellow-magenta palette denote L waves and T waves, respectively.

Figure 5 .

 5 Figure 5. Incident L-wave, concentration φ = 42 %, initial instant (a) and after 3000 time steps (b). The vertical columns denote the positions of the receivers.

Figure 6 .

 6 Figure 6. Incident plane T -wave in a medium with 24 % inclusion concentration: (a) signals recorded along an array, (b) coherent signals obtained after summation.

1 Figure 7 .

 17 Figure 7. Evolution of α T in terms of Na at f = 300 kHz, φ = 36 %.

Figure 8 .

 8 Figure 8. Optimum number Nopt of configuration of scatterers for determining c e,β or α e,β at all densities.

Figure 10 .

 10 Figure 10. Effective properties of the longitudinal incident waves obtained with the Independent Scattering Approximation (ISA), Waterman-Truell (WT) and Conoir-Norris (CN) models.

6. 3 .Figure 11 .

 311 Figure 11. Comparison between numerical simulations and modeling predictions of the phase velocity in the case of longitudinal (left) and transverse (right) incident waves at the lower inclusion concentrations (φ ≤ 30%).

Figure 12 .

 12 Figure 12. Comparison between the effective attenuation obtained in numerical simulation and modeling predictions at the lower inclusion concentrations (φ ≤ 30%). Results obtained in the case of longitudinal (left) and transverse (right) incident waves versus the Waterman-Truell (top) and Conoir-Norris predictions (bottom). Solid line: results of numerical simulations; dashed line: modeling predictions. mentally. We will therefore examine this question numerically. Let us assume k2 e,β = k 2 0,β + n 0 d 1,β + (n 0 ) 2 d 2,β + O(n 3 0 ),(8)

Figure 13 .

 13 Figure 13. Error ε β (β = L, T ) versus φ at various dimensionless frequencies (top: WT; bottom: CN; left: longitudinal wave; right: transverse wave). ×: implicit method; •: explicit method.

Table 1 .

 1 Parameters for the RDF calculations (1): surface concentration φ, number of configurations N , number of inclusions N .

	6% 12% 18% 24% 30% 36% 42% 48% 54% 60%
	N 100	30	30	20	10	10	10	3	3	3
	N 573 1145 1719 2292 2865 3438 4011 4584 5157 5730

Waves in Random and Complex Media

Acknowledgements

The authors wish to thank Dr Emilie Franceschini for her helpful comments, and Jessica Blanc for her careful reading of the manuscript.

Appendix A. Effective wavenumbers

As established in [START_REF] Conoir | Effective wavenumbers and reflection coefficients for an elastic medium containing random configurations of cylindrical scatterers[END_REF][START_REF] Varadan | Acoustic, Electromagnetic, and Elastic Wave Scattering. Focus on the T-Matrix Approach[END_REF], the scattered field ψ s recorded at r and due to an inclusion centered at r 1 can be expressed through a linear operator T applied to the exciting field ψ e : ψ s (r; r 1 ) = T (r 1 ) ψ e (r; r 1 ). If the positions of the N inclusions are known, closed-form solution of the problem can be obtained. If, on the contrary, the inclusions are randomly distributed, it is generally attempted to determine the effective field ψ e corresponding to the ensemble average of the positions of all the inclusions. The effective field at one representative inclusion (say the first one) is expressed in terms of the scattering induced by another representative inclusion (say the second). Based on the quasicristalline approximation (QCA) [START_REF] Lax | Multiple scattering of waves. II. The effective field in dense systems[END_REF], the latter scattered field is assumed to be excited by the same effective field as the first inclusion:

where ψ i is the incident field. The probability density p(r 2 |r 1 ) in (A1), which is a pair-correlation function, expresses the probability of finding an inclusion at r 2 , given that an inclusion is placed at r 1 .

Various models have been developed for determining the effective wavenumbers analytically, such as the Waterman-Truell (WT) [START_REF] Waterman | Multiple scattering of waves[END_REF] and Lloyd-Berry (LB) [START_REF] Lloyd | Wave propagation through an assembly of spheres IV. Relations between different multiple scattering theories[END_REF] models, to cite but a few. As pointed out by Linton-Martin [START_REF] Linton | Multiple scattering by random configurations of circular cylinders: second order corrections for the effective wavenumber[END_REF], these models differ mainly in the hole correction, i.e. the assumption made about p in the integration of (A1). The LB model is based on the following assumption:

where a and ξ are defined in section 2.1. To obtain the effective wavenumbers in L and T waves, the effective field can be decomposed into a modal sum, where the modal amplitude depends on the effective wavenumber [START_REF] Conoir | Effective wavenumbers and reflection coefficients for an elastic medium containing random configurations of cylindrical scatterers[END_REF]. Introducing this form into (A1) and using (A2), the problem reduces to searching for the non-trivial solution of the infinite linear system [7, equation (31)]:

where I is the identity matrix, T is a matrix defined in [7, equation (37)], and

(A4) In practice, the modal sum in (A3) is truncated. The effective wavenumbers k e satisfying (A3) are associated with k e,L and k e,T . However, searching k e,L and k e,T is an intricate and time-consuming process. Another explicit but approximate form can be obtained using Taylor expansions:

where δ 1,β and δ 2,β are defined in [7, equations (62a) and (62b)] in the Conoir-Norris (CN) model. Let us now consider the WT model, which is based on

where (r 2r 1 ) • e y is the distance between two inclusions in the direction of propagation (e y here) of the incident wave. Based on this hypothesis, the system will have the same form as (A3), but M will be different, and this parameter is deduced from [8, equations (12-13)]:

The approximate form (A5) can still be used with the WT model. The same expression for δ 1,β in (A5) can be used as with the LB model, but the definition for δ 2,β in this case is that given in [5, equation (4)]. Note that ISA matches with (A5) to the first order: k 2 e,β = k 2 β + n 0 δ 1,β .


C Couratin

Ph Velten

X Fléchard

E Liénard

G Ban

A Cassimi

P Delahaye

D Durand

D Hennecart

F Mauger

A Méry

O Naviliat-Cuncic

Z Patyk

D Rodríguez

K Siegień-Iwaniuk

J-C Thomas

First measurement of pure electron shake-off in the β decay of trapped 6 He + ions

Keywords: numbers: 23.40.-s, 31.15.ac, 34.50.Fa, 37.10.Ty

The electron shake-off probability of 6 Li 2+ ions resulting from the β -decay of 6 He + ions has been measured with high precision using a specially designed recoil ion spectrometer. This is the first measurement of a pure electron shake-off following nuclear β decay, not affected by multielectron processes such as Auger cascades. In this ideal textbook case for the application of the sudden approximation, the experimental ionisation probability was found to be P exp so = 0.02339(36) in perfect agreement with simple quantum mechanical calculations.

Electron shake-off (SO) and shake-up (SU) are fundamental atomic processes in which a bound electron is excited into the continuum or in a new orbital, resulting from a sudden change of the central potential. This monopole ionization or excitation may be due to a modification of the nuclear charge, like in nuclear β decay, nuclear electron capture, internal conversion and alpha decay [1], or to the creation of a vacancy in an atomic inner shell induced by collisions with charged particles [2] or by photo-ionization [3,4]. The probabilities of these processes can be calculated in the framework of the sudden approximation (SA). The accuracy of the calculation depends on how fast the central potential changes as compared with the relaxation time of electrons in the new core potential. Nuclear β decay offers ideal conditions to test such calculations since the change in the central potential occurs in less than 10 -18 s, which is the transit time of the emitted β particles through the orbital electron cloud.

The first calculations of SO and SU probabilities following β decay [5][6][7] used hydrogen like wave functions. More sophisticated calculations were performed using numerical self-consistent wave functions for many-electron atoms [8]. The comparison between calculations and experiments is usually difficult [9] since secondary processes, like the emission of Auger electrons, contribute to the final charge state of the daughter ions. The simplest case investigated so far for such comparisons was the SO following the β decay of 6 He atoms [10]. The single ionization probability of the daughter 6 Li + ions was there found in good agreement with a former measurement [11] but the double ionization probability was overestimated by almost one order of magnitude. This stresses the difficulties of treating such systems, even with only two active electrons.

For radioactive species with one active electron, such as 6 He + ions, electron-electron correlations and secondary relaxation processes are absent, leaving only two possible mechanisms for the daughter ionization. The dominant one, the electron SO, is caused by both the rapid change of the nuclear charge and the sudden recoil velocity acquired by the daughter nucleus resulting from β decay. In 6 He decay, the recoil energy can reach 1.4 keV so that the recoil effect can have a sizable impact on the SO probabilities [8]. The second ionization mechanism is a direct collision, in which the β particle knocks out the bound electron. Its probability depends on the energy of the β particle, (E max = 3.5 MeV), as compared to the 54.4 eV electron binding energy [5]. The direct collision is thus expected to have a very small contribution in the 6 Li 2+ ionization following the 6 He + β decay. 6 Li 2+ ions produced after β decay offer therefore an ideal case in which simple quantum mechanical calculations can be performed in the SA with analytic wave functions. Recent progress in the production and manipulation of exotic nuclei associated with new trapping techniques allow refined investigations of the β decay of stored ions [12][13][14]. We present here the first measurement of the SO probability following the β decay of hydrogen-like 6 He + ions. This is a unique system for a comparison with theoretical predictions, which has motivated a careful look at possible systematic effects that could impact the detection of recoiling ions.

The experimental setup has been described in detail elsewhere [13,15,16]. Only the main features of the apparatus are presented here along with the modifications which were necessary for the purpose of this measurement. The experiment has been carried out at GANIL, Caen, France. The radioactive 6 He nuclei were produced at the SPIRAL target-ECR ion source system. After mass separation the 6 He + ions were guided at 10 keV through the LIRAT low energy beam line up to the entrance of the LPCTrap apparatus. At this point, the typical 6 He + beam intensity was 10 8 s -1 . The first stage of the apparatus is a Radio Frequency Cooler and Buncher (RFQCB) [17] for the beam preparation. This linear Paul trap is mounted on a high voltage platform whose voltage was set to decelerate the ions down to 50 eV. The system was filled with H 2 buffer gas, at a pressure of 7 × 10 -3 mbar, to cool down the ions below 1 eV. The 6 He + ions from LIRAT were continuously injected in the RFQCB and accumulated into bunches close to the exit. The cooled bunches were then extracted at a repetition rate of 5 Hz and reaccelerated towards the measurement transparent Paul trap using a pulsed cavity located 12 cm downstream from the exit of the RFQCB. The ions were transported between the two traps with a kinetic energy of about 1 keV and were decelerated down to 100 eV by a second pulsed cavity located at the entrance of the measurement Paul trap. The ions were confined in this trap (Fig. 1) by a 1.15 MHz RF voltage of 120 V pp applied continuously to the two inner rings. The intermediate rings were set to ground potential and the outer rings were set at a voltage of 12 V to minimize trapping losses. During the experiment, up to 2×10 4 6 He + ions were successfully trapped in the measurement trap at each injection cycle, which corresponds to an overall transport and trapping efficiency of 10 -3 . H 2 buffer gas, at a pressure of 4 × 10 -6 mbar, was also used in the trapping chamber to further cool down the trapped ions. The β particles and the recoiling ions resulting from the β decay of the trapped 6 He + ions were detected in coincidence using detectors located around the trap (Fig. 1). The β telescope, composed of a 300 µm thick double sided silicon strip detector (DSSSD) followed by a plastic scintillator, provides the position and the energy of the incoming β particles. The signal from the plastic scintillator triggers the acquisition system and defines the reference time for a decay event. Recently, a new recoil ion spectrometer has been built to separate the charge states of the recoiling ions. Ions emitted towards the recoil spectrometer cross a first collimator through a 90% transmission grid set at ground potential. They are then accelerated by a -2 kV potential applied to a second 90% transmission grid mounted at the entrance of the free flight tube (Fig. 1). Inside this tube, an electrostatic lens at -250 V allows a 100% collection efficiency of the ions, which are detected with a micro-channel plate position sensitive detector (MCPPSD). A -4 kV voltage applied on the front plate of the MCPPSD ensures a maximum detection efficiency for both charge states, independently of the recoil ion initial kinetic energy [19].

For each detected event, the energy and position of the β particle, the time of flight (TOF) and position of the recoil ion, the time of the event within the trapping cycle, and the phase of the trap RF voltage were recorded. In about 25 ms after injection, the trapped ion cloud has reached thermal equilibrium with a final thermal energy kT ∼ 0.1 eV [18]. After a trapping interval of 150 ms (Fig. 2), the ions were extracted towards a second MCPPSD, located downstream from the trap (Fig. 1) that serves as a monitor of the ion cloud [15]. The following 50 ms period within the cycle (Fig. 2), with no ions in the trap, enables to measure the background for each cycle.

FIG. 2: Coincidence events as a function of time within the trapping cycle. The regions 1 and 2 correspond respectively to the selection of "in-trap" and "out-trap" events.

The procedure applied for the detector calibrations was identical to that described in Ref. [16]. Several conditions then applied to the data: (1) the energy deposited in the plastic scintillator had to be larger than 0.4 MeV;

(2) the signals in the DSSSD must have had a valid conversion corresponding to a minimum ionizing particle; and (3) the signals from the delay lines of the MCPPSD should provide an unambiguous determination of the ion impact position. Two different sets of events were selected using the detection time, T cycle , of events within the trapping cycle (Fig. 2). The 50 ≤ t cycle ≤ 149 ms interval contains decay events from trapped ions having reached thermal equilibrium whereas, in the interval 151 ≤ t cycle ≤ 200 ms, the Paul trap was emptied so that only background events were recorded. Such events, called "out-trap" events, result from the decay of neutral 6 He atoms in the detection chamber. This background contribution was less than 1% of the recorded data. After normalization, the background events were subtracted from those in the first time interval, which contains the events of interest. The resulting TOF spectrum is shown in Fig. 3. The regions labeled 1 and 4 contain only acci- dental events associated with uncorrelated signals from the β detector and from the MCPPSD. These events are used for the accidentals subtraction from region 2, where the distributions corresponding to 6 Li 3+ and 6 Li 2+ recoil ions are clearly visible. The tail observed in region 3 arises from recoil ions neutralized and scattered by metallic surfaces prior to entering the acceleration region of the spectrometer. These events have a minimum TOF of 3.2 µs and are excluded from the selected data.

The data analysis is similar to the one detailed in Ref. [16], and is based on the comparison between the experimental TOF spectrum and two sets of Monte-Carlo (MC) simulated spectra obtained for 6 Li 2+ and for 6 Li 3+ recoil ions. For both sets, the β decay kinematics was accurately incorporated, using the βν angular correlation coefficient predicted by the Standard Model, including radiative corrections terms [20]. The initial positions and velocities of the decaying ions were sampled accordingly to distributions obtained from simulations of the motion of trapped ions in the Paul trap, in the presence of H 2 buffer gas. The 6 Li 2+ and 6 Li 3+ paths from the decay point to the MCPPSD were computed using the SIMION8 software. The paths of the β particles were computed with the GEANT4 toolkit to account for the scattering on the electrodes of the Paul trap and on the detectors. The response functions of the β telescope and of the MCPPSD were also included in the MC simulation. It has previously been shown [21] that the shape of the TOF spectrum strongly depends on the ion cloud temperature and on the ions flight distances. The whole simulation was therefore carried out assuming different cloud temperatures around kT = 0.1 eV, and for different positions of the MCPPSD relative to the trap. After subtraction of background and accidental events, the experimental TOF spectrum was fitted with a linear combination of the simulated spectra corresponding to the two charge states (Fig. 4). The free parameters of the fit are: a global normalization, the electron SO probability, P exp so , the ion cloud temperature, and the distance of the MCPPSD from the trap center. The best fit (Fig. 4) has a chi-square at the minimum of χ 2 = 533 for 523 degrees of freedom, which corresponds to a P-value of 0.37. This indicates a very good statistical consistency between the data and the model. The fit leads to an electron SO prob- ability P exp so = 0.02339 ± 0.00035, where the quoted error is purely statistical at one standard deviation. Comparisons between the simulations and the experimental data have been performed using other observables, like the position profiles in both detectors and the β energy spectra, and also show a good agreement. The main sources of systematic effects, including the calibration and position accuracy of the detectors, the high-voltage power supply accuracy, background subtraction, and charge exchange on H 2 buffer gas, have been investigated. The detection efficiency of the MCPPSD for the two charge states was determined from the gaussian shapes of the associated measured charges [19]. The efficiency for 3 Li 3+ was found to be (0.38 ± 0.05)% larger than for 3 Li 2+ and has been included. Effects giving a contribution larger than 10 -5 to the uncertainty on P exp so are listed in Table I. The methodology used to estimate these effects is detailed in Ref. [16]. The total systematic uncertainty, 7.0 × 10 -5 , is very small compared to the statistical error.

The electron SO probability in an hydrogen-like system is a standard application of the SA and can be calculated following e.g. Ref. [START_REF] Landau | Quantum Mechanics Non-Relativistic Theory[END_REF]

P so = 1 - n,l,m | 1, 0, 0, Z | exp (-i K r) | n, l, m, Z ′ | 2 , (1)
where n, l and m are the radial, orbital momentum and orbital momentum projection numbers respectively, Z = 2 (Z ′ = 3) is the number of protons in the initial (final) system and K is the wave vector of the final TABLE I: Dominant sources of systematic effects along with the size of the correction of P exp so if any (second column), the impact on the error on P exp so (third column) and the method used to estimate the parameters (fourth column).

Source

Corr. (10 -5 ) Error (10 a The size of the corrections are given for indication. The instrumental effects were actually incorporated in the MC fit, so that the value of P exp so obtained from the fit includes these corrections.

system. For the SO following the β decay of 6 He + ions, this leads to

P so = (2.33810 + 0.00412E rec ) × 10 -2 , (2)
where E rec is the recoil energy of the daughter nucleus expressed in keV. The calculations were performed applying 10 4 non-relativistic hydrogen-like wave functions with different values of n, l and m. The reached numerical accuracy was smaller than 10 -8 . Next, a correction to the SA probability was estimated to for the finite duration of the potential change, which also includes the direct collision mechanism. It was assumed that the β particle travels at the speed of light and that its wave function is distributed on a thin spherical shell with radius R = ct. Therefore, the effective perturbation interaction, δV (r, R), between the β particle and the orbital electron has a constant value, e 2 /R, inside the sphere and the Coulomb form, e 2 /r, outside the sphere of radius R. The time dependent Schrödinger equation was solved perturbatively, with an unperturbed Hamiltonian taken for the hydrogen-like 6 Li 2+ ion [START_REF] Landau | Quantum Mechanics Non-Relativistic Theory[END_REF]. The variation in the SO probability was estimated from

δP so = - 1 (hc) 2 R0 0 dR 1 R0 0 dR 2 1s, Z|δV (r, R 1 )δV (r, R 2 )|1s, Z , (3)
where R 0 = 2 1s, Z|r|1s, Z is the radius of a sphere containing 93% of the electron charge. This leads to the value δP so = -20 × 10 -5 , which is a 1% relative correction to the dominant value obtained in the SA. Such a small correction was indeed expected for a fast process like nuclear beta decay. Substituting in Eq.( 2) the mean recoil energy of the events selected in the experimental data analysis (whose error contribution is negligible) together with the correction from Eq.(3) leads to

P so = 0.02322. (4)
The final experimental result obtained in this work, P exp so = 0.02339 ± 0.00036, (5) is in perfect agreement with the theoretical result above and with a previous calculation [10]. The measured value is however inconsistent, by 10σ, with the ionization probability per electron of the K-shell estimated by Feinberg for ǫ = Z ′ /Z = 3/2, W K = 0.01983 [5].

In conclusion, we have reported the first measurement of the SO probability following the β decay of hydrogenlike ions. The decay of 6 He + ions is a unique system that fulfills all conditions to perform a precision comparison between the experimental result and simple, albeit complete, quantum mechanical calculations. Such conditions are: a single active bound electron, a very fast change of the central potential and a resultant pure shake-off process not affected by secondary ionizations like Auger emissions. At the present level of precision, the experimental result was found to be consistent with the theoretical prediction.

We warmly thank Bernard Pons for fruitful discussions that eventually led to the submission of this article. We express our gratitude to G. Darius and M. Herbane for their contributions at an earlier stage of the setup. We are indebted to the LPC staff for their strong involvement in the LPCTrap project and to the GANIL staff for the preparation of a high quality ion beam. K.S.-I. and Z.P. acknowledge the financial support from Narodowe Centrum Nauki (Poland) grant No. 2011/01/B/ST2/05131.

FIG. 1 :

 1 FIG. 1: Top view of the experimental setup. The insert shows the structure of the six stainless steel rings of the Paul trap. See text for details.

FIG. 3 :

 3 FIG. 3: TOF spectrum of recoil ions. The four regions used in the data analysis are delimited by vertical lines. The horizontal gray line indicates the average level of accidental events.

FIG. 4 :

 4 FIG. 4: Upper panel: Fit of the experimental spectrum with the MC simulations. The range selected for the fit is indicated by the vertical lines. Lower panel: Normalized residuals of the fit.


Arnaud Réveillère

Daniel Monfort

Sophie Lecacheux

Ludovic Grisanti

Héloise Muller

Didier Bertil

Jeremy Rohmer

Olivier Sedan

John Douglas

Audrey Baills

Comparative risk assessments for Guadeloupe: earthquakes and storm surge

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Risk division

Comparative risk assessments for Guadeloupe: earthquakes and storm surge > 17

  Adapted fromUhlemann et al., 2011 Grünthal et al., 2006, for the city of Cologne > Incomplete but quantitative measure of the disaster Assets estimationconstruction cost / m 2

>••

  Overestimated cost. Hypotheses: no reimbursement of light damage (no declaration, insurance excess), signification of the CCR number, %loss -DS relation … The marine submersion of the coastal areas results from the conjugated effects of: • the tide • the atmospheric surge (due to wind and low atmospheric pressure) • the waves set-up (local elevation of the mean sea level due to wave breaking) The maximum level reached by the water is the run-up Probabilistic storm surge hazard: Validation of the economic losses (Les Saintes) • Finalization of the probabilistic storm surge hazard • Mapping inundated areas and estimating storm surge losses > The research leading to these results has been carried out in the frame of the MATRIX Project, funded by the European Commission's Seventh Framework Program [FP7/2007-2013] under grant agreement n° 265138. The BRGM research project RISCOTE has also cofunded the storm surge part.

Building value per net floor area per construction type, Construction type Construction cost (€ per net floor area) Source (incl. year and location)

	>
	Individual housing	1127	EPTB 2010 for Overseas
			Territories
	Collective housing	1000	Based on the average social
			housing price in Guadeloupe
	Shelter	600	Assumption
	Industrial / large business	1390	Based on market price

based on the construction value rather than on the market value of

the building similarly to

Kleist & al., 2006;[START_REF] Fema | HAZUS Earthquake Loss Estimation Methodology: User's Manual -Federal Emergency Management Agency[END_REF][START_REF] Dutta | A mathematical model for flood loss estimation[END_REF]

2003 > Local and recent data are used, if possible > 3 Assets estimation -average surface Living Surface per occupency type Occupency type Average living space per dwelling Source (incl. year and location)

	Traditional housing	71 m 2	INSEE, 2006, for Guadeloupe
	Recent private housing	101 m 2	INSEE, 2006, for Guadeloupe
	Villa	150 m 2	Assumption
	Collective housing	65 m 2	INSEE, 2009 for France & DGAFP,
			for France
	Makeshift shelters	50 m 2	Assumption
	Industrial buildings	300 m 2	CCI, for Guadeloupe

> 4

Assets estimationoverall methodology

Structure replacement value per SDRS type Living Surface per occupency type Building value per construction type, per net floor area Occupency type per SDRS type Ratio of living space on net floor area Net floor area per SDRS type Content replacement value per SDRS type Content estimation relatively to the occupation type

> 5 « SDRS type »: building vulnerability typology defined by the « Regional Scenario for Seismic Risk » study & surveys. Cf.

[START_REF] Bertil | Scénarios Départementaux de Risque Sismique en Guadeloupe[END_REF]

Assets estimation -results > 6

SDRS

type Building stock in Guadeloupe Replacement value

	Loss estimation per Damage State
	> Per vulnerability type
						Building	Content
	Name	Description	Nb of dwellings Share	€ per dwelling)	(€ per dwelling)
	HABFOR	makeshift shelter	6 424	3%	37 500	18 750
	MCPIER	stone houses	609		0.3%	100 021	50 011
	CASTRA … > Structural repairs traditionnal houses (wood) 15 710 … …	7% …	100 021 … > Content replacement 50 011 …
	Economic sector Private housing Industry > Total assets EMS-98 DS Structure damage ratio 0 0% 1 0-1% 2 1-20% 3 20-60% 4 60-100% 5 100%	Guadeloupe exposed assets Total (G €) k€/hab share 36.5 91.2 77% 1.1 2.7 2% Central damage EMS-98 DS factor 0% 0 0.5% 1 10% 2 40% 3 80% 4 100% 5	Replacement ratio 0% 1 % 2 % 12 % 25 % 50 %
	Commerce & service From Tyagunov et al. (2006)	5.7	14.3	12% Adapted from FEMA (2003)
	Schools & hospitals for German buildings typology 4.4	11.0	9% for US buildings typology
	Others (roads, energy &		0.0	-	0%
	water supply, etc.)
	Total		47.6	119.1	100%
							> 7