Raphaël F Garcia

Jeannine Gagnepain-Beyneix

Sébastien Chevrot

Philippe Lognonné

Very Preliminary Reference Moon Model

Keywords: Moon, Seismology, Internal Structure, Core, Seismic Body Waves

The deep structure of the Moon is a missing piece to understand the formation and evolution of the Earth-Moon system. Despite the great amount of information brought by the Apollo passive seismic experiment (ALSEP), the lunar structure below deep moonquakes, which occur around 900 km depth, remains largely unknown. We construct a reference Moon model which incorporates physical constraints, and fits both geodesic (lunar mass and polar moment of inertia, and Love numbers) and seismological (body wave arrivals measured by Apollo network) data. In this model, the core radius is constrained by the detection of S waves reflected from the core. In a first step, for each core radius, a radial model of the lunar interior, including P and S wave velocities and density, is inverted from seismic and geodesic data.

In a second step, the core radius is determined from the detection of shear waves reflected on the lunar core by waveform stacking of deep moonquake Apollo records. This detection has been made possible by careful data selection and processing, including a correction of the gain of horizontal sensors based on the principle of energy equipartition inside the coda of lunar seismic records, and a precise alignment of SH waveforms by a non-linear inversion method. The Very Preliminary REference MOON model (VPREMOON) obtained here has a core radius of 380±40 km and an average core mass density of 5200±1000 kg/m 3 . The large error bars on these estimates are due to the poorly constrained S-wave velocity profile at the base of the mantle and to mislocation errors of deep moonquakes. The detection of horizontally polarized S waves reflected from the core and the absence of detection of vertically polarized S waves favour a liquid state in the outermost part of the

Introduction

The radius of Earth's core has been relatively precisely determined more than one century ago from the analysis of seismic data [START_REF] Oldham | The Constitution of the Interior of the Earth as Revealed by Earthquakes[END_REF], only 17 years after the first measurement on a seismogram [START_REF] Von Rebeur-Paschwitz | The Earthquake of Tokio[END_REF]. But more than 40 years of negative attempts and processing efforts have been necessary before the publication of the first reported observations of seismic core phases in the Apollo seismic data [START_REF] Weber | Seismic Detection of the Lunar Core[END_REF]. However, the results of this pioneering paper in which some inconsistency between the discontinuity radii deduced from different boday waves remains, must be confirmed by a more precise analysis which integrates the trade-off between the mantle seismic velocities and the core size, and the constraints on the body wave amplitudes. Note also that Love numbers were only computed a posteriori and not integrated in the model inversion. As the radius of the core is a key parameter for constraining scenarios of Moon formation by a giant impact [START_REF] Benz | The origin of the moon and the single impact hypothesis[END_REF][START_REF] Canup | Origin of the Moon in a giant impact near the end of the Earth's formation[END_REF], its precise determination will greatly improve our understanding of the geometry and dynamics of such giant and catastrophic impact, providing crucial constraints on the state of primordial Earth's mantle. Prior to [START_REF] Weber | Seismic Detection of the Lunar Core[END_REF], none of the published seismic models [START_REF] Toksoz | Structure of the moon[END_REF][START_REF] Nakamura | Structure of the lunar mantle[END_REF]Nakamura et al., , 1982;;[START_REF] Khan | An inquiry into the lunar interior: A nonlinear inversion of the Apollo lunar seismic data[END_REF][START_REF] Lognonné | A new seismic model of the Moon: implications for structure, thermal evolution and formation of the Moon[END_REF][START_REF] Gagnepain-Beyneix | A seismic model of the lunar mantle and constraints on temperature and mineralogy[END_REF][START_REF] Khan | Joint inversion of seismic and gravity data for lunar composition and thermal state[END_REF] have put direct constraints on the lunar core radius, and only indirect constraints were achieved. Those come from the observations of magnetic signals related to the lunar's core magnetic field rejection [START_REF] Russell | Measurements of the lunar induced magnetic moment in the geomagnetic tail -Evidence for a lunar core[END_REF][START_REF] Hood | Initial measurements of the lunar induced magnetic dipole moment using Lunar Prospector magnetometer data[END_REF] and from the inversion of geodetic data [START_REF] Hood | Geophysical constraints on lunar bulk composition and structure -A reassessment[END_REF][START_REF] Mueller | Lunar composition -A geophysical and petrological synthesis[END_REF][START_REF] Kuskov | Constitution of the Moon5. Constraints on composition, density, temperature, and radius of a core[END_REF][START_REF] Kuskov | Geochemical Constraints on the Seismic Properties of the Lunar Mantle[END_REF], including the Love numbers [START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF][START_REF] Khan | Does the Moon possess a molten core? Probing the deep lunar interior using results from LLR and Lunar Prospector[END_REF][START_REF] Khan | Further constraints on the deep lunar interior[END_REF]. Among these constraints, the difference between the spin of the polar axis and the Cassini spin state measured by [START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF] is central for suggesting a liquid core. So far, the suggested values of this lunar liquid core radius range from 250 km to 450 km.

Constraints on the core composition are even much weaker, as the end member values for its size lead to either pure iron core for the smallest values, to ilmenite core for the largest ones (see for example Lognonné and Johnson (2007) for the trade-off of core models verifying both the mantle seismic constraints and the global geodetic ones). The seismically inverted temperature at the bottom of the mantle provided by [START_REF] Gagnepain-Beyneix | A seismic model of the lunar mantle and constraints on temperature and mineralogy[END_REF] or [START_REF] Khan | Joint inversion of seismic and gravity data for lunar composition and thermal state[END_REF] are all lower than the pure iron liquidus temperature, sugesting that some light elements must be present in the core if it is liquid.

These two models however neither include the Love numbers nor core seismic phases, even though these two types of information could directly constrain the size of the core in their inversions.

In this study, we investigate the seismic structure of the lunar interior following a two-step approach. First, we construct a set of acceptable seismic models of the Moon constrained by mass, moment of inertia, Love numbers and arrival times of P and S waves measured by the Apollo passive seismic experiment. Then, we estimate the core radius by detecting core reflected S wave arrivals from waveform stacking methods. The final result is a preliminary reference model for the Moon, that includes the size and average density of the core.

Construction of seismic moon models

This section describes the construction and selection of the best radial seismic lunar models. These models are described by a small number of parameters, respecting a set of a priori geophysical equations and fitting seismological and geodesic observations. The aim is to obtain the best physical models for each core radius and use them to detect core reflected phases.

A priori information

The radial models are constructed in order to fit exactly the lunar mass M o = 7.3458 × 10 22 kg [START_REF] Konopliv | Recent Gravity Models as a Result of the Lunar Prospector Mission[END_REF][START_REF] Goossens | Lunar degree 2 potential Love number determination from satellite tracking data[END_REF]. The crustal model is extracted from [START_REF] Gagnepain-Beyneix | A seismic model of the lunar mantle and constraints on temperature and mineralogy[END_REF] seismic model (see figure 1). Crustal density (ρ c ) below the 1 km thick regolith layer is assumed to be constant and allowed to vary in order to fit geodesic observations. The strong degree one crustal thickness variation is taken into account by assuming a Moho depth of 28 km below the Apollo network (on the near side) in order to fit seismic observations, and an average Moho depth of 40 km [START_REF] Chenet | Lateral variations of lunar crustal thickness from the Apollo seismic data set[END_REF] in order to fit geodesic observations. Studies of lunar gravity and topography constrain both the average crustal thickness and the density contrast between the crust and the mantle.

However, these parameters are correlated. Assuming an average Moho depth of 40 km, the density contrast between crust and mantle is predicted to be about 0.55 [START_REF] Chenet | Lateral variations of lunar crustal thickness from the Apollo seismic data set[END_REF][START_REF] Wieczorek | The constitution and structure of the lunar interior[END_REF]. This value is taken as an a priori for the construction of the radial Moon models.

The lunar mantle is assumed to be homogeneous and without any phase transition. This assumption is strong, but if the composition of the lunar mantle is close to Earth's mantle, no phase change of major silicate constituents is expected in the pressure and temperature ranges of the Moon mantle. In addition, previous investigations of lunar mantle chemistry based on seismic data did not detect any strong chemical composition variation with depth [START_REF] Khan | Joint inversion of seismic and gravity data for lunar composition and thermal state[END_REF], suggesting that the 550 km depth discontinuity, reported in earlier models [START_REF] Nakamura | Seismic velocity structure of the lunar mantle[END_REF][START_REF] Khan | An inquiry into the lunar interior: A nonlinear inversion of the Apollo lunar seismic data[END_REF], is weakly resolved by the Apollo data set (Lognonné and Johnson, 2007). In addition, we assume an adiabatic temperature gradient in the lunar mantle. Even if not fully valid in the middle and upper mantle, where the temperature profile is mainly conductive [START_REF] Gagnepain-Beyneix | A seismic model of the lunar mantle and constraints on temperature and mineralogy[END_REF], the adiabatic hypothesis is expected to generate no more than a few percent of differences in seismic velocities [START_REF] Bina | Seismological Constraints upon Mantle Composition[END_REF] which is smaller than the typical error of lunar seismic models. Assuming homogeneity and adiabaticity of the lunar mantle, we can use the Adams-Williamson equation to derive the density profile of the Moon mantle:

dρ dr = - ρ(r)g(r) Φ(r) (1)
where r is the radius, ρ(r) the volumic mass, g(r) gravity, and Φ(r) the seismic parameter given by:

Φ(r) = V 2 P (r) - 4 3 V 2 S (r) (2)
where V P (r) and V S (r) are respectively the P and S wave velocities. If we assume a state of hydrostatic equilibrium for the lunar interior, pressure variations are described by:

dP dr = -ρ(r)g(r) (3)
where P (r) is pressure.

Consequently, if seismic velocities are known, the above equations can be integrated from top to bottom in order to construct gravity, pressure, and density profiles, using the additional relation between gravity and density:

g(r) = r o r 2 g o - 4πG r 2 o r o r 2 ro r ρ(u)u 2 du (4)
where G is the gravitational constant, r o the average Moon radius (1737.1 km) and g o = GMo r 2 o the surface gravity. This procedure has been used since the early sixties in order to produce density models of Earth's interior assuming a starting value of density ρ o [START_REF] References Alterman | Oscillations of the Earth[END_REF].

In order to link seismic wave velocities to density, a Birch law [START_REF] Birch | Density and Composition of Mantle and Core[END_REF] is used:

V P (r) = a + bρ(r) (5)
where a and b are assumed to be constant because the mantle is assumed to be homogeneous. In addition, because the P-wave to S-wave velocity ratio increases with depth in previous Moon seismic models (Nakamura et al., 1982;[START_REF] Khan | An inquiry into the lunar interior: A nonlinear inversion of the Apollo lunar seismic data[END_REF][START_REF] Lognonné | A new seismic model of the Moon: implications for structure, thermal evolution and formation of the Moon[END_REF][START_REF] Gagnepain-Beyneix | A seismic model of the lunar mantle and constraints on temperature and mineralogy[END_REF], it is assumed to vary linearly with radius according to:

V P (r) V S (r) = A + Br (6)
For given values of density at the top of the mantle (= ρ c + 0.55), Birch law parameters (a and b), and P to S velocity ratio parameters (A and B), a seismic model of the mantle predicting V P (r), V S (r) and ρ(r) can be constructed down to the core.

For a given value of core radius (R core ), and using the lunar Mass, the average core density (ρ core ) is deduced.

The procedure described above allows us to build seismic and density profiles of the lunar interior obeying simple physical laws with only 6 parameters: crust density (ρ c ), Birch law parameters (a and b), P to S velocity ratio parameters (A and B) inside the mantle, and core radius (R core ).

These parameters are gathered in a model vector m. Average core density is computed in order to fit exactly the mass of the Moon. Only models with core density in the range [4.0-10.0] kg/cm 3 , and fitting both Love numbers k 2 =0.0213±0.0025 and h 2 =0.039±0.008 [START_REF] Williams | A scheme for lunar inner core detection[END_REF][START_REF] Goossens | Lunar degree 2 potential Love number determination from satellite tracking data[END_REF] and polar moment of inertia ratio IR=0.3932±0.0002 [START_REF] Konopliv | Recent Gravity Models as a Result of the Lunar Prospector Mission[END_REF] within their error bars will be selected.

In order to construct physical Moon models, the 6 parameters are explored following ranges and steps given in table 1. The range of crustal density (ρ c ) is defined according to previous estimates [START_REF] Gagnepain-Beyneix | A seismic model of the lunar mantle and constraints on temperature and mineralogy[END_REF][START_REF] Chenet | Lateral variations of lunar crustal thickness from the Apollo seismic data set[END_REF]. Birch law parameters are explored between values corresponding to the Earth upper mantle as defined in PREM (a = -7.4146 and b = 4.5872) [START_REF] Dziewonski | Preliminary reference Earth model[END_REF] and IASP91 (a = -10.7346 and b = 5.5743) [START_REF] Kennett | Traveltimes for global earthquake location and phase identification[END_REF] reference Earth models. P wave to S wave velocity ratio is explored around the value obtained for a

Poisson solid ( √ 3 ≈ 1.73). The ranges for origin (A) and slope (B) of its variation with radius are chosen in order to allow either decrease or increase with depth, and variations around values obtained in previous models. Core radius is explored from 250 km to 490 km.

Parameter

Range

Step Crust density (kg/cm 3 ) (ρ c ) [2.6 3.0] -Birch law origin (a) [-10.74 -7.41] -Birch law slope (b) [4.58 5.57] -

V P V S law origin (A) [1.0 2.95] - V P V S law slope (km -1 ) (B) [-0.000662 0.0004] - Core radius (km) (R core )
[250 490] 5

Best radial Moon models

The next step is to select the best radial Moon models that fit both P and S arrival times and geodesic observations. In order to let the core radius as a free parameter, a best radial Moon model will be obtained for each core radius.

Then, for each core radius, an exploration of the five remaining parameters is performed with a Neighbourhood Algorithm (NA) [START_REF] Sambridge | Geophysical inversion with a neighbourhood algorithm-I. Searching a parameter space[END_REF] with ten randomly selected starting models and performing twelve iterations of NA with parameters ns=30 and nr=3. These parameters mean that at each iteration the neighbourhoods of the three best models are explored with 10 new random models. The cost function to minimize is defined as the sum of χ 2 functions of seismic travel times and geodesic observations (Love numbers and polar moment of inertia ratio): However, in order to be fully consistent, the computation of χ 2 seismo requires relocation of all natural seismic events, and correction for P and S wave shallow structure below the seismic stations. Therefore, for each physical model tested, the P and S wave traveltimes given by [START_REF] Lognonné | A new seismic model of the Moon: implications for structure, thermal evolution and formation of the Moon[END_REF] are inverted in order to relocate the events and to obtain P and S wave station corrections of zero average. Only arrival time data with errors smaller than 10 s are kept. A total of 343 P and S travel times from 64 events (8 artificial impacts, 19 meteor impacts, 10 shallow events and 27 deep moonquakes) are used. The inversion is performed with a damped Gauss-Newton algorithm [START_REF] Tarantola | Inverse Problem Theory[END_REF] with starting event location parameters given by [START_REF] Gagnepain-Beyneix | A seismic model of the lunar mantle and constraints on temperature and mineralogy[END_REF] and zero station corrections. The χ 2 seismo value of P and S wave travel times obtained after inversion is computed for each physical model tested.

J(m) = 1 Ns N s i (t obs i -t calc i ) 2 σ 2 i + 1 3 (h obs 2 -h calc 2 ) 2 σ 2 h 2 + (k obs 2 -k calc 2 ) 2 σ 2 k 2 + (IR obs -IR calc ) 2 σ 2 IR (7)
Results of the inversion are summarized in figures 2, 3 and 4. Figure 2 gives an example of the parameter space sampling by NA algorithm for a core radius of 380 km. The sampling by NA algorithm concentrates around a single minimum value inside the volume of parameter space explored. Figure 3a plots the minimum values of the cost function J(m), χ 2 geod and χ 2 seismo obtained for each core radius. The best models for core radius between 300 km and 400 km give a similar fit of the seismological and geodesic observations.

For core radius smaller than 300 km, the fit of seismological observations is slightly degraded. For core radius larger than 400 km, χ 2 geod at minimum value of the cost function increases with radius, reflecting the difficulty to fit geodesic observations with large core radii. Figure 3b plots the variation with radius of the P and S wave velocities and density for the best models obtained for each core radius. A comparison with the starting model of [START_REF] Gagnepain-Beyneix | A seismic model of the lunar mantle and constraints on temperature and mineralogy[END_REF] shows that the fit of geodesic observations favours S wave velocities at the bottom of the mantle higher than previously inferred. Figure 4 presents the correlation between parameters obtained for the best radial models. Birch law parameters are strongly correlated between themselves, but not with core radius. V P V S ratio presents an increase with depth consistent with previous investigations. The parameters describing the variations of V P V S ratio are correlated between themselves, but there is no significant correlation with core radius. The crust density is almost constant around 2.76 kg/cm 3 . It means that, despite large variations of Birch law parameters, imposing the average crustal thickness (40 km) and the density jump at the crust/mantle interface (0.55 kg/cm 3 ) strongly constrains the average crust density in a radial model consistent with geodesic observations. However, an average crust density of 2.76 kg/cm 3 is consistent with our assumptions and with previous estimates [START_REF] Chenet | Lateral variations of lunar crustal thickness from the Apollo seismic data set[END_REF]. The core density is correlated with core radius with high densities for small core sizes and low densities for large core sizes. Overall, geodesic and travel time data poorly constrain the size of the core, which motivated us to look for additional constraints that could come from core reflected phases.

Modelling core reflected phases

This section demonstrates that the core radius can be constrained by the detection of horizontaly polarized S waves reflected on the core (ScSH) and precise the data set that be can be used for such a detection. The data and signal processing methods allowing ScSH wave detections are described in the next section.

Waveform modelling

We focus on SH waves reflected on the core because, if the upper part of the lunar core is fluid [START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF], their energy is completely reflected back to the surface. Ray theoretical amplitudes of body waves are computed including geometrical spreading and attenuation (see last table of the paper for values of quality factor). As shown in figure 5a, ScSH amplitude at the recording station is of the order of 10 to 20% of the direct SH wave amplitude depending on the attenuation model and assuming similar amplitudes radiated at the source. Deep moonquakes with best signal to noise ratio (mainly A01, A06 and A07 events) are assumed to be driven by tides [START_REF] Cheng | Tidal stresses in the moon[END_REF][START_REF] Araki | Focal processes of deep moonquakes[END_REF]. A recent analysis of the focal mechanism of these events, which is based on a model that fits only occurence times of the events, suggests fault plane orientations with a dip angle between 60 • and 70 • with large error bars [START_REF] Weber | Constraints on deep moonquake focal mechanisms through analyses of tidal stress[END_REF]. However, studies of these events based only on observations of P and S wave amplitudes and polarities favour focal mechanisms with a horizontal focal plane [START_REF] Nakamura | A1 moonquakes -Source distribution and mechanism[END_REF][START_REF] Koyama | Focal mechanism of deep moonquakes[END_REF]. As shown in figure 5b, these focal mechanisms have a maximum of SH energy radiated along the vertical direction, which correspond to the take-off angles of ScSH wave. In contrast, vertically polarized S waves will lose part of their energy when reflected at the core surface due to conversion into P waves reflected and transmitted at this interface. Smaller amplitudes for ScSV waves than for ScSH waves are thus expected, as observed on Earth.

Error propagation

The deep moonquake records present a low signal to noise ratio. Records allowing P and S arrival detections are obtained only after stacking many individual records from the same deep moonquake cluster. Owing to this low signal to noise ratio, to the low number of stations, and to the uncertainty on the radial seismic model, the error bar on the event location coordinates are usually quite large (see table 2). In order to minimize the effect of mislocation are used. These differential times present a low sensitivity to the structure below the station because S and ScS waves have similar crustal delays. In addition, these differential times also present a low sensitivity to event mislocations along latitude and longitude coordinates, even if, as we will see below, this sensitivity is non-zero. Because the largest error is on the quake depth, and because this parameter is strongly influencing the ScS-S differential time (hereafter called t d ), the event depth will be determined simultaneously with the core radius. However, ScS-S differential times alone do not allow us to relocate the lateral position of the event. Therefore, the stacking process can be effective only if the errors on latitude and longitude of the events generate relative ScS-S differential time variations between stations smaller than half the dominant period of the records, in order to ensure coherent (in phase) stack of the body waves. The differential time ScS-S at one station for an event location different from the theoretical one may be written:

t j d (θ 0 + δθ, φ 0 + δφ) = t j d (θ 0 , φ 0 ) + ∆t d + δt j d (8)
where (θ 0 , φ 0 ) are the (latitude,longitude) theoretical coordinates of the event,

∆t d = 1 N N j t j d (θ 0 + δθ, φ 0 + δφ) -t j d (θ 0 , φ 0 ) is the average differential time
shift residual over all the stations due to the event mislocation, and δt j d = t j d (θ 0 + δθ, φ 0 + δφ) -t j d (θ 0 , φ 0 ) -∆t d is the differential time shift residual of station j relative to this average. The second term on the right hand side produces a shift of all the ScS-S differential times that is displacing the ScS stack relative to its theoretical position, whereas the last term on the right hand side produces relative shift between the stations that may generate incoherent stacks if it is larger than half the dominant period of the body wave. The relative ScS-S differential time variations between stations i and j is defined by 2 presents

Rδt ij d = δt i d -δt j d . Table

Correction of the gain of horizontal sensors

The horizontal components of ground velocity are recorded by two different sensors (X and Y) of Apollo LP seismometers. Because these sensors may have slightly different gains, the rotation of X and Y components in order to reconstruct radial and transverse components may be biased. In order to correct for this effect, X/Y amplitude ratios were computed for each station in the coda of seismic events with high signal to noise ratios. Because of the high level of seismic scattering in the lunar crust, the seismic wavefield in the coda of lunar quakes is approximately in equipartion state [START_REF] Larose | Lunar subsurface investigated from correlation of seismic noise[END_REF][START_REF] Sens-Schönfelder | Temporal changes in the lunar soil from correlation of diffuse vibrations[END_REF]. Thus, if the gains of the two sensors were identical, and assuming only radial variations of seismic scattering properties, we should observe the same energy level on the two horizontal components [START_REF] Margerin | Energy partition of seismic coda waves in layered media: theory and application to Pinyon Flats Observatory[END_REF], and the X/Y ratio in the coda should be equal to one. If not, the X/Y ratio gives an estimate of the gain ratio of the two sensors. A limitation to this analysis is a small correlation between the horizontal components of the Apollo 12 seismometer discovered by [START_REF] Vinnik | First seismic receiver functions on the Moon[END_REF]. However, this study also indicates that this correlation is not observed on the other Apollo stations. S14, the codas of artificial impacts created by Lunar modules and stage IV of Saturn V rocket were used because of high signal to noise ratio for these events. For stations S15 and S16, deep moonquake codas with high signal to noise ratios were added to artifical impacts. The X/Y ratio of different events is almost constant for artificial impacts (before 1973), but it presents strong variations for deep moonquakes (after 1973) possibly due to both a lower signal to noise ratio and a less diffuse wavefield. However, the relative gain of X and Y components can be reasonably well estimated for all the stations. For stations S12 and S14, only artificial impacts are taken into account. The Y component is then corrected for this instrumental effect before performing the rotation in radial and transverse components for each deep moonquake event.

Correction of relative frequency responses of the stations

The horizontal component records of seismometers from the Apollo 14 and Apollo 16 missions present frequency contents significantly different from those of Apollo 12 and Apollo 15 missions. Figure 8a Figure 8: Logarithmic average of power spectral densities over all deep moonquake transverse records before (a) and after (b) correction of S14 and S16 transfer functions as a function of frequency for seismic stations S12 (in blue), S14 (in green), S15 (in black) and S16 (in red). On the left, curves for sensors X (dashed green) and Y (dotted green) of station S14 are also presented. All the curves have been smoothed by a moving average filter over 251 points in order to clarify the figures. Amplitude differences are due to differences of energy between the different deep moonquake records.

nal spectral peak at a frequency of 0.87 Hz, and Apollo 16 seismic station presents lower energy at low frequency compared to Apollo 12 and 15 stations. These anomalous features observed on the two horizontal sensors are not present on the vertical components of the records (not shown). It suggests that this effect is mainly due to the response of the ground just below the stations (site response). However, in order to be able to use the phase of the records in stacks, these records must have a similar frequency content.

The spectral amplitudes of S14 and S16 transverse components have thus been corrected by a spectral amplitude ratio computed between these stations and an average value between stations S12 and S15. The power spectral densities obtained after correction are presented in Figure 8b. These corrections of S14 and S16 records give a power spectral density similar to stations S12 and S15.

Data filtering

S waves usually have more energy at lower frequencies. Moreover, we expect both the crustal scattering and the attenuation to be reduced at low frequencies. Consequently, the Apollo data recorded in peaked mode were converted into the long period mode by using theoretical responses of these two acquisition modes, and filtered by an order 3 butterworth band pass filter with corner frequencies at 0.3 Hz and 0.9 Hz. This operation allows us to enlarge and equalize the available frequency band. The frequency limits have been chosen by trial and error in order to minimize the high frequency signal, and to reduce as much as possible the lower cut off frequency. Below 0.3 Hz and above 0.9 Hz, the noise dominates the signal in stacks of Apollo deep moonquake data. Figure 9 presents power spectral densities of the signals before and after the filtering process. The polarisation filtering described in [START_REF] Weber | Seismic Detection of the Lunar Core[END_REF] is a non linear method allowing to enhance linearly polarized seismic arrivals. Assuming V S V P ratio equal to 0.5 in the lunar regolith [START_REF] Gagnepain-Beyneix | A seismic model of the lunar mantle and constraints on temperature and mineralogy[END_REF], P and S body waves are linearly polarized for incidence angles smaller than 60 • = sin -1 V S V P [START_REF] Nuttli | The effect of the earth's surface on the S wave particle motion[END_REF]. Therefore, this method can be applied to most teleseismic P and SV wave arrivals on the Moon. However, assuming that SH waves are present only on the transverse component, such filtering is useless for these waves.

About polarisation filtering

The radial (R) and transverse (T ) components are computed from X and Y records by the following relation, assuming that X and Y sensors are aligned respectively along East and North directions (the real orientations of the sensors is given in Table 3):

R T = sin θ cos θ -cos θ sin θ × X Y ( 9
)
where θ is the azimuth of the event at the receiver. In the following we will note the real (true) displacements with superscripts T and the measured displacements without superscript, and assume that Y = Y T and X = γX T .

The gain corrections for the horizontal components has been presented above, and γ is equal to the X/Y ratio determined previously. If this correction is not applied before the rotation of X and Y components, the relations between measured radial and transverse components and the true ones is the following:

R T = γ sin θ cos θ -γ cos θ sin θ × sin θ -cos θ cos θ sin θ × R T T T (10)
In particular, the radial component can be expressed as R = a 11 (γ)R T + a 12 (γ)T T with a ij (γ) being the terms of the matrix A corresponding to the matrix product described above.

The polarisation filter [START_REF] Weber | Seismic Detection of the Lunar Core[END_REF] can be written as OZ j = Z j M j and OR j = Z j M j respectively along the vertical and radial components, with

M j = n i=-n Z j+i R j+i ,
where j is the time step and n = 6. Now assuming that SV body wave has a motion W S j and an incidence angle j o , and that ǫ z and ǫ r are respectively the noise on the vertical and radial component, the true vertical and radial components can be written as Z T j = W S j sin(j o ) + ǫ z and R T j = W S j cos(j o ) + ǫ r . Reporting these equations in (10), and assuming Z = Z T , we obtain

Z j R j = (W S j sin(j o ) + ǫ z ) a 11 (γ)W S j cos(j o ) + a 11 (γ)ǫ r + a 12 (γ)T T (11)
Consequently, the non linear polarisation filter is subject to two different Because the polarisation filtering is able to enhance the coherent SV wave signal on both radial and vertical components, we can use this information in order to refine the arrival time picks of SV waves and to infer the length of the SV pulse. Figure 11 presents a focus on the polarisation filtering of radial components around the SV arrival. The SV wave arrivals are clearly enhanced by this process on both the raw and filtered records. More interestingly, the SV pulse length can be estimated in the 5 to 10 seconds range, and it looks to be varying more with event number than with station number.

This observation suggests that the pulse length reflects more the properties of the signal coming from the source region than the variation of scattering properties below the stations.

However, P and S waves reflected back from discontinuities at radius smaller than 500 km inside the Moon have incidence angles smaller than 2 • . Therfore, the amplitude of the SV wave on the vertical component is predicted to be at best 4% of the one on the radial component. Due to the high coda energy and to the low signal to noise ratio of deep moonquake stacks, the SV wave signal coming from deep reflectors inside the Moon will be smaller than the noise on the vertical component (W S j sin(j o ) << ǫ z ). The same argument also holds for deep reflections arriving as P waves at the station. In this case, the polarisation filtering will only enhance time intervals for which the noise (or coda signal) is in phase on both components.

In conclusion, the polarisation filtering is useless for ScSH wave detection.

It gives good results for SV and P waves when the incidence angle is not too small, and allows us to estimate the length of the S pulse in the 5-10 s range. However, it cannot be used for the detection of waves reflected by deep discontinuities inside the Moon because in this case the noise dominates the signal on one component, and both OZ and OR filtered signals present enhancements for time intervals for which Z and R codas are in phase. For all these reasons, we decided not to use this method for ScSH wave detection. Figure 12: Results of the SH wave optimal waveform and delay times inversion. From bottom to top, optimal SH waveform (thick line) and filtered records of stations S12, S14, S15 and S16 after optimal alignment. From left to right, records for events A01, A06 and A07. Vertical dashed lines delimit the time interval with maximum coherent energy which is used below to define the optimal SH waveform.

S wave alignment and differential travel time computations

The arrival times of S waves on the transverse components of deep moonquake stacks are hand picked by visual inspection of vertical, radial and transverse components. During this process, the polarity of some records has been changed in order to ensure similar polarities of S waveforms for the different records of the same deep moonquake. The stacks of deep moonquake records strongly depend on the S wave arrival time. In order to obtain a very precise alignment of SH wave between the different stations, the records are inverted by a non-linear method searching for optimal SH waveform and relative delay times between the stations [START_REF] Chevrot | Optimal waveform and delay time analysis by simulated annealing[END_REF]. Results of this waveform matching procedure are presented in figure 12. Interestingly, the best matching waveform presents significant energy on a length that is consistent with the SV pulse lengths obtained by polarisation filtering (≈5 s for A01, ≈7 s for A06, ≈10 s for A07). The SH hand picks are corrected by delays deduced from this analysis.

The arrival times of S and ScS waves are computed inside the best radial models for each core radius using the deep moonquake relocations obtained for these models. These computations are performed by a new implementation of the Tau-P algorithm [START_REF] Buland | The computation of seismic travel times[END_REF][START_REF] Calvet | Traveltime sensitivity kernels for PKP phases in the mantle[END_REF].

Stacking process

All the waveforms are aligned on the S wave arrivals and normalized to maximum amplitude of the first one hundred seconds of the coda. Again, events with less than four recording stations and events for which horizontal location errors produce relative differential times between stations larger than half the dominant period (≈ 1.25 s) are excluded. Only high quality events A01, A06 and A07 are kept. The stacks of the waveforms for each event are performed by computing the differential times ScS-S for a given moonquake depth D i and a given core radius R core using the corresponding best radial model. After alignement and stack, the ScSH waveform energy is computed in a 10 second window after the predicted arrival time on the stacked trace. This process is repeated for all the best radial models obtained previously corresponding to core radii between 250 km and 490 km with a 5 km step, and for event depth in the ±400 km range around the epicenter location with a 2 km step. Such a process allows us to compute, for each event i, NRJ ScS i (R core , D i ) the energy of stacked ScS waveforms depending on core radius R core and event depth D i . The semblance of the stacked records is also defined as the ratio between the stacked ScS waveform energy and the average energy of the records used to perform the stack [START_REF] Neidell | Semblance and other coherency measures for multichannel data[END_REF]. This parameter is varying between zero (complete destructive interference) and one (fully constructive interference), providing a measure of stack efficiency normalized to the energy of the records. As for the energy, the semblance depends on core radius R core and event depth D i :

SEMB ScS i (R core , D i )
In order to include additional a priori information on the events depth, stacked energies and semblances are multiplied by a gaussian function of the form G i (D i ) = exp(-

(D i -D 0 i ) 2 σ 2 i ) with D 0
i the a priori event depth and ∆D 0 i its error. The standard deviation of the gaussian is defined by σ i = 10∆D 0 i , in order to allow large variations around the a priori event depth. A priori deep moonquake locations and related errors on the parameters are extracted from the study by [START_REF] Gagnepain-Beyneix | A seismic model of the lunar mantle and constraints on temperature and mineralogy[END_REF]. This weigth function

G i (D i )
is only a smooth way to impose event depth in the ±5∆D 0 i range around theoretical value.

Figure 13 shows the ScSH stack energies obtained for the three events A01, A06 and A07. On these plots, a correlation between core radius and event depth is clearly seen. However, the ScSH energy presents a peak in the 380-400 km core radius range for the three events.

Once these energies are computed, their maximum value for each core radius are defined by

NRJmax ScS i (R core ) = max D i NRJ ScS i (R core , D i ) .
These functions are plotted in figure 14(a-c) for the three events selected. They present a global or local maximum in the 380-400 km core radius range.

However, the statistical test performed below demonstrate that, when taken separately, the maxima of these events are statistically significant only for events A06 and A07.

To investigate the statistical significance of the results, we consider the average of maximum ScS stack energies of individual events:

NRJsum(R core ) = 1 N ev Nev i NRJmax ScS i (R core ) ( 12
)
with N ev the number of events. However, even if the records are scaled before the stacking process, NRJmax ScS i (R core ) functions present different amplitudes for the different events. In our case, NRJsum(R core ) is dominated by events A01 and A06. In order to give an equal weight to all the events, these functions are scaled to their standard deviation and their mean is substracted before summation in order to produce energy plots including all the events with a similar contribution. The final result is a scaled energy depending only on core radius which is defined for ScS waveform stacks by:

NRJ sum(R core ) = 1 N ev Nev i=1 (NRJmax ScS i (R core )-< NRJmax ScS i (R core ) > std(NRJmax ScS i (R core ))
(13) with <> and std() respectively the average and standard deviation operators.

The same operations are performed for the semblance by computing for each event the maximum semblance at different core radius (SEMBmax ScS i (R core ))

and the average of these curves over all the events (SEMBsum(R core )). The scaled version of semblance curve does not have to be computed because the semblance measure is already normalized to the energy of the records.

Functions NRJsum(R core ), NRJsum(R core ) and SEMBsum(R core ) are presented in panels (d), (e) and (f) of figure 14. This figure demonstrates that when the three events are taken together the energy maximas at 380 and 395 km core radius are statistically significant (at three standard deviations level).

In order to test the dependence of the results on the filtering process, the whole process (including S wave alignment) is repeated on the data before acquisition mode conversion and bandpass filtering. The results are presented in panels (g), (h) and (i) of figure 14. The peaks of ScSH wave energy around 380 and 395 km core radius are also present but less statistically significant.

The consistency of the results between the two different frequency ranges strongly argue in favour of a ScSH detection. Moreover, the lower ScSH energy at high frequencies (in peaked mode) than at low frequencies (in long period bandpass mode) justifies the search for this body wave at frequencies as low as possible.

Because the ScSH wave stacks may be slightly shifted inside the ten seconds time window, the core radius obtained may be also shifted. For example, if the energy of ScSH wave stack inside the window is significant only 2 or 3 s after the beginning of the window, S-ScS travel time is under-estimated and core radius over-estimated. In order to remove this effect, the records obtained after filtering have been deconvoled from the optimal source time functions presented on figure 12. After deconvolution, SH energy is focused on a 2 s window. The stacking process is reproduced by computing energy on a 2 s window centered on the ScS arrival. The results are presented in figure 15, for the three different moonquakes. Clear maximas of ScSH energy are obtained for events A06 and A07, but at slightly smaller core radius (370-390 km), which suggests that the over-estimation of the core radius due to the 10 s window may be real. However, the small size of stacking window has the drawback to focus the stack energy in a very narrow core radius range.

Due to quake mislocations and also possibly to lateral heterogeneities at the base of the lunar mantle, the real ScS-S differential travel times vary from one event to the other and different core radius are obtained for the different events. As a consequence, the three curves do not interfere constructively when summed together.

The test presented above and the increasing difficulty to explain geodesic observations with increasing core radius above 380 km, favour a core radius value at the lowest range of the broad energy peak observed on figure 14.

Optimal core radius estimate is 380 km, with a 30 km error bar estimated from the range of core radius for which the energy is above the two standard deviation level.

A similar exercise is also done for the radial components of the records in order to detect ScSV waves. The results obtained on filtered data are presented on figure 16. As expected, ScSV wave is not detected inside the S coda because its amplitude is smaller than ScSH one. This exercise also demonstrates that the bootstrap method described below is pertinent for the validation of ScSH detections. Figure 14: Plots of stacked maximum energies N RJmax i (R core ) for events A01 (a), A06 (b), A07 (c), and the average of these energies N RJsum(R core ) (d), the scaled average of these energies N RJsum(R core ) (e) (red curves), and the average semblance SEM Bsum(R core ) (f) (red curve). Panels (g), (h) and (i) are identical respectively to panels (d), (e) and (f), but for data before acquisition mode conversion and bandpass filtering. The average value (blue plain line), the one standard deviation around this value (blue dashed lines) and the two standard deviation value (blue dot-dashed line) obtained for the bootstrap ensemble are shown on each plot. On panels (d), (e), (f), (g), (h) and (i) curves for individual events (thin black lines) are also presented.

Validation of the results with bootstrap

We applied a bootstrap method to test the statistical significance of our results [START_REF] Efron | An introduction to the bootstrap[END_REF]. The SH arrival picks are perturbed randomly in a ±30s range. We repeat the stacking process for each new data set. We then compute the mean, maximum and standard deviation over 150 random realizations of the resulting stacks of energy. This method gives an estimate of the background noise corresponding to random waveform alignment. A peak can be considered significant if it is above two standard deviations of the bootstrap ensemble.

As described previously, the two standard deviation level is passed in the 380-400 km core radius range for filtered data using a 10s window (Figure 14), and in the 370-390 km core radius range after deconvolution of optimal SH waveforms for events A06 and A07. The statistical validation is obtained for ScSH detection and not for ScSV detection. Moreover, an error estimate of ±40 km is provided. This error will be evaluated more precisely below. This shift correspond to the value of ∆t d due to the error on the (longitude, latitude) coordinates of the quake, and it should be below 2 s for the events processed in this study. A good correlation between ScSH and SH stacks is obtained for events A06 and A07, but it is low for event A01. This waveform similarity, already enhanced by the success of stacking after deconvolution of SH optimal waveform, gives further support for a ScSH detection.

Stacked waveforms comparison

Discussion

Error estimates

Figure 14 provides a formal error estimate on the core radius which takes into account the uncertainty in the stacking process, and errors on the depth of the events. However, the error on the seismic model must also be taken into account in order to estimate its effect on the core radius. From the set of models generated by the NA algorithm at 380 km radius, we compute the maximum and minimum values of seismic velocities and densities at each depth inside the ensemble of models within the half width of the probability peak. These values define the error bars on the best radial model. They are plotted on figure 19a. These error bars are almost constant inside the mantle for density and P-wave velocity with a relative error of respectively 0.2% and 3.3%. However, the error on S-wave velocity presents a minimum in the mid-mantle (3%) and a maximum (11%) at the base of the mantle. The strong increase in S-wave velocity errors at the base of the mantle translates to a significant error on the absolute value of core radius. Considering a vertically incident ScS wave from a deep moonquake at a depth of 900 km, the integrated error on the S-wave velocity profile below the event produces an error on the ScS travel time of about 9% on the two-way travel time from the event to the core, or about 19 s. This travel time error can be translated directly into an error on the core radius of approximately 33 km. This rough estimate is probably an upper estimate because the vertical incidence of ScS rays has the highest sensitivity to core radius. However, it highlights the difficulty to obtain an absolute value of core radius without S-wave arrival time measurements at large epicentral distances. Taking into account the true ray geometries, and the bootstrap analysis of ScSH wave detection, the error bar on the core radius is about 40 km.

The error on the average core density for a 380 km core radius is about 0.3 kg/cm 3 . By looking at figure 4c, varying the core radius in a range of ±40 km around its best estimate gives approximately two times larger average core density variation than the error bar estimate at 380 km core radius.

Moreover, because the core density is inferred from mass and moment of inertia, and due to the low value of the core moment of inertia (of the same order as the error bar on the moment of inertia of the whole planet), small variations of crustal thickness and density may generate large variations of core density. The average crustal thickness is not inverted here, but the error on this parameter requires to increase our estimate of the error on the core density. Consequently, a conservative estimate of average core density is 5.2±1.0 kg/cm 3 .

Comparison with previous results

Our estimates of core radius at 380±40 km and average core density at 5.2±1.0 kg/cm 3 are respectively at the upper and lower limits of previous estimates from both lunar laser ranging observations [START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF] and induced magnetic moment of the Moon [START_REF] Hood | Initial measurements of the lunar induced magnetic dipole moment using Lunar Prospector magnetometer data[END_REF]. This seismic estimate of core radius is significantly larger than the 330 km proposed by [START_REF] Weber | Seismic Detection of the Lunar Core[END_REF], which in some way indicates a large sensitivity to uncertainties in the seismic model at the base of the mantle and to error propagation. However their large error bars of these core radius estimates cover the domain of previously inferred values. Moreover, the detection of ScSH waves and not ScSV waves supports a liquid core beneath the core mantle boundary, as already suggested by the dissipation of lunar rotation [START_REF] Williams | Lunar rotational dissipation in solid body and molten core[END_REF]. The average core density is relatively low and suggests either a high level of light elements in a completely fluid core, or a large core temperature. For example, if the average core density is explained by the presence of sulfur alone, more than 10% are required for the nominal value [START_REF] Sanloup | Density measurements of liquid Fe-S alloys at high-pressure[END_REF][START_REF] Balog | Equation of state of liquid Fe-10 wt % S: Implications for the metallic cores of planetary bodies[END_REF]. But, when the error bar is taken into account the percent of sulfur content can vary from almost zero to large values. However, because an inner core is expected for thermodynamical reasons [START_REF] Wieczorek | The constitution and structure of the lunar interior[END_REF] and its size is not constrained by this study, the interpretation of the average core density into core composition is strongly limited. 5 gives station corrections. These corrections are small and their relative error is large. [START_REF] Nakamura | Seismic Q of the lunar upper mantle[END_REF]Nakamura et al., 1982) and arbitrarily fixed in the deep mantle. Seismic velocities and attenuation parameters inside the core are arbitrarily fixed by the authors. However, the internal structure of the lunar core remains largely unknown, and the model is still characterized by rather strong uncertainties on the different parameters owing to the paucity of available data.

The constraints on the core density are strongly related to geodesic parameters such as the polar moment of inertia of the Moon, but also to the average crustal structure. The average crustal thickness is fixed in our model, but its variation can change significantly the polar moment of inertia budget of the planet, and consequently the core density obtained. Futher constraints brought by SELENE and GRAIL missions on these parameters will indirectly strongly constrain the average core density, and the radius of the inner core which is expected for thermodynamical reasons [START_REF] Wieczorek | The constitution and structure of the lunar interior[END_REF], and for which a seismic signature has been suggested by [START_REF] Weber | Seismic Detection of the Lunar Core[END_REF].

Our model relies on the hypotheses of homogeneity and adiabaticity of the lunar mantle. Consequently, radial and lateral deviations from this model, that may eventually be detected by future lunar seismometers, will provide additional constraints on the internal dynamics of the lunar interior. In particular, a high temperature gradient or partial melt inside the upper and lower boundary layers of the mantle may strongly modify seismic velocities and density in these regions. Both the very broad band planetary seismometers developed in the last decades [START_REF] Lognonné | Ultra broad band seismology on InterMarsNet[END_REF][START_REF] Lognonné | The NetLander very broad band seismometer[END_REF][START_REF] Lognonné | Planetary Seismology[END_REF] and future lunar geophysical stations, such as SELENE2 [START_REF] Tanaka | The science objectives of the SELENE-II mission as the post SELENE mission[END_REF] and Lunette/ILN [START_REF] Neal | LUNETTE: Establishing a Lunar Geophysical Network Without Nuclear Power Through a Discovery-Class Mission[END_REF] will provide crucial additional data to further improve our knowledge of the lunar interior.

Figure 1 :

 1 Figure 1: A priori crust model extracted from Gagnepain-Beyneix et al. (2006) assuming a constant density inside the crust below the 1 km thick regolith layer.

Figure 2 :

 2 Figure 2: Example of parameter space sampling by NA algorithm for a core radius of 380 km. (a) sampling projected into the plane of a and b Birch law parameters. (b) sampling projected into the plane of A and B parameters of VP VS ratio law. (c) sampling projected into the plan of α and ρ c parameters. Each dot represents a model vector and best model values are indicated by a red star.

Figure 3 :

 3 Figure 3: (a) Minimum values of cost function J(m) (in red), χ 2 geod (in black) and χ 2 seismo (in blue) obtained for each core radius. Values of χ 2 geod (black dashed line) and χ 2 seismo (blue dashed line) corresponding to minimum values of the cost function J(m) are also plotted. (b) P-wave (in red) and S-wave (in blue) velocities (in km/s), density (in green) and VP VS ratio (light blue) of the best models obtained for each core radius (dashed lines). The model of Gagnepain-Beyneix et al. (2006) is plotted as plain black lines for comparison.

  errors and seismic structure just below the stations, ScS-

Figure 4 :Figure 5 :

 45 Figure 4: Correlations between parameters when considering the best radial models. (a) best values of a and b Birch law parameters. (b) best values of A and B parameters of VP VS ratio law. (c) crust and core densities as a function of core radius (in km) for the best radial models. Color scale of circles goes from black to white with increasing core radius parameter.

  an estimate of the maximum values (Max(Rδt ij d )) and standard deviations (Std(Rδt ij d )) of relative ScS-S differential time variations between stations inside the one standard deviation error ellipse of event positions. These values are controlled by the error on the event location, but also by the event position relative to the stations.When the back azimuth of the event is approximately the same for all the stations (event outside the network area), ∆t d may be large, but Rδt ij d values remain small. In contrast, for an event located in the middle of the network, ∆t d is small, but Rδt ij d values present large variations.

328 4 .

 4 Data and methods329This section describes the data and method used to detect SH waves 330 reflected from the lunar core surface. First, the stacking process of individual 331 deep moonquake records along X and Y horizontal sensors is described. Then, 332 a new method, based on the energy equipartition inside the coda waves, is 333 applied to correct for relative gain variation between X and Y horizontal 334 sensors. In addition the site response below the station is corrected, and data are filtered in the frequency range with highest signal to noise ratio.336Finally, the stacking process is described, and bootstrap validations of the 337 results are presented.4.1. Stacks of deep moonquake individual recordsIndividual records of deep moonquakes are first aligned by cross-correlation of the vertical component and the horizontal component with the best S/N ratio. The part of the individual records presenting spikes have been removed by hand from the stacks. As a result, the number of individual records stacked varies with time along the stack. However this method is probably the best one to ensure that no biases are introduced by a spike correction method.The number of individual records stacked depends on the seismic station and on the deep moonquake nest. This number is usually comprised between 20 and 100. An example of deep moonquake stack is presented in figure6which shows the stacks of deep moonquake cluster A01 on the vertical component of Apollo 12 station.

  Figure 6: (a) examples of individual records of events belonging to deep moonquake cluster A01 and recorded on the vertical component of station Apollo 12 after alignement by crosscorrelation. (b) evolution of the stack as a function of the number of records stacked. From top to bottom each line includes a new record in the stack. The numbers on the left of each trace give the maximum amplitude of the trace expressed in Apollo digital units.

  presents a logarithmic average of the power spectral densities of all the deep moonquake signals selected per seismic station. Apollo 14 seismic station presents an addition-

Figure 9 :

 9 Figure9: Logarithmic average of power spectral densities over all deep moonquake transverse records before (dotted lines) and after (plain lines) modes conversion and filtering as a function of frequency for seismic stations S12 (in blue), S14 (in green), S15 (in black) and S16 (in red). All the curves have been smoothed by a moving average filter over 251 points in order to clarify the figures.

  [START_REF] Weber | Seismic Detection of the Lunar Core[END_REF] used polarisation filtering in order to detect body wave arrivals reflected by deep Moon discontinuities. We evaluate the performance of this filter for picking S waves and detecting reflections from deep Moon discontinuities.

  noise sources. First, the term a 12 (γ)T T perturbs the radial component. Second, the SV wave signal should be over the noise level along the vertical component (W S j sin(j o ) >> ǫ z ) in order to ensure enhancement of the SV wave signal. Our tests demonstrate that the first noise source due to wrong rotation in the horizontal plane is generally negligible because γ is close to one. Figure 10 presents the polarisation filtering performed on raw deep moonquake stacks after rotation for event A01 and station S15. This figure can be compared directly to figure S1 of Weber et al. (2011). SV wave signal is clearly enhanced by the polarisation filtering both on OZ and OR because the SV wave signal is clearly dominating noise on both vertical and radial components. Such a clear enhancement of SV pulse by polarisation filtering is not observed in Weber et al. (2011) study. This comparison suggests that the stacks of individual deep moonquake signals coming from the same deep moonquake cluster (A01 in this case) was performed differently between the two studies.

  Figure 10: (top three traces): Raw deep moonquake stacks components (cluster A1) at Apollo station 15. (bottom three traces): Polarization filtered components of the same stack. S wave arrival time predicted by our best model (blue line) is presented.

Figure 11 :

 11 Figure11: Polarisation filtering of radial components (OR j ) aligned on the hand picked SV wave arrival times for raw records (dashed lines) and records after filtering (plain lines). From top to bottom, records of stations S12, S14, S15 and S16. From left to right, records for events A01, A06 and A07. The vertical component of station S14 for event A07 is not available in our data base, consequently the polarisation filtering cannot be performed. The vertical dotted lines indicate the S wave arrival time predicted by our best moon model.

Figure 13 :

 13 Figure 13: Stacked energies N RJ ScS i (R core , D i ) in color, as a function of core radius (in km) and event depth (in km), for events A01 (a), A06 (b) and A07 (c). The energies are multiplied by a gaussian function G i (D i ) centered around a priori depth D 0 i . Color goes from blue (low energy) to red (high energy), and the maximum value is indicated by a black diamond.

Figure 15 :

 15 Figure15: Plots of stacked maximum energies N RJmax i (R core ) after deconvolution of optimal SH waveforms for events A01 (a), A06 (b), A07 (c). The average value (blue plain line), the one standard deviation around this value (blue dashed lines) and the two standard deviation value (blue dot-dashed line) obtained for the bootstrap ensemble are shown on each plot.

Figure 16 :

 16 Figure 16: Plots of N RJsum(R core ) (a), N RJsum(R core ) (b) and SEM Bsum(R core ) (c) (red curves) for stacks of the radial components of the records. The average value (blue plain line), the one standard deviation around this value (blue dashed lines) and the two standard deviation value (blue dot-dashed line) obtained for the bootstrap ensemble are shown. Curves obtained for individual events (thin black lines) are also presented.

Figure 17

 17 Figure 17 presents record sections of transverse components aligned on predicted ScSH arrival for a 380 km core radius. After alignment on ScSH, the waveforms are stacked in order to produce a ScSH stacked waveform for each event. Similar processing is also done for the SH wave. The comparison of these stacked waveforms is shown in figure 18. In this figure, the SH and ScSH waveforms have been cross-correlated on a 10 s window, and are presented for the best correlation within a ±2 s time shift of ScSH stack.

  Figure 17: Record sections of transverse components plotted as a function of epicentral distance (in degrees) and time relative to predicted ScSH arrival time for a 365 km core radius (in seconds). Red marker indicate the SH wave arrival time pick. Vertical dashed lines indicate the 10 second window used to stack ScSH waveforms. Different colors are used for the different quakes: A01 in black, A06 in blue and A07 in green.

Figure 18 :

 18 Figure 18: Comparison between SH (black plain line) and ScSH (blue dashed line) stacks as a function of time and event number. The ScS stacks are aligned for best correlation in the ±2 seconds range. Blue diamonds indicate the start of ScS stack window before alignment. Events numbers are 1 for A01, 2 for A06 and 3 for A07. Correlation coefficient over the 10 seconds window is indicated along the stack before the stack window.

6.=

  Description of the best model This section gives seismic event location parameters, station corrections, and the reference model. The reference model is denominated VPREMOON for Very Preliminary REference MOON model. The parameter values, and related error bars at fixed core radius, for this model are: ρ c = 2.762±0.0048, a = -8.0783 ± 0.790, b = 4.728376 ± 0.245, A = 1.816595 ± 0.313, B = -0.000054 ± 0.000251km -1 and R core = 380 km. In this model, calculated values of geodesic observations are k calc 2 = 0.0223, h calc 2 = 0.0394, IR calc = 0.3932 and l calc 2 = 0.0106 (this last parameter is not inverted here). The values obtained at best fit are χ 2 seismo = 1.4814 and χ 2 geod = 0.0575. For comparison, our starting model derived from Gagnepain-Beyneix et al. (2006) model gives k calc 2 = 0.0231, h calc 2 = 0.0408, IR calc = 0.3931, l calc 2 the relocations of the events used in this study inside the VPREMOON model. These locations are close to the locations published by Gagnepain-Beyneix et al. (2006) because the seismic model is very similar in most of the mantle. The errors on these locations are underestimated, because the uncertainty on the velocity model is not formally taken into account in their computation. Table

Figure 19 :

 19 Figure19: VPREMOON seismological (a) and geodesic models (b). On the left (a), P and S wave velocities (in km/s) and density (in kg/cm 3 ) are plotted as a function of radius (in km) for the seismological model. The error bars are indicated for each parameters by dot-dashed lines. On the right (b), Pressure (in GPa), gravity (in 0.1m/s 2 ) and density (in kg/cm 3 ) are plotted as a function of radius (in km) for the geodesic model.

  a preliminary reference model of the Moon based on a priori crustal structure, physical constraints on density and seismic velocities variations with depth, by fitting both seismological and geodesic data. The core radius is determined from the detection of transversely polarized core reflected S wave from a low number of deep moonquakes properly located and presenting high signal to noise ratio. The VPREMOON model constrains the core size to 380±40 km radius and average density to 5.2±1.0 kg/cm 3 , and favours a liquid outer core. It constitutes the first reference model including simultaneously physical constraints, geodesic and seismological observations, and detection of S waves reflected on the lunar core.

Table 1 :

 1 Summary of parameter ranges explored in the construction of physical Moon models.

  IR obs , IR calc , σ h 2 , σ k 2 and σ IR respectively observed and calculated values of h 2 and k 2 Love numbers, polar moment of inertia ratio and corresponding

	equal weight to seismological and geodetic observations in order to construct
	a reference model constraining both density and seismic velocities.
	The h calc 2 , k calc 2	and IR calc values are computed during model construction.
	with t obs i S phases), t calc the observed travel time of the seismic phase (including both P and i the predicted travel time inside model m, σ 2 i the travel time error, Ns the number of travel time measurements; and h obs 2 , h calc 2 , k obs 2 , k calc 2 ,
	errors. The first part of the right hand side of equation (7) will be referred
	to as χ 2 seismo , whereas the second part as χ 2 geod . This cost function gives an

Table 2 :

 2 Summary of data used in this study. Parameters are the following: latitude, longitude and depth of the quake, and corresponding errors[START_REF] Gagnepain-Beyneix | A seismic model of the lunar mantle and constraints on temperature and mineralogy[END_REF], initial number of stations available, number of stations selected after removing bad quality records, parameters Max(Rδt ij d ) and Std(Rδt ij d ).

	Event Name	Lat. (deg)	Long. (deg)	Depth (km)	∆lat. (deg)	∆lon. (deg)	∆depth (km)	Init. # stat.	# stat. sel.	Max(Rδt (s)	ij d )	Std(Rδt (s)	ij d )
	A01	-17,44	-38,37	917	1,1	0,6	11	4	4	1,82		1,1
	A05	-31,04	-44,92	902	3	3	80	2	2	1,34		-
	A06	49,7	54,69	860	1	0,7	11	4	4	1,57		0,9
	A07	23,97	53,7	900	0,8	0,7	12	4	4	1		0,94
	A08	-27,97	-28,08	940	2	1,2	21	1	1	-		-
	A09	-37,8	-30,85	975	4	2,6	43	3	2	8,56		-
	A10	34,01	-28,04	1139	3	3	80	2	2	0,57		-
	A11	9,6	19,49	1233	0,8	0,7	12	3	2	0,73		-
	A14	-28,7	-33,94	880	1,7	1,3	22	3	2	4,34		-
	A15	-0,94	-2,96	885	3	3	80	3	3	13,39		8,8
	A16	6,79	5,14	1105	1,2	0,7	18	3	3	2,26		1,55
	A17	23,08	-17,97	861	2,6	0,4	15	3	3	7,51		5,27
	A18	18,56	34,72	882	1,8	0,9	24	3	2	3,4		-
	A19	15,96	37,91	841	3	3	80	3	2	9,46		-
	A20	21,72	-41,01	1055	2,4	0,8	13	3	3	4,32		2,83
	A21	-15,8	-43,49	1060	3	3	80	2	1	-		-
	A24	-36,85	-38,9	980	2,1	1,7	32	4	3	4,91		2,49
	A25	34,4	59,3	898	2	1,7	26	4	4	2,57		2,1
	A26	12,19	10,22	1135	1,5	0,7	20	3	2	1,5		-
	A27	22,51	18,53	1059	1,9	1,4	14	4	2	3,49		-
	A28	32,69	11,64	720	3	3	80	3	1	-		-
	A30	11,78	-34,26	921	1,5	1	23	4	4	1,87		2,03
	A33	6,89	117,75	887	1,5	1,3	30	3	2	0,46		-
	A34	7,04	-9,28	932	1,2	0,6	26	3	3	3,8		2,58
	A41	13,9	-26,79	953	5,6	2,2	84	3	1	-		0
	A42	22,68	-53,46	1004	1,8	1,4	24	3	2	3,94		-
	A44	51,85	57,08	956	5,8	1,9	20	3	3	1,98		0,78

Owing to the instrument response of Apollo stations, the largest dominant 322 period ensuring a good signal to noise ratio is about 2.5 s. Therefore, the only 323 events for which the stacking process is expected to work properly according 324 to the error on lateral event coordinates are events A01, A06, A07 and A44.

325

In addition, detection by stacking process is more reliable if 4 stations are 326 available. Consequently, only deep moonquake stacks for events A01, A06 327 and A07 can possibly allow ScSH reliable detections.

Table 3

 3

		:	Azimuth	(in	degree)	of	X	and	Y	sensors	extracted
	from	Apollo	Scientific		Experiments		Data	Handbook	available	at
	http : //nssdc.gsfc.nasa.gov/planetary/online books.html
	Station	X		Y
	name	azimuth (in deg.) azimuth (in deg.)
	S12		180		270
	S14		0		90
	S15		0		90
	S16		334.5		64.5

Table 6

 6

	and figure 19 present the

Table 5 :

 5 P and S wave station corrections and associated error bars.

	Station P cor.	P cor.	S cor.	S cor.
	name	(in s) error (in s) (in s) error (in s)
	S12	-1.034	0.250	-0.229	0.462
	S14	1.162	0.230	-0.584	0.413
	S15	-0.298	0.279	-0.950	0.468
	S16	0.171	0.758	1.763	1.344

Table 6 :

 6 VPREMOON model: on the left the seismic model and on the right the geodesic model in which Moho depth as been corrected to 40 km. Attenuation parameters are extracted from previous studies

acknowledgments

We acknowledge Vadim Monteiller for helpful discussions concerning quake relocation methods. This study was supported by Centre National d'Etudes Spatiales (CNES) research projects, Programme National de Planétologie (PNP) of CNRS-INSU, PRES "Université de Toulouse" (Toulouse University) and Campus Spatial Paris Diderot (IPGP contribution number XXX).

We thank the two anonymous reviewers for improving the paper by their constructive comments.


Selda C ¸alkavur
email: selda.calkavur@kocaeli.edu.tr

Patrick Solé
email: sole@enst.fr

Vector space multisecret-sharing scheme based on Blakley's method

Keywords: Secret sharing, multisecret-sharing scheme, vector space. 2010 AMS Mathematics Subject Classifications: 94A62, 15A03, 15A18

Secret sharing schemes were introduced by Adi Shamir and George Blakley, independently, in 1979. In a (k, n)threshold secret sharing scheme, any set of at least k out of n participants can retrieve the secret but no set of (k -1) or less can. Shamir's secret sharing scheme is more popular than Blakley's, even though the former is more complex than the latter. The practical reason is that Blakley's scheme lacks determined, general and suitable matrices. In this paper, we present a multisecret-sharing scheme based on vector spaces over R n . We use Blakley's method. This scheme is ideal in the sense that the size of each secret equals the size of any share.

I. INTRODUCTION

The secure storage of the private keys of a cryptosystem is an important problem in cryptography and security. The possession of a highly sensitive key by an individual may not be desirable as the key can easily be lost or as the individual may not be fully trusted. Giving copies of the key to more than one individual increases the risk of compromise. A solution to this problem is to give shares of the key to several individuals, forcing them to cooperate to find the secret key. This not only reduces the risk of losing the key but also makes compromising the key more difficult. In threshold cryptography, secret sharing deals with this problem. The secret is shared among a group of n users so that only when a sufficient number k of them come together can the secret be reconstructed. In the literature well-known secret sharing schemes include Shamir [START_REF]How to share a secret[END_REF] based on polynomial interpolation, Blakley [START_REF] Blakley | Safeguarding Cryptographic Keys[END_REF] based on hyperplane geometry, and Asmuth-Bloom [START_REF] Asmuth | A modular approach to key safeguarding[END_REF] based on the Chinese Remainder Theorem. These schemes are called (k, n)-threshold secret sharing.

In this paper we use Blakley's scheme. Blakley's method [START_REF] Blakley | Safeguarding Cryptographic Keys[END_REF] uses principles of geometry to share the secret. In [START_REF] Blakley | Linear algebra approach to secret sharing schemes[END_REF], Blakley et al. only provided a guideline on how to design a matrix of linear systems for perfect secrecy, and no actual matrix was given. Recently, researchers began to use Blakley's geometry-based secret sharing approach in the area of secret image sharing [START_REF] Chen | A geometry-based secret image sharing approach[END_REF], [START_REF] Tso | Sharing secret images using Blakley's concept[END_REF]. Chen et al. [START_REF] Chen | A geometry-based secret image sharing approach[END_REF] and Tso [START_REF] Tso | Sharing secret images using Blakley's concept[END_REF] independently applied Blakley's scheme for secret image sharing.

Ulutas ¸et al. [START_REF] Ulutas | Improvements in Geometry-Based Secret Image Sharing Approach with Steganography[END_REF] proposed an enhanced scheme for secret image sharing, which adopts Blakley's secret sharing method and Steganography together to share the secret and create meaningful shares. As for threshold cryptography, Bozkurt et al. [START_REF] Bozkurt | Threshold Cryptography Based on Blakley Secret Sharing[END_REF] proposed the first threshold RSA signature scheme as the underlying secret sharing scheme. Hei, Du and Song [START_REF] Hei | Two Matrices for Blakley's Secret Sharing Scheme[END_REF] present two matrices that can be used for Blakley's secret sharing scheme.

Multisecret-sharing schemes form another family of secret sharing schemes. This scheme was proposed [START_REF] Horn | Comment: "Multistage secret sharing based on one-way function[END_REF], [START_REF] He | Multistage secret sharing based on one-way function[END_REF], [START_REF] Li | A New (t, n)-Threshold Multisecret Sharing Scheme[END_REF], [START_REF] Pang | A New (t, n)multisecret sharing scheme based on Shamir's secret sharing[END_REF], [START_REF] Yang | A New (t, n)multisecretsharing scheme[END_REF], [START_REF] ¸alkavur | Multisecret-sharing schemes and bounded distance decoding of linear codes[END_REF].

Karnin et al. [START_REF] Karnin | On secret sharing systems[END_REF] examined the following situation. There are k secrets s 1 , s 2 , ..., s k to be shared. Such systems are called [k, m, n] (multisecret-sharing) threshold schemes. Each [k, m, n]-threshold scheme for multisecret-sharing gives k single-secret (m, n)-threshold schemes [START_REF] Ding | Linear multisecret-sharing schemes and error correcting codes[END_REF].

In this work we propose a new multisecret-sharing scheme based on vector spaces over R n . We consider any basis elements of R n . The proposed scheme is ideal and secure by means of the information rate. We use Blakley's method to share and recover the secret. We determine the access structure of this scheme. We analyse its security.

The rest of the paper is organized as follows. In Section II we introduce Blakley's secret sharing scheme. In Section III we present a new multisecret-sharing scheme using Blakley's method and analyse its security. Section IV collects concluding remarks.

II. BLAKLEY'S SECRET SHARING SCHEME

Blakley's secret sharing scheme uses hyperplane geometry to solve the secret sharing problem. The secret is a point in a t-dimensional space and n shares are affine hyperplanes that pass through this point. An affine hyperplane in a tdimensional space with coordinates in a commutative field F can be described by a linear equation of the following form.

a 1 x 1 + a 2 x 2 + ... + a t x t = b.
The intersection point is obtained by finding the intersection of any t of these hyperplanes. The secret can be any of the coordinates of the intersection point or any function of the coordinates. We take the secret to be the first coordinate of the point of intersection.

A. Dealing Phase

Let F be the field we are working on. The dealer generates a secret point x in F , where the first coordinate x [START_REF] Anton | Elementary Linear Algebra[END_REF] is set to the secret value and sets the values of the other coordinates randomly from the field F . The ith user will get a hyperplane equation over F ,

a i1 x 1 + a i2 x 2 + ... + a it x t = y i (1)
For a (t, n)-threshold scheme there will be n such hyperplane equations, and hence we will have an n×t linear system,

AX = Y (2)
The dealer then sends the secret value of y i along with a i1 , ..., a it to user i. The coefficients a ij are not sensitive and can be made publish if needed.

B. Share Combining Phase

Share combining step is simply finding the solution of a linear system of equations. Suppose that a coalition S = {i 1 , ..., i t } of users come together. They form a matrix A s using their hyperplane equations and solve

A s X = Y s , (3)
where Y s is the vector of the secret shares of the users. The secret is found as the first coordinate of the solution [START_REF] Li | A New (t, n)-Threshold Multisecret Sharing Scheme[END_REF].

III. MULTISECRET-SHARING SCHEMES BASED ON VECTOR SPACES

Consider a finite-dimensional vector space. We construct a multisecret-sharing scheme based on this space.

• Let the vector space F n be the secret space.

• Let any vector of F n be the secret. We know that if vector addition is defined to be matrix addition and vector scalar multiplication is defined to be matrix scalar multiplication, then the set W of all n×n matrices with real entries is a vector space.

Definition: If W is any vector space and

V = { v 1 , v 2 , ..., v n } is a set of vectors in W , then V is called a basis for W if the following two conditions hold: i) V is linearly independent, ii) V spans W [1]. Theorem 1. ([1]) If V = { v 1 , v 2 , .
.., v n } is a basis for a vector space W , then every vector v in W can be expressed in the form

v = c 1 v 1 + c 2 v 2 + ... + c n v n
in exactly one way, where c 1 , c 2 , ..., c n ∈ R .

A. Proposed Method

Now we construct a multisecret-sharing scheme based on the vector space F n . Let any vector of F n be the secret S = (s 1 , s 2 , ..., s n ). First we find any basis of F n is called V = { v 1 , v 2 , ..., v n } and calculate the shares y i , i = 1, 2, ..., n.

V.S T = Y T , where V is an element of W . We write the following linear equation system for each participant.

     v 1 v 2 . . . v n      .      s 1 s 2 . . . s n      T =      y 1 y 2 . . . y n     
It can reach by solving this equation system. The dimension of matrix V is n × n and the elements are the basis elements of F n . The transpose of S is denoted by S T . The dimension of S is n × 1. Y is a sharing vector.

Theorem 2. If the matrix V is non-singular, then the secret can be recovered. Otherwise it cannot be reached.

Proof. If the matrix V is singular, then the vectors of V will be linearly dependent. In this case these vectors do not consist of a basis of F n . Theorem 3. In this multisecret-sharing scheme we have the following facts.

1) The qualified set in the access structure is the set of all n participants. 2) No subset of less than n participants can recover the secret.

Proof. 1) The secret is recovered thanks to the basis elements of F n and their number is n.

2) The number of basis elements cannot be less than n.

Only n elements can be used to recover the secret but (n -1) cannot.

Corollary 1. The multisecret-sharing scheme satisfying the hypothesis of the above theorem is also a (n, n)-threshold secret sharing scheme.

Proof. It is clear that the participants are the basis elements of F n and their number is n. These n participants can be reached the secret by combining their shares.

Example 1. Consider the vectors [START_REF] Blakley | Safeguarding Cryptographic Keys[END_REF][START_REF] Blakley | Linear algebra approach to secret sharing schemes[END_REF].

v 1 = (1, 2, 1), v 2 = (2, 9, 0), v 3 = (3,
The set V = {v 1 , v 2 , v 3 } is a basis for F 3 .
• To show that the set V spans

where c 1 , c 2 , c 3 ∈ F or (b 1 , b 2 , b 3 ) = (c 1 + 2c 2 + 3c 3 , 2c 1 + 9c 2 + 3c 3 , c 1 + 4c 3 )
on equating corresponding components

c 1 + 2c 2 + 3c 3 = b 1 , 2c 1 + 9c 2 + 3c 3 = b 2 , c 1 + 4c 3 = b 3 . (4
) Thus, to show that V spans F 3 , we must demonstrate that the system (4) has a solution for all choices of b = (b 1 , b 2 , b 3 ).

• To prove that V is linearly independent, we must show that the only solution of

c 1 v 1 + c 2 v 2 + c 3 v 3 = 0 (5)
is

c 1 = c 2 = c 3 = 0. As above, if (5)
is expressed in terms of components, the verification of independence reduces to showing that the homogeneous system

c 1 +2c 2 +3c 3 = 0, 2c 1 +9c 2 +3c 3 = 0, c 1 +4c 3 = 0 (6)
has only the trivial solution. Observe that systems (4) and ( 6) have the same coefficient matrix. Thus we can prove that in systems ( 4) and ( 6) the matrix of coefficients

A =   1 2 3 2 9 3 1 0 4   .
has a nonzero determinant, that is |A| = -1. So V is a basis for F 3 . Now we try to construct a multisecret-sharing scheme based on F 3 .

Let the secret vector be S = (s 1 , s 2 , s 3 ) = (-3, 1 -1) ∈ F 3 . We calculate the shares as follows.

y T 1 = v 1 S T = (1, 2, 1).(-3, 1, -1) T = -2 y T 2 = v 2 S T = (2, 9, 0).(-3, 1, -1) T = 3 y T 3 = v 3 S T = (3, 3, 4).(-3, 1, -1) T = -10
The participants (the basis elements) can reach the secret by combining their shares as follows.

  1 2 1 2 9 0 3 3 4   .   s 1 s 2 s 3   =   -2 3 -10   .
It is seen that the secret S = (s 1 , s 2 , s 3 ) = (-3, 1, -1) by solving the above linear system. This multisecret-sharing scheme is also a (3, 3)threshold scheme.

B. Information throughput

Since the security of a system depends on the amount of information that must be kept secret, the size of the shares given to the participants is an important point in the design of secret sharing schemes. Besides, if the shares are too large, the memory requirements for the participants will be too strong, and the algorithms used to compute the shares will become inefficient. So one of the basic parameters in secret sharing is the information rate ρ of the scheme, which is defined in [START_REF] Padro | Robust vector space secret sharing schemes[END_REF] to be the ratio between the length (in symbols) of the secret and the maximum length of the shares given to the participants.

A secret sharing scheme is said to be ideal if its information rate is equal to one, which is the maximum possible value. So our scheme is ideal since the secret S [START_REF] Anton | Elementary Linear Algebra[END_REF] is one element of F, like the share Y i of the generic user i.

C. Security Analysis

In a secret sharing scheme, there exists the possibility that some participants lie about the value of their shares in order to obtain some illicit benefit. Therefore, the security against cheating is a key point in the implementation of secret sharing schemes. If F = F q a finite field of order q, one possible attack would be for s participants to collude together, and to guess the shares of the remaining n -s users. It is readily seen that the probability of success of such an attack is 1 q n-s . Thus a safe system would require q to be large.

IV. CONCLUSION

In the present article we have introduced a new multisecretsharing scheme based on vector spaces over the F vector space F n for some field F. The reconstruction algorithm is based on Blakley's method. We determine the access structure, and calculate the information rate of this scheme. We have analysed its information flow. The new system is ideal in a precise sense in terms of information troughput. It is reasonably secure if a sufficiently large field is used.

The choice of a field depends on the practical application and machine implementation. Using a finite field F = F q might require a large size q in view of probabilistic attacks. Using an infinite field like the reals removes this security problem at the price of a possible computational burden and precision issues. So we have preferred to expose our system for a general computation field F.

F 3 ,

 3 we must show that an arbitrary vector b = (b 1 , b 2 , b 3 ) can be expressed as a linear combination b = c 1 v 1 + c 2 v 2 + c 3 v 3 of the vectors in V . Expressing this equation in terms of component gives (b 1 , b 2 , b 3 ) = c 1 .(1, 2, 1) + c 2 .(2, 9, 0) + c 3 .(3, 3, 4),