Théo Cavalieri
email: theo.cavalieri@gmail.com

Jean Boulvert

Vicent Romero-García

Gwénaël Gabard

Marie Escouflaire

Josselin Regnard

Jean-Philippe Groby

Rapid additive manufacturing of optimised anisotropic metaporous surfaces for broadband absorption

Keywords: numbers: 61.43.Gt, 81.05.Rm, 46.40.Ff, 52.38.Dx, 43.20.Gp

Porous materials and metamaterials play a key role in sound absorbing and insulation solutions in acoustics. With the growing interest in additive manufacturing techniques, recent work has focused on the printing of porous and resonant structures for acoustic purposes. Usual metaporous surfaces/interfaces are generally built by periodically inserting resonant elements in an existing porous layer. This complex manufacturing process can be significantly simplified by using additive manufacturing techniques, which also eases the design and optimisation of the metaporous surface.

In this work, the acoustic properties of the metaporous surface are controlled by simple geometric parameters defining both the anisotropic porous layer and the shapes of the resonators. Hence, we focus on optimising split-ring resonators embedded in a micro-treillis porous layer, which are built in a single part using additive manufacturing techniques. A finite-element method together with the Bloch wave decomposition provides a numerical model used to predict the reflection and absorption coefficients under normal incidence. The geometric parameters of the anisotropic metaporous surface are then optimised by non-linear minimisation techniques to maximise acoustic absorption. An optimal metaporous surface is 3D printed by fused-deposition modeling and its acoustic properties are measured in an impedance tube. The measurements are in good agreement with the predicted optimal broadband absorption coefficient. This work demonstrates the benefits of additive manufacturing for designing metaporous acoustic surfaces.

I. INTRODUCTION

Porous materials are commonly used in the context of noise reduction, particularly for sound absorption over wide frequency ranges [START_REF] Allard | Propagation of Sound in Porous Media: Modelling Sound Absorbing Materials[END_REF][START_REF] Cox | Acoustic Absorbers and Diffusers: Theory, Design and Application[END_REF] . When rigidly backed, porous layers provide their strongest absorption for frequencies higher than that of the so-called quarterwavelength resonance and perfect absorption is obtained when critical coupling occurs [START_REF] Jiménez | Perfect absorption of sound by rigidlybacked high-porous materials[END_REF] .

Although efficient for broadband absorption, the use of homogeneous porous layers implies bulky and heavy structures to cope with practical applications. The absorption efficiency of porous layers can be improved at lower frequencies by grading the layer [START_REF] Boulvert | Optimally graded porous material for broadband perfect absorption of sound[END_REF] , or by exploiting the anisotropy of the material for diffuse-field absorption [START_REF] Cavalieri | Graded and anisotropic porous materials for broadband and angular maximal acoustic absorption[END_REF] , but they still do not display deep-subwavelength properties.

Deep-subwavelength properties are on the contrary the specific features of metamaterials or metasurfaces [START_REF] Jiménez | Ultra-thin metamaterial for perfect and quasi-omnidirectional sound absorption[END_REF][START_REF] Romero-García | Perfect and broadband acoustic absorption by critically coupled sub-wavelength resonators[END_REF][START_REF] Groby | The use of slow waves to design simple sound absorbing materials[END_REF][START_REF] Ma | Acoustic metasurface with hybrid resonances[END_REF][START_REF] Christensen | Extraordinary absorption of sound in porous lamella-crystals[END_REF][START_REF] Donda | Extreme low-frequency ultrathin acoustic absorbing metasurface[END_REF] . Nevertheless, deep-subwavelength perfect absorption of metamaterials is usually narrow band, because of the high quality factor of resonators at low frequencies. Detuned resonators [START_REF] Yang | Optimal sound-absorbing structures[END_REF][START_REF] Jiménez | Rainbow-trapping absorbers: Broadband, perfect and asymmetric sound absorption by subwavelength panels for transmission problems[END_REF] are thus used to widen the absorption efficiency of metasurfaces, but they usually do not display absorption as broadband as porous layers.

Metaporous surfaces thus appear as an excellent compromise between the broadband efficiency of porous layer and the deep-subwavelength features of metasurfaces [START_REF] Groby | Enhancing the absorption coefficient of a backed rigid frame porous layer by embedding circular periodic inclusions[END_REF] . They consist in periodically embedding possibly resonant inclusions in a porous layer. Their absorbing properties rely on a complex balance between the thermal and viscous losses in the porous layer and in the resonant elements, the resonances of the inclusions [START_REF] Lagarrigue | Absorption of sound by porous layers with embedded periodic arrays of resonant inclusions[END_REF][START_REF] Groby | Enhancing the absorption properties of acoustic porous plates by periodically embedding helmholtz resonators[END_REF] , and trapped modes [START_REF] Groby | Enhancing the absorption coefficient of a backed rigid frame porous layer by embedding circular periodic inclusions[END_REF][START_REF] Groby | Using simple shape three-dimensional rigid inclusions to enhance porous layer absorption[END_REF] that confine the energy between the inclusions and the rigid backing, and extension of the Wood anomaly to periodic layers [START_REF] Groby | Enhancing the absorption coefficient of a backed rigid frame porous layer by embedding circular periodic inclusions[END_REF] . Although the resonances of the inclusion lead to the enhancement of the porous layer absorption (irrespective of the frequency range), perfect absorption is limited to the viscous regime, except if the resonator is directly connected to the surrounding medium, i.e., the air half space [START_REF] Groby | Enhancing the absorption properties of acoustic porous plates by periodically embedding helmholtz resonators[END_REF] . Initially proposed for transmission problems [START_REF] Groby | Acoustic response of a periodic distribution of macroscopic inclusions within a rigid frame porous plate[END_REF][START_REF] Groby | Acoustic response of a rigid-frame porous medium plate with a periodic set of inclusions[END_REF][START_REF] Deckers | A wave based method to predict the absorption, reflection and transmission coefficient of two-dimensional rigid frame porous structures with periodic inclusions[END_REF] , metaporous surfaces are particularly efficient for subwavelength broadband absorption in reflection problems and have thus seen an increasing interest these last years [START_REF] Zhou | Perfect acoustic absorption by subwavelength metaporous composite[END_REF][START_REF] Yang | Metaporous layer to overcome the thickness constraint for broadband sound absorption[END_REF][START_REF] Xiong | Sound attenuation optimization using metaporous materials tuned on exceptional points[END_REF] . They have been studied extensively and coupled with other resonant elements [START_REF] Groby | Absorption of a rigid frame porous layer with periodic circular inclusions backed by a periodic grating[END_REF] or combined with poroelastic layers, thus enabling the use of elastic resonators [START_REF] Weisser | Acoustic behavior of a rigidly backed poroelastic layer with periodic resonant inclusions by a multiple scattering approach[END_REF] .

Nevertheless, the manufacturing of metaporous surfaces is still tedious and requires inserting resonant elements in porous layer in which the volume of the resonator has been previously removed. The porous material should also be available, thus its properties are subjected to constraints and cannot be freely chosen.

Over the past decade, rapid additive manufacturing methods have made significant advances and complex structures can now be printed with ever finer levels of accuracy. In this article, we optimise a metaporous surface that is then manufactured by fused-deposition modelling (FDM). The complete metasurface (i.e. the porous layer with the embedded resonators) is described by a small number of geometric parameters and numerically modelled by a finite-element method (FEM). Fine-tuning both the micro-geometry of the porous layer and the macro-geometry of the resonant inclusions can lead to an optimal design that provides high levels of acoustic absorption over a chosen frequency range. While the macrogeometry of the metaporous surface is directly optimised, the micro-geometry of the porous layer is optimised via effective acoustic properties calculated using the two-scale asymptotic homogenisation theory [START_REF] Auriault | Homogenization of coupled phenomena in heterogenous media[END_REF][START_REF] Sanchez-Palencia | Non-homogeneous media and vibration theory[END_REF][START_REF] Bensoussan | Asymptotic Analysis of Periodic Structures[END_REF] applied to the representative unit-cell. The optimisation itself is performed using the downhill simplex method with a cost function that relies on the FEM model [START_REF] Nelder | A simplex method for function minimization[END_REF][START_REF] Press | Numerical recipes: the art of scientific computing[END_REF] . The porous layer in which the inclusions are embedded is composed of a transverse isotropic porous material, which prevents the use of multiple-scattering theory (MST), as proposed initially [START_REF] Groby | Acoustic response of a periodic distribution of macroscopic inclusions within a rigid frame porous plate[END_REF] .

The present article is organised as follows: we first introduce the studied geometry and the manufacturing techniques in Sec. II. We then describe the solution procedures and the numerical methods used to model and optimise the metaporous surfaces in Sec. III. In Sec.

IV the validation procedure with experimental techniques is presented. Finally, multiple numerical and experimental results are presented and compared in Sec. V. Internal acoustic fields are also provided to give more physical insight into the role of the anisotropy of the porous material in the resonator coupling. Additionally, detailed validation procedures of the numerical solutions are given in Appendix A.

II. METAPOROUS SURFACE

In this section, we introduce the metaporous surface and its geometry. The additive manufacturing technique used to produce the samples is then presented.

A. Description of the configuration

The design of the metaporous surface is illustrated in Fig. 1. It consists in acoustically rigid 2D split-ring resonators (SRR) periodically embedded in a transverse isotropic porous material. Several split-rings are embedded in the porous layer so as to form a super-cell. The whole medium is described in a two-dimensional Cartesian coordinate system R = (O, e 1 , e 2 ) with position vector x = (x 1 , x 2 ) ∈ R 2 . The metaporous surface has a thickness L 2 = 25 mm in the normal direction e 2 and is of infinite extent in the transverse direction e 1 . The resonant inclusions are placed periodically within the porous layer, with a regular spacing along e 1 .

The periodic medium is described by its lattice constant (or periodicity) L 1 = 84 mm. This corresponds to twice the transverse dimension of the impedance tube used for the experimental validation presented in Sec. IV.

ψ (1) a ψ (2) a h (1) h (2) R (1) out R (1) in R (2) out R (2) in (1) (2) ψ (1) r ψ
(2) r -L 2 0 L 1 /2 L 2 x 1 x 2 Γ 0 Γ L Γ S Γ - Γ + Figure 1. ( Color
online) Schematic of the split-ring resonators embedded in anisotropic porous material. In green is the air, in blue is the porous material Ω e with transverse isotropic properties, set between interfaces Γ 0 and Γ L , with periodic conditions on Γ ± .

As shown in Fig. 1, four resonators are embedded within the porous material, which make up the periodic layer. The unit cell of this periodic medium is composed of two pairs of distinct resonators, and displays a reflection symmetry with respect to the vertical plane Γ S along e 2 . Therefore, we will only focus on one half of the unit cell presented in Fig. 1. Each resonator j ∈ 1, 2 is described by the following geometric parameters: the inner and outer radii R out respectively, the depth h (j) (or center position along e 2 ), the horizontal position (j) = (2j -1)L 1 /8 and finally both aperture and rotation angles, respectively ψ (j) a and ψ (j) r . In this work, the domains will be denoted as follows: the porous material is Ω e , the jth rigid resonator is Ω (j) r , the cavity of the jth resonator is Ω (j) c , and the neck of the jth resonator is Ω (j) n . The metaporous surface is then considered as

Ω := x = (x 1 , x 2 ) | x ∈ Ω e j Ω (j) c ∪ Ω (j) n ∪ Ω (j) r .

B. Manufacturing techniques

Recent manufacturing techniques have been employed to build porous materials [START_REF] Boulvert | Acoustic modeling of micro-lattices obtained by additive manufacturing[END_REF][START_REF] Zieliski | Reproducibility of sound-absorbing periodic porous materials using additive manufacturing technologies: Round robin study[END_REF] . Such materials can exhibit anisotropic properties, which can play in our favour when looking for optimal structures. In order to manufacture such metasurfaces, we employ a fused-deposition modelling (FDM) technique. It consists in the deposition of molten polylactic acid (PLA), in the form of thin filaments. The raw material initially comes in the form of a solid wire of diameter 1.75 mm. It is then heated to approximately 473 K within the nozzle, and deposited in thinner filaments of 0.4 mm in diameter. The structure is built layer after layer, which makes it possible to embed a rigid inclusion in the anisotropic porous material. To do so, the geometry of the optimal structure is first described in the standard tessellation language (STL). It is then translated in series of computer numerical control (CNC) instructions, known as G-code. This operation is readily performed by slicing softwares. The geometry of the porous material is controlled by the filling fraction F of the micro-trellis, which is a common infill pattern in slicing procedures. With a fixed filament diameter, the filling fraction is controlled by the lattice constant of the micro-trellis. The micro-trellis directly provides transverse isotropic properties, which can thus be tuned by changing the value of F. As seen in Fig. 2(a) the sample is placed on the heatbed under the nozzle. The cavities and necks of the resonators embedded in the porous layer can be seen in detail in Fig.

2(b-e).

In order to print the solid resonators, the filling fraction is set to F = 1. The metaporous super-cell is manufactured using the printer Raise3D Pro2. As the structure is build layer-by-layer, the actual manufactured geometry differs slightly from the theoretical geometry shown in Fig. 1. Optical microscopy is employed to provide qualitative analysis of the additive manufacturing defects. In Fig. 3(a) we observe the porous micro-trellis made of orthogonal filaments, with lattice constant 1 mm. The porous material itself presents some defects, mainly the filament striction and surface roughness, which have been studied in previous works [START_REF] Boulvert | Acoustic modeling of micro-lattices obtained by additive manufacturing[END_REF] . A crosssectional view of one of the SRRs is shown in Fig. 3(b), where we can observe surface roughness in the cavity of the resonator, linked to the layer-by-layer manufacturing process.

This induced roughness was found to strongly affect the thermal and viscous losses when of dimension comparable to the spacing between adjacent surfaces [START_REF] Huang | 3d-printed sound absorbing metafluid inspired by cereal straws[END_REF] . A close-up image of the SRR neck is shown in Fig. 3(c), illustrating the limits of FDM when dealing with parts of small dimensions. The combination of these discrepancies has an impact on the overall acoustic behaviour of the metaporous surface.

III. MODELLING AND OPTIMISATION

In order to optimise the metaporous surface, we need a predictive model for its acoustic properties. This is done using a finite-element method combined with the Bloch wave decomposition. The porous material is described as an equivalent fluid to account for viscous and thermal dissipation phenomena. The optimal geometric parameters are identified using a non-linear constrained minimisation of the reflection coefficient.

A. Equivalent fluid model

The porous structure itself can be decomposed into two domains, Ω s and Ω f the solid and fluid phases, respectively. Hence, we denote Ω e ≡ Ω f ∪ Ω s the total domain of the equivalent fluid, the open porosity of this two-phase medium is

φ = |Ω f |/|Ω f ∪ Ω s | = 1 -F.
As the propagation medium consists in a periodic arrangement of a micro-geometric unit cell, the proposed structure can be studied under scale separation [START_REF] Auriault | Homogenization of coupled phenomena in heterogenous media[END_REF][START_REF] Sanchez-Palencia | Non-homogeneous media and vibration theory[END_REF][START_REF] Bensoussan | Asymptotic Analysis of Periodic Structures[END_REF] . At the microscopic scale (i.e. the scale of the micro-geometric unit cell), the rheology of the fluid phase is driven by thermo-acoustic and visco-inertial phenomena. Under a sufficiently large scale separation factor, the macroscopic effective properties of the medium can be described by homogenisation of the local fields. The Johnson-Champoux-Allard-Lafarge (JCAL) model gives an approximation for both thermal and viscous solutions [START_REF] Lafarge | Dynamic compressibility of air in porous structures at audible frequencies[END_REF] .

More recent works on additive manufacturing and sound propagation in graded porous layers have used the same description [START_REF] Boulvert | Optimally graded porous material for broadband perfect absorption of sound[END_REF][START_REF] Cavalieri | Acoustic wave propagation in effective graded fully anisotropic fluid layers[END_REF][START_REF] Terroir | General method to retrieve all effective acoustic properties of fully-anisotropic fluid materials in three dimensional space[END_REF] . The aforementioned model relies on multiple parameters, which estimate the dynamic thermal and viscous permeabilities from an asymptotic description. Three of them are scalars [φ; Λ ; Θ 0 ] and related to thermal effects, respectively the open porosity, the thermal characteristic length and the thermostatic permeability. Three others are direction-dependent [τ ∞ jj ; Λ j ; K 0 jj ] and related to viscous effects, respectively the high-frequency limit of tortuosity, the characteristic viscous length and the viscostatic permeability in the direction e j . Although the transport parameters can be estimated by FEM problems at the scale of the unit cell [START_REF] Terroir | General method to retrieve all effective acoustic properties of fully-anisotropic fluid materials in three dimensional space[END_REF][START_REF] Zieliński | Benchmarks for microstructure-based modelling of sound absorbing rigid-frame porous media[END_REF] , it has been shown that manufacturing defects play an important role in the measured effective properties [START_REF] Boulvert | Acoustic modeling of micro-lattices obtained by additive manufacturing[END_REF][START_REF] Huang | 3d-printed sound absorbing metafluid inspired by cereal straws[END_REF] . In order to account for the internal geometrical defects, we rely on databases obtained from inverse characterisation techniques [START_REF] Boulvert | Optimally graded porous material for broadband perfect absorption of sound[END_REF][START_REF] Niskanen | Deterministic and statistical characterization of rigid frame porous materials from impedance tube measurements[END_REF] . Thus, the JCAL parameters are obtained by fitting experimental results with respect to the filling fraction F.

B. Finite-Element Model

We consider acoustic wave propagation under ambient conditions and in the linear harmonic regime, using the time convention e -iωt . In the frequency domain, the homogeneous anisotropic Helmholtz equation in the two-dimensional Cartesian coordinate system is given by B∇

• ([H]∇p) + ω 2 p = 0 , (1)
with excess acoustic pressure p and angular frequency ω. The effective bulk modulus of the equivalent fluid is denoted B, and

[H] = [ρ] -1
is the inverse mass density tensor. These properties encapsulate the thermo-acoustic and visco-inertial properties of the equivalent fluid, and are complex and frequency-dependent functions.

x 1 In the domain Ω 0 , we define an impinging plane wave in the form of

x 2 Ω (1) c Ω (2) c Ω (1) r Ω
(2) r Ω (1) n Ω (2) n L air L pml L 2 Ω pml Ω 0 Ω e L θ Γ 0 Γ L Γ - Γ + Γ S
p inc (x, ω) = |A|e iϕ 0 e ik•x , (2)
with the wavevector k = {k 1 , k 2 } T where the superscript T denotes the transposition. We have k 1 = k 0 sin(θ) and k 2 = k 0 cos(θ), with k 0 = ω/c 0 the wavenumber in free field, c 0 the speed of sound, and θ the angle of incidence with respect to e 1 . For simplicity, we set the amplitude |A| = 1 and the phase offset ϕ 0 = 0. As the metaporous surface is rigidly backed, the normal particle velocity at the interface Γ 0 must be zero, so v • n| Γ 0 = 0 where v is the vector of acoustic particle velocity. On the interface Γ L between the porous layer and the ambient air, the excess pressure p and the normal velocity v 2 are continuous. The periodicity of the acoustic field along e 1 is enforced between the interfaces Γ + and Γ -by

writing    p| Γ + = p| Γ -e ik B •x B v • n| Γ + = v • n| Γ -e ik B •x B (3)
where the vector for Bloch periodicity is k B = {k 0 sin(θ), 0} T and x B = x| Γ + -x| Γ -. The Sommerfeld radiation condition [START_REF] Sommerfeld | Partial differential equations in physics[END_REF][START_REF] Schot | Eighty years of sommerfeld's radiation condition[END_REF] is given by lim

x 2 →+∞ x 1/2 2 ∂ ∂x 2 -ik 0 p s (x 2 ) = 0 , (4)
so that the scattered pressure p s vanishes far from the interface Γ L . This is implemented under the form of a perfectly matched layer (PML) [START_REF] Berenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF] . In the domains Ω (j) n , we model the thermal and viscous losses in the necks of the resonators. The necks are seen as slits of

dimension d (j) = 2R (j) in sin ψ (j)
a /2 , a chord on the jth circle of inner radius R (j) in and aperture angle ψ (j) a , see Fig. 1. Assuming d (j) λ in the slits, the acoustic wave propagation is described by effective parameters, namely the mass density ρ n and characteristic impedance Z n which are both complex-valued and frequency-dependent [START_REF] Stinson | The propagation of plane sound waves in narrow and wide circular tubes, and generalization to uniform tubes of arbitrary crosssectional shape[END_REF] . The air and PML domains Ω 0 and Ω pml are meshed using quadrilateral elements while the porous layer, the necks and cavities are meshed with triangular elements. In order to capture finer details of the geometry, we use a total of 20 × 10 3 elements so that no element is larger than a fifth of the smallest dimension.

C. Numerical estimation of reflection and absorption coefficients

The scattered pressure field emerging from the interaction with the metaporous surface is quantified by its reflection coefficient. The following procedure enables the calculation of the reflection coefficient from a Bloch wave decomposition. As the scattered pressure field is periodic along the direction x 1 , it is expressed as a sum of Bloch waves:

p s (x, ω) = |A|e iϕ 0 q∈Q R q e ik (q) •x , (5)
where the Bloch wave number q is an integer belonging to

Q = -Q, Q ⊂ Z.
The components of the wavevector k (q) are given along both directions and for each value of q by

k (q) =      k 1 + 2πq/L 1 k 2 0 -k (q) 1 2 1/2      . ( 6
)
The coefficient R q corresponds to the reflection coefficient of the qth Bloch wave. The scattered pressure field that emerges is numerically calculated at every point of the domain Ω 0 . The pressure field is then recorded on the line L shown in Fig. 4 and defined by

x 2 = L 2 + L air /2 and 0 < x 1 < L 1 . The Bloch modes satisfy the following orthogonality relation ∀(q, q ) ∈ Q 2 :

L e i k (q) 1 -k (q ) 1 x 1 dx 1 = |L|δ q,q , (7)
where δ q,q is the Kronecker delta and |L| = L 1 . The reflection coefficient of the qth Bloch wave is given by

R q = 1 |L| e -ik (q) 2 x 2 (L) L p s (x, ω)e -ik (q) 1 x 1 dx 1 , (8)
and the modulus of the total reflection coefficient reads

|R(ω)| 2 = q∈Q |R q | 2 Re k (q) 2 /k 2 . (9)
The acoustic absorption coefficient, which represents the fraction of energy dissipated within the metaporous surface, is given by

α(ω) = 1 -|R(ω)| 2 ∈ R . (10)

D. Optimisation routines

In order to optimise the geometric parameters of the metaporous surface, a constrained non-linear minimisation routine is employed. We use the Nelder-Mead (or downhill simplex) method [START_REF] Nelder | A simplex method for function minimization[END_REF][START_REF] Press | Numerical recipes: the art of scientific computing[END_REF] to look for the optimal design. From an initial set of geometric parameters P 0 , the method takes a series of steps to minimise a cost function. The scheme stops when the convergence is reached (with respect to absolute tolerances), or when the maximum number of iterations is attained. The necks of the SRRs are chosen to be in opposite directions along e 2 with ψ

(1) r = 0 rad and ψ

(2) r = π rad. Also, the center position h (j) of each resonator on the vertical axis e 2 is dependent on the resonator's outer radius. In order to easily manipulate this geometric constraint, we introduce the relative depth h(j) , such as

h (j) = R (j) out + h(j) (L 2 -2R (j)
out ). The vector of parameters to be optimised is written

P := F, ψ (1) a , ψ (2) a , R (1)
in , R

in , R

out , R

out , h(1) , h(2

) T ⊂ R 9 , (11)
the bounds of which are described by P ± :

                   F ∈ ]0, 1[ , ψ (j) a ∈ [0, π/2] , R (j) in ∈ [1 × 10 -3 , max(L 1 /8, L 2 /2)] , R (j) out ∈ [1 × 10 -3 , max(L 1 /8, L 2 /2)] , h(j) ∈ [0.1, 0.9] . (12)
Moreover, some constraints are imposed over the geometric parameters, in order to satisfy the overall topology of the system shown in Fig. 1. The constraints are established as the outer perimeter of the resonator is contained in Ω e , and the wall thickness of the resonators must be greater than 800 microns. This translates into

         ∀j = 1, 2 , R (j) out -R (j) in ≥ 8 × 10 -4 , R (j) out ≤ min(L 1 /8, L 2 /2) . (13)
The acoustic absorption coefficient is maximised within a certain frequency range, when the cost function gets closer to zero. In order to do so, we define a cost function that translates how far the current vector of parameters is from the optimal value α = 1. The cost function is the arithmetic average of the absorption within the frequency range:

J(P) := 1 - 1 N f N f i=1 α(f i ) ∈ ]0, 1[ , (14)
with N f the number of frequencies. The minimisation procedure requires an initial guess, which is chosen by setting the vector P 0 with random parameters, following the uniform distribution

P (i) 0 := U P (i) -, P (i)
+ . In this work, we focus on the frequency range between f min = 1000 Hz and f max = 4000 Hz, with N f = 10.

IV. EXPERIMENTAL SETUP

As shown in Fig. 1 the unit cell of the periodic system is symmetric, hence the structure behaves as if half the unit cell was bounded by rigid walls (acting as perfect mirrors) when excited by an plane wave at normal incidence [START_REF] Groby | Total absorption peak by use of a rigid frame porous layer backed by a rigid multi-irregularities grating[END_REF] . Thus, the Bloch modes on the total unit cell are identical to those within rigid boundaries assuming one half of the symmetric unit cell is placed in an impedance tube. In this section, we describe the experimental setup where the proposed metaporous surface is rigidly backed in an impedance tube. We give details on the collecting and processing of experimental data used to retrieve the absorption coefficient.

A. Impedance tube measurement

The waveguide has a square cross-sectional area of 42×42 mm 2 and is made of aluminum.

The cut-off frequency of the first transverse mode is thus f c = c 0 /0.084 ≈ 4060 Hz. A copper plug is backing the metaporous surface at position x 2 = -L 2 ensuring the Neumann boundary condition.

The pressure field is evaluated at N = 11 non-equidistant points located along x 2 . We use a GRAS 46BD 1/4" CCP Pressure Standard Microphone Set, together with the conditioner GRAS Power Module Type 12AQ. The microphone is moved along the e 2 axis, driven by the stepping motor Nanotec TYPE 4H8618M1408, and controlled using CHARLYROBOT F4350 CERNEX.

Robot x

(j) 2 The frequency range is spanned linearly with the emission of N f = 200 monochromatic plane waves. The pressure field is averaged over M = 2 9 = 512 periods, where the period is T = 2π/ω. We use the analyser Stanford SRS MODEL SR785 for the input electrical signal, which is finally transformed into a plane wave acoustic field by a loudspeaker. Since the rod used to move the microphone is at the end of the tube, the source is positioned to the side as shown in Fig. 5. This change is not of predominant importance as only low-frequency plane-waves are considered, below the cut-off frequency of the tube.

0 x 2 x 1 Γ 0 Γ L L 2 Π - Π + GPIB RS-232 s out

B. Data processing

For frequencies below the cut-off frequency f c , the pressure at point j and angular frequency ω is the sum of the forward and backward plane waves:

p j (ω) = Π + (ω)e ik 0 x (j) 2 + Π -(ω)e -ik 0 x (j) 2 , ( 15
)
where Π ± are the complex amplitudes of the forward and backward propagating waves. The measurements at the different points j can be combined in matrix form:

               p 1 p 2 . . . p N                =         e ik 0 x (1) 2 e -ik 0 x (1) 2 e ik 0 x (2) 2 e -ik 0 x (2) 2 . . . . . . e ik 0 x (N ) 2 e -ik 0 x (N ) 2            Π + (ω) Π -(ω)    ⇔ p = [E]Π . (16)
The complex amplitudes Π + (ω) and Π -(ω) can be retrieved by inverting the system of Eq. ( 16),

Π = [E]p . ( 17
)
As the system is of dimension dim [E] = 2 × N , its pseudo-inverse is computed using the generalised Moore-Penrose inverse [START_REF] Dresden | The fourteenth western meeting of the american mathematical society[END_REF][START_REF] Penrose | A generalized inverse for matrices[END_REF] . It is based on the singular value decomposition the left-and right-singular vectors respectively. The superscripts and ˜denote the complex conjugation and pseudo-inverse respectively. The complex reflection coefficient is obtained from the ratio of the two complex amplitudes:

(SVD) 30 so that [E] = V[S] -1 U ,
R(ω) = Π + /Π -∈ C , (18)
and the absorption coefficient is given by Eq. ( 10).

V. RESULTS AND DISCUSSIONS

In this section, results are presented in terms of absorption, comparing numerical predictions with experimental data. Moreover, the internal acoustic fields for the pressure p and particle velocity v are discussed based on the FEM simulations.

A. Maximised broadband absorption

We optimise the aforementioned macro-and micro-geometries, as defined by Eq. ( 11), so as to maximise the acoustic absorption between f min and f max . The predicted and measured acoustic absorption of the optimal metaporous surface are given in Fig. 6(c). For comparison, we also show the absorption results for two other cases: first, the porous material alone in the absence of inclusion, and a case where the SRRs are replaced by rigid circular inclusions (non-resonant). These configurations are investigated numerically and experimentally in Fig. 6(a) and Fig. 6(b) respectively. The geometric parameters of the optimised metaporous surface are given in Tab. I. Looking at Fig. 6(a,b), the experimental results and simulated predictions are in very good agreement. We notice small discrepancies which can emerge from a combination of distinct phenomena in Fig. 6(c). The viscous and thermal losses inside the SRR cavities are not accounted for in the numerical simulations, and as discussed in Sec. II the manufacturing process induces defects to the micro-structure of the porous material and to the geometry of the resonators. When the physics of the system is driven at the local scale (in the presence of SRR), it is very sensitive to small manufacturing defects. This is the reason why the discrepancies are much more visible in Fig. 6(c) for the SRR, than in Fig. 6(a,b) for the porous layer and the non-resonant inclusions.

F ψ (1) a (deg) ψ (2)
a (deg) R

in (mm) R

in (mm) R

out (mm) R

0.40 33 10 5.3 6.1 6.1 6.9 0.35 0.21 Table I. Geometric parameters of the optimised metaporous surface.

Experimentally the system with non-resonant circular inclusions reaches quasi-perfect absorption α ≈ 0.985 at the frequency f = (2100 ± 20) Hz in Fig. 6(b). This absorption peak reveals the presence of a trapped mode that confines the acoustic energy between the rigid backing and the periodic inclusions. For the SRR configuration in Fig. 6(c), three distinct peaks of absorption are observed. Over the frequency range of optimisation, we experimentally reach a mean absorption of α avg ≈ 0.82, with quasi-perfect absorption α ≥ 0.999 at f = (1720 ± 20) Hz. This corresponds to a sub-wavelength regime as L ≈ λ/8.

To gain more insight into the physical mechanisms responsible for these three peaks, it is useful to analyse the corresponding fields of acoustic pressure and particle velocity provided by the FEM model and presented in Fig. 7.

B. Internal fields and couplings

When excited at normal incidence and in the absence of inclusion, the transverse particle velocity v 1 is zero. In this case, only the properties in the normal direction e 2 of the equivalent fluid are contributing to the acoustic behaviour of the homogeneous porous layer. In our case, the periodic inclusions introduce resonances and trapped modes that couple the velocity field in both its normal and transverse directions. The FEM provides the complete solution to the anisotropic Helmholtz equation within the metaporous surface and one can retrieve the internal fields. We use these fields to better understand the behaviour of the metaporous surface illustrated in Fig. 6. More precisely, we focus on the first peaks at f 1 = (1700±50) Hz and f 3 = (3400 ± 50) Hz which correspond to the resonances of the SRRs. The peak at f 2 = (2850 ± 50) Hz is the trapped mode which increases the energy located between the inclusions and the rigid backing. From the results in the aforementioned Fig. 6, the absolute scattered pressure |p s | is minimised in Ω 0 at those frequencies (since the reflection is minimised). Along with the pressure field, the normal and transverse components of the acoustic particle velocity shown in Fig. 7 exhibit the couplings introduced by the SRRs. Regarding the mode at f 1 , we see in Fig. 7(a) that the largest values of pressure are confined within Ω

0 L 1 /4 L 1 /2 0 -L 2 /2 -L 2 Mode at f 1 Pressure |p| 0 L 1 /4 L 1 /2 0 -L 2 /2 -L 2 Normal velocity |v 2 | 0 L 1 /4 L 1 /2 0 -L 2 /2 -L 2 Transverse velocity |v 1 | 0 L 1 /4 L 1 /2 0 -L 2 /2 -L 2 Mode at f 2 0 L 1 /4 L 1 /2 0 -L 2 /2 -L 2 0 L 1 /4 L 1 /2 0 -L 2 /2 -L 2 0 L 1 /4 L 1 /2 0 -L 2 /2 -L 2 Mode at f 3 0 2 4 5 Pressure (Pa) 0 L 1 /4 L 1 /2 0 -L 2 /2 -L 2 10 -4 10 -3 10 -2 10 -1 Particle velocity (m.s -1 ) 0 L 1 /4 L 1 /2 0 -L 2 /2 -L 2 (a) (b) (c) (d) (e) (f) (g) (h) (i)
c , which correspond to the resonance of this SRR at that particular frequency. This resonance induces a strong particle velocity in the neck Ω in Fig. 7(g). Again, the high values of normal particle velocity in Fig. 7(h) highlight this resonance phenomena at frequency f 3 .

Although the metaporous surface is excited by a monochromatic plane wave at normal incidence, the resonators introduce a strong coupling in the porous domain Ω e between the normal and transverse directions. The Figs. 7(b,c,e,f,h,i) show the decimal logarithm of normal and transverse velocity fields within the metaporous surface. We notice that both transverse and normal velocity fields are of the same orders of magnitude, between 10 -4 m.s -1 and 10 -1 m.s -1 . Thus, the normal and transverse properties of the equivalent fluid, which are encapsulated in the mass density tensor [ρ], are both contributing to the overall acoustic behaviour of the metaporous surface.

VI. CONCLUSIONS

Porous materials and metamaterials are commonly used in acoustics, for absorption and insulation purposes. In this work, we propose an optimal design of resonant metaporous surface for acoustic absorption purposes. The proposed solution consists in the periodic embedding of multiple split-ring resonators in a transverse isotropic porous layer. By tuning the macro-geometry of rigid resonators as well as the micro-geometry of the porous medium, we maximise the absorption in the frequency range between 1 and 4 kHz. The metaporous surface provides strong broadband absorption capabilities and subwavelength quasi-perfect absorption.

The numerical simulations based on FEM provide a predictive model for the reflection and absorption coefficients. This predictive model is used in the cost function for the minimisation algorithm, whose goal is to maximise the absorption in a target range of frequencies.

Experiments have been conducted in an impedance tube with a square cross section.

They provide a validation of the numerical model, as the experimental data are in good agreement with the simulated results. Additive manufacturing techniques have proved their ability to build resonant metaporous surfaces.

The present work demonstrates the capabilities of metaporous surfaces for absorbing sound waves, and highlights the use of additive manufacturing for such acoustic treatments.

While the design and manufacturing of materials incorporating inclusions is usually tedious, the solutions presented here are fabricated in a single stage. In addition to offering a solution that is easy and quick to implement, this first prototype was developed at low cost, although this is slightly higher than for conventional porous materials.

Several points of this work provide perspectives for future work. Accounting for the oblique incidences would make it possible to obtain optimal structures for applications in diffuse acoustic fields, which however requires particular attention to the experimental validation. Finally, the optimisation of the shape and number of inclusions could offer more possibilities regarding the compromise between structural dimensions and frequency range of absorption.

APPENDIX A: VALIDATION OF THE NUMERICAL MODEL

The FEM model is compared to two other approaches in order to validate the simulation results when considering anisotropic porous materials and in presence of inclusions.

A. Transfer matrix method

The acoustic wave propagation in a homogeneous porous medium can be derived using the well known transfer matrix method (TMM) [START_REF] Brouard | A general method of modelling sound propagation in layered media[END_REF] . Using the state-vector formalism, applied to the acoustic pressure and the normal component of velocity we introduce W = {p, v 2 } T .

The acoustic fields satisfy the following differential system of equations

∂ x 2 W -iω   0 H -1 22 B -1 eq 0   W = 0 , (19)
with B eq the equivalent bulk modulus, possibly accounting for out-of-plane rotations and oblique incidence. H 22 is the apparent inverse mass density in the normal direction [START_REF] Cavalieri | Acoustic wave propagation in effective graded fully anisotropic fluid layers[END_REF][START_REF] Terroir | General method to retrieve all effective acoustic properties of fully-anisotropic fluid materials in three dimensional space[END_REF] .

As the normal particle velocity v 2 is imposed to be zero on the interface Γ 0 , we have the boundary condition W| Γ 0 = {p, 0} T . The state vector can then be calculated using matrix exponentials:

   1 + R (R -1)/Z    Γ L = [M ]    p(0) 0    Γ 0 , (20)
with Z = (ρ 0 c 0 )/ sin(θ) and [M ] the matricant [START_REF] Pease | Methods of Matrix Algebra[END_REF][START_REF] Horn | Matrix analysis[END_REF] . From the differential system of equations, the reflection coefficient is obtained as

R = M 11 M 21 Z + 1 / M 11 M 21 Z -1 , (21)
and the absorption coefficient is α = 1 -|R| 2 . The JCAL parameters of the homogeneous isotropic porous material are given in Tab. II. For this case of numerical validation, we consider a porous layer of thickness L 2 = 25 mm.

φ Θ 0 [×10 -8 m 2 ] τ ∞ 1 τ ∞ 2 Λ 1 [×10 -4 m] Λ 2 [×10 -4 m] Λ p [×10 -4 m] K 0 1 [×10 -9 2 ] K 0 2 [×10
The results in Fig. 8 show the comparison of FEM and TMM results for different angles of incidence. The results are in excellent agreement, showing that both methods can be to describe the propagation of waves in an homogeneous and anisotropic porous material.

B. Multiple scattering theory

In the particular case where the porous layer is isotropic and the inclusions are cylindrical, the multiple scattering theory (MST) provides an analytical description of the acoustic fields. As this theory is not the centrepiece of this work, it is only briefly presented in this appendix for validation purposes. The total acoustic field is a superposition of the undisturbed (meaning incident, in absence of inclusions) and scattered fields [START_REF] Romero-García | Fundamentals and Applications of Acoustic Metamaterials[END_REF] . The scattered field emerges from the collective effect produced by the interactions of the undisturbed field with every obstacle. In response to the undisturbed acoustic excitation, every obstacle contributes to the global scattered field. For every obstacle, a local description of the scattered field in polar coordinates system leads to cylindrical wave functions, i.e. Bessel and Hankel functions. The Bloch wave decomposition provides a solution for the total field as previously explained in Sec. III. In order to compare the FEM results to the MST solutions, we consider a simple geometry shown in Fig. 9. Table III. JCAL parameters of the homogeneous isotropic porous material.

L 2 L 1 2r L 1 /2 L 2 /2 Ω j,0 Ω j,1 Ω j,-1
The results in Fig. 8 show the comparison of FEM and MST results for different angles of incidence. The results provided by the FEM and the MST are in excellent agreement, thus demonstrating the Wood anomaly at high frequencies, in Fig. 10(c,d).

Figure 2 .

 2 Figure 2. (Color online) (a) Photography of printing machine in operation. (b,c,d,e) Samples of the resonant porous structures. The red arrow in (e) represents the direction of incidence of the acoustic wave.

Figure 3 .

 3 Figure 3. Optical microscopy images of (a) top-view the porous micro-structure, (b) the split-ring resonator and (c) side view of the SRR neck.

Figure 4 .

 4 Figure 4. (Color online) Schematic of the finite-element model. The porous material is shown in blue, the air in green, the PML in yellow, the slits in red, and the rigid SRR in white. The interfaces are shown in dashdotted lines and the probing line-segment in dashed black line.

Figure 5 .

 5 Figure 5. (Color online) Schematic of the experimental setup.

  where [S] is the matrix of singular values, [U ] and [V ]

Figure 6 .

 6 Figure 6. Measured and predicted absorption coefficients. (a) Porous material alone, (b) porous material with circular inclusions, (c) porous material with SRR.

Figure 7 .

 7 Figure 7. (Color online) Internal fields for acoustic pressure (a,d,g), transverse (b,e,h) and normal (c,f,i) particle velocity. Top row (a,b,c) shows the fields for the Mode 1 at f = 1650 Hz, middle row (d,e,f) shows the fields at Mode 2, with f = 2850 Hz and bottom row (g,h,i) shows the fields at Mode 3, f = 3400 Hz. The colorbar on the bottom-left corner relates the pressure fields(a,d,g) while the colormap on the bottom-right corner is attributed to velocity fields (b,c,e,f,h,i). The latter is intentionally saturated as the lower bound is limited to three orders of magnitude below the maximal value of the field.

( 2 )n

 2 of the resonator, as shown in Fig. 7(b). At frequency f 2 , the excess pressure is trapped between the rigid backing at Γ 0 and the rigid SRR Ω (1) r . This trapped mode is clearly visible on the pressure field in Fig. 7(d) as the strong values of pressure indicates the presence of this localised mode. For the mode at frequency f 3 , we observe large values of pressure within the SRR cavity Ω

Figure 8 .

 8 Figure 8. (Color online) Comparison of FEM and TMM results, for angles of incidence θ = π/2 (a), θ = π/3 (b), θ = π/4 (c), and θ = π/6 (d).

Figure 9 .

 9 Figure 9. (Color online) Geometry of a one-dimensional periodic lattice with a cylindrical inclusion embedded in an isotropic homogeneous porous material.

Figure 10 .

 10 Figure 10. (Color online) Comparison of FEM and MST results, for angles of incidence θ = π/2 (a), θ = π/3 (b), θ = π/4 (c), and θ = π/6 (d).

  ABBREVIATIONS

Table II .

 II JCAL parameters of the homogeneous transverse isotropic porous material in the directions e 1 and e 2 .

								-8 m 2 ]
	0.800	2.538	1.139 1.864	2.758	3.710	3.799	4.669	1.812

ACKNOWLEDGMENTS

This research was funded by the Chaire Industrielle MACIA (ANR-16-CHIN-0002) funded by the Safran Group and the ANR (Agence Nationale de la Recherche, France).

DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.


Sheng Gou

Zongyu Yue

Kaichang Di

Roberto Bugiolacchi

Hua Zhu

Patrick C Pinet

Zhanchuan Cai

Mare basalt flooding events surrounding Chang'e-4 landing site as revealed by Zhinyu crater ejecta

Keywords: Zhinyu crater, Von Kármán crater, Mare flooding, Homogenous composition, Single magma chamber, Highlights

Von Kármán crater's floor was flooded with mare basalts during the Imbrian period. This site is the target of China's ongoing Chang'e-4 mission that includes the Yutu-2 rover and its instrumental payload. The Zhinyu crater, one of the largest craters within a few tens of kilometers from the landing site, is the product of a fresh impact that excavated subsurface, basaltic materials onto the surface. The compositional characteristics of the continuous ejecta around the crater vary radially with distance, suggesting possible mineral heterogeneity at depth. At least three main mare flooding phases could have occurred within Von Kármán, producing a basalt layer of at least 320 m in thickness. However, a broadly consistent olivine composition (Fo#: 55) of the three proposed basalt layers suggests that the sources of the mare basalts would have been similar in composition, perhaps even consisting of a single magma chamber where magma had time to evolve. The Yutu-2 rover that is edging its way westwards, could in principle test this hypothesis, thus further constraining the thermal history of the Von Kármán crater.

Introduction

The Von Kármán crater (~186 km in diameter), lying in the northwest of South Pole-Aitken (SPA) basin on the lunar farside, was formed 3.97 Ga ago in the pre-Nectarian period [START_REF] Stuart-Alexander | Estimating modal abundances from the spectra of natural and laboratory pyroxene mixtures using the modified Gaussian model[END_REF][START_REF] Yingst | Geologic Mapping of the Planck Quadrangle of the Moon[END_REF]. During the Imbrian period, the floor of Von Kármán crater was flooded with several basaltic lava flows (or mare basalts) [START_REF] Huang | Geological Characteristics of Von Kármán Crater, Northwestern South Pole-Aitken Basin: Chang'E-4 Landing Site Region[END_REF][START_REF] Pasckert | Lunar farside volcanism in and around the South Pole-Aitken basin[END_REF][START_REF] Ling | Composition, mineralogy and chronology of mare basalts and non-mare materials in Von Kármán crater: Landing site of the Chang'E-4 mission[END_REF]). China's

Chang'e-4 probe touched down successfully in the Von Kármán crater at 177.588°E, 45.457°S on January 3, 2019, and released the Yutu-2 rover to conduct in situ measurements on this previously unexplored region (Di et al., 2019a;[START_REF] Pinet | The Moon's mantle unveiled[END_REF].

The surface of the Chang'e-4 landing site is generally thought to be mantled by ejecta from the Finsen crater (e.g., [START_REF] Gou | Lunar deep materials observed by Chang'e-4 rover[END_REF][START_REF] Hu | Mineral abundances inferred from in-situ reflectance measurements of Chang'E-4 landing site in South Pole-Aitken basin[END_REF][START_REF] Gou | Forsteritic olivine and magnesium-rich orthopyroxene materials measured by Chang'e-4 rover[END_REF][START_REF] Huang | Diverse rock types detected in the lunar South Pole-Aitken Basin by the Chang'E-4 lunar mission[END_REF][START_REF] Lin | Olivine-norite rock detected by the lunar rover Yutu-2 likely crystallized from the SPA impact melt pool[END_REF], though the estimation of ejecta thickness from Finsen crater at the landing site is debated [START_REF] Huang | Geological Characteristics of Von Kármán Crater, Northwestern South Pole-Aitken Basin: Chang'E-4 Landing Site Region[END_REF]Di et al., 2019b;[START_REF] Fu | The subsurface structure and stratigraphy of the Chang'E-4 landing site: orbital evidence from small craters on the Von Kármán crater floor[END_REF]. Zhinyu (176.15 °E,45.34 °S), a fresh, simple crater with a diameter of ~3.8 km and a depth of ~ 800 m (Fig. 1) is one of the largest impacts near the Chang'e-4 landing site, lying about 30 km away to the west. The Zhinyu impact penetrated the mare basalts within the Von Kármán crater, thus it opens an exploratory window into the composition and stratigraphy of the subsurface. Hypervelocity impacts exhume materials from the subsurface and deposit them both within and outside the newly formed crater, allowing for investigation. For instance, [START_REF] Fu | The subsurface structure and stratigraphy of the Chang'E-4 landing site: orbital evidence from small craters on the Von Kármán crater floor[END_REF] targeted 40 dark-haloed craters and 77 bright ray craters around the Chang'e-4 landing area and constructed a possible stratigraphy (including a low-Ti mare basalt layer) based on the excavation depths of these craters and their distances from the Finsen crater. Further, observed differences in the ejecta layers would represent geological variations with depth at the target point. For example, [START_REF] Ding | Fragments Delivered by Secondary Craters at the Chang'E-4 Landing Site[END_REF] proposed that several fresh craters surrounded by blocky fragments encountered by the Yutu-2 rover might be secondaries of the Zhinyu crater. Since the primary source for the regolith at the Chang'e-4 landing site is reckoned to represent Finsen crater ejecta (compositionally different from mare basalt), caution must be applied in interpreting the in situ measured spectra around these secondary craters. Therefore, the Zhinyu crater ejecta provides a rare opportunity to explore the sub-regolith mare basalts, and should also contribute to the correct interpretation of the in situ measured data by the Yutu-2 rover.

Data

Moon mineralogy mapper spectra

The Moon Mineralogy Mapper (M 3 ), which flew onboard India's Chandrayaan-1 lunar orbiter, was a push-broom imaging spectrometer operating from the visible into the near-infrared (0.42 -3.0 μm). M 3 observed the lunar surface both in a highresolution target mode (80 m/pixel; 260 spectral channels) and in a low-resolution global mode (140 m/pixel; 85 spectral channels) from a 100 km orbit [START_REF] Pieters | The Moon mineralogy mapper (M3) on Chandrayaan-1[END_REF]. The M 3 Level 2 spectral reflectance images (REFIMG), which have been photometrically calibrated and thermal emission corrected [START_REF] Clark | Thermal removal from near-infrared imaging spectroscopy data of the Moon[END_REF][START_REF] Besse | One Moon, many measurements 2: Photometric corrections[END_REF], are used in this study to extract spectra from regions of interest (ROIs) for further detailed analysis (Supplementary Fig. 1).

Multiband imager products and images

The multiband imager (MI) was a high-resolution imaging camera onboard the Japanese Kaguya lunar orbiter that included separate visible and near infrared sensors [START_REF] Ohtake | Scientific objectives and specification of the SELENE Multiband Imager[END_REF]. It acquired push-broom images in five bands for the visible and four bands for the near-infrared spectral ranges [START_REF] Ohtake | Scientific objectives and specification of the SELENE Multiband Imager[END_REF]. The spatial resolution of the MI images is 20 m/pixel for the visible bands and 62 m/pixel for the near-infrared bands taken at a nominal orbit altitude of 100 km [START_REF] Ohtake | Scientific objectives and specification of the SELENE Multiband Imager[END_REF].

The mineral abundance products derived from topographically-corrected MI reflectance data [START_REF] Lemelin | Global map products from the Kaguya Multiband Imager at 512 ppd[END_REF] are used in this study to evaluate compositional variabilities of the mare basalts. As concentrations of TiO2 and FeO are useful in classifying major lunar rock types [START_REF] Lucey | FeO and TiO2 concentrations in the South Pole-Aitken basin: Implications for mantle composition and basin formation[END_REF], derived abundances around the Zhinyu crater were produced using the algorithm of [START_REF] Otake | Lunar Iron and Titanium Abundance Algorithms Based on SELENE[END_REF] to investigate possible differences.

Methods

Basalt layer identification

There is evidence that the lava infill of the Von Kármán crater occurred during the Imbrian (Im) period [START_REF] Huang | Geological Characteristics of Von Kármán Crater, Northwestern South Pole-Aitken Basin: Chang'E-4 Landing Site Region[END_REF]. Inferred TiO2 content is commonly used as a classifier for basalts: very low Ti (<1.5 wt% TiO2), low-Ti (1.5-6 wt% TiO2), and high-Ti (>6 wt% TiO2) [START_REF] Hiesinger | Chapter 23 -The Moon[END_REF]. The overall TiO2 content of the continuous ejecta of Zhinyu crater is relatively low (see section 4.1), making it not suitable for subdivision. Consequently, to evaluate the mineral composition of the ejecta around Zhinyu crater, mineral maps (Fig. 2) of clinopyroxene (CPX), orthopyroxene (OPX), olivine (OL) and plagioclase (PLG) were extracted from the MI derived global map products [START_REF] Lemelin | Global map products from the Kaguya Multiband Imager at 512 ppd[END_REF]. It emerges that three rings with varied mineral abundance, in particular for the OL (Fig. 2a), are distributed concentrically around the Zhinyu crater. The boundaries of these three rings are broadly outlined by the OL variation [START_REF] Lai | First look by the Yutu-2 rover at the deep subsurface structure at the lunar farside[END_REF] and defined in this study as Im1, Im2 and Im3, in the order from the rim outward (Fig. 2). Previous studies have proposed buried mare basalts fields inside the Von Kármán crater [START_REF] Pasckert | Lunar farside volcanism in and around the South Pole-Aitken basin[END_REF][START_REF] Ling | Composition, mineralogy and chronology of mare basalts and non-mare materials in Von Kármán crater: Landing site of the Chang'E-4 mission[END_REF]. Here the spectral/mineralogical variations suggest there might have been at least three distinct mare flooding events with a degree of compositional differences.

Olivine composition analysis

The olivine fraction content of igneous rocks is highly diagnostic for interpreting the magma source and the degree of petrologic evolution [START_REF] Isaacson | Remote compositional analysis of lunar olivine-rich lithologies with Moon Mineralogy Mapper (M 3 ) spectra[END_REF]. Based on the established relation between laboratory spectral absorption features and OL composition [START_REF] Sunshine | Determining the composition of olivine from reflectance spectroscopy[END_REF][START_REF] Isaacson | Deconvolution of lunar olivine reflectance spectra: Implications for remote compositional assessment[END_REF], the relative Mg/(Mg + Fe) ratio, as known as the forsterite number (Fo#), can be estimated from the remotely-sensed spectra. Therefore, the average spectrum of each layer was extracted from M 3 images (Supplementary Fig. 1) and deconvolved by the Modified Gaussian Model (MGM) (Supplementary Figs. 234) [START_REF] Stuart-Alexander | Estimating modal abundances from the spectra of natural and laboratory pyroxene mixtures using the modified Gaussian model[END_REF] with different mineral combination configurations (Supplementary Table 1) to estimate mineral modes [START_REF] Clénet | A new systematic approach using the Modified Gaussian Model: Insight for the characterization of chemical composition of olivines, pyroxenes and olivine-pyroxene mixtures[END_REF][START_REF] Clénet | A systematic mapping procedure based on the Modified Gaussian Model to characterize magmatic units from olivine/pyroxenes mixtures: Application to the Syrtis Major volcanic shield on Mars[END_REF]. The OL composition is determined by fitting the MGM-derived OL-related absorption centers to a laboratoryderived compositional trend line [START_REF] Sunshine | Determining the composition of olivine from reflectance spectroscopy[END_REF][START_REF] Isaacson | Deconvolution of lunar olivine reflectance spectra: Implications for remote compositional assessment[END_REF]. The uncertainty of the predicted Fo# values that arises from model initial conditions and other error sources, e.g., the low signal-to-noise ratio and low resolution of M 3 spectrum, is estimated to be about 10% Fo# # and possibly more in the case of ternary mineral mixtures [START_REF] Isaacson | Deconvolution of lunar olivine reflectance spectra: Implications for remote compositional assessment[END_REF][START_REF] Pinet | Characterization of the Olivine/Plagioclase Mineralogy at Copernicus Crater from MGM Deconvolution of M3 Observations[END_REF][START_REF] Pinet | Reassessing the relationship between olivine composition and reflectance Spectroscopy from advanced MGM deconvolution[END_REF].

Basalt layer thickness estimation

For a lunar simple crater with a diameter smaller than 15 km, the maximum excavation depth (Hexc) is approximately 1/10 of the transient crater diameter (Dt), which is about 0.84 times the rim-to-rim diameter (D), i.e., Dt = 0.84D [START_REF] Melosh | Impact Cratering: A geologic process[END_REF].

The ejecta apron of an impact crater is composed of target materials excavated during the excavation stage [START_REF] Melosh | Impact Cratering: A geologic process[END_REF]. They are deposited around the crater following an inverted pre-impact stratigraphy of the underlying bedrock: the materials from the shallower depths tend to be thrown farther away from the crater rim than materials excavated from deeper layers [START_REF] Thomson | A new technique for estimating the thickness of mare basalts in Imbrium Basin[END_REF]. The ejecta becomes thinner and more mixed with local mature regolith with increasing radial distance [START_REF] Melosh | Impact Cratering: A geologic process[END_REF].

The average radius of continuous ejecta blanket (Rce) can be estimated from crater radius (R) using the relation R ce = (2.348 -0.454 +0.564 )R 1.006 [START_REF] Moore | Multiringed basins -illustrated by Orientale and associated features[END_REF].

Although many laboratory experiments and numerical simulations have studied impacts into layered targets (e.g., [START_REF] Stöffler | Experimental hypervelocity impact into quartz sand: Distribution and shock metamorphism of ejecta[END_REF][START_REF] Senft | Modeling impact cratering in layered surfaces[END_REF][START_REF] Prieur | Formation of Simple Impact Craters in Layered Targets: Implications for Lunar Crater Morphology and Regolith Thickness[END_REF], little work has been done to specifically address the relationship between the radial position of ejecta on the upper surface of the continuous ejecta blanket and depth of excavation in layered targets (e.g., [START_REF] Baratoux | Mineralogical structure of the subsurface of Syrtis Major from OMEGA observations of lobate ejecta blankets[END_REF][START_REF] Thomson | A new technique for estimating the thickness of mare basalts in Imbrium Basin[END_REF].

As a first approximation in this study, it's assumed that continuous ejecta at the farthest

Results and discussions

Composition variations

Results of the MI-derived mineral abundances [START_REF] Lemelin | Global map products from the Kaguya Multiband Imager at 512 ppd[END_REF] for the three layers (in the order Im1, Im2, and Im3) show that the PLG abundance increases from 38.2% to 43.1%, the CPX decreases from 26.3% to 18.8%, the OPX increases from 22.6% to 25.2%, and the OL decreases from 12.9% to 10.5% and then increases to 13.2%. The mineral abundances almost fall within each other's standard deviation (Fig. 3), suggesting variations to be marginal. For example, the average CPX abundances for Im1, Im2 and Im3 are ~26.3±5.3%, ~24.8±6.9% and ~18.8±6.2%, respectively. Also, the average OL abundances are ~12.9±2.6%, ~10.5±4.7% and ~13.2±3.4%, respectively. MGM fitting results reveal that Fo# values derived from OL-related absorptions are ~55, suggesting that the OL to be relatively homogenous within the sampled units. , 3.38 Ga, 3.57 Ga, and 3.72 Ga (Ling et al., 2019).

These AMAs indicate there might be older buried lava flows in the Von Kármán crater [START_REF] Pasckert | Lunar farside volcanism in and around the South Pole-Aitken basin[END_REF][START_REF] Ling | Composition, mineralogy and chronology of mare basalts and non-mare materials in Von Kármán crater: Landing site of the Chang'E-4 mission[END_REF], and the composition variations (Fig. 3) highlighted in this study indicate the recent erupted basalt (Im3) might have a lower CPX/OPX ratio than the older basalts (Im1). The kink (deflection diameter) between the youngest (3.38 Ga) and oldest (3.72) isochron curves is ~1 km, suggesting the thickness of buried lava flows to be >100 m [START_REF] Ling | Composition, mineralogy and chronology of mare basalts and non-mare materials in Von Kármán crater: Landing site of the Chang'E-4 mission[END_REF]. The estimated total basalt thickness (≥320 m) in this study agrees with this estimate. As a result, a refined stratigraphic sequence for the low-Ti mare basalts [START_REF] Ling | Composition, mineralogy and chronology of mare basalts and non-mare materials in Von Kármán crater: Landing site of the Chang'E-4 mission[END_REF][START_REF] Qiao | Geological characterization of the Chang'e-4 landing area on the lunar farside[END_REF][START_REF] Fu | The subsurface structure and stratigraphy of the Chang'E-4 landing site: orbital evidence from small craters on the Von Kármán crater floor[END_REF] is proposed in this study, from top to bottom: Im3 (~107 m), Im2 (~84 m), and Im1 (≥129 m).

A Lunar Penetrating Radar (LPR), with two work frequency channels (channel 1, 60 MHz; channel 2, 500 MHz), is part of the instrumentation on the Yutu-2 rover to study the subsurface structure of its traversing area [START_REF] Fang | Lunar Penetrating Radar onboard the Chang'e-3 mission[END_REF][START_REF] Jia | The scientific objectives and payloads of Chang'E-4 mission[END_REF].

Initial interpretations of the data indicate that there are several subsurface layers/units, pointing to a complex thermal history. For example, the high-frequency LPR data reveal that the shallow subsurface at the landing area is roughly composed of three layers [START_REF] Lai | Comparison of Dielectric Properties and Structure of Lunar Regolith at Chang'e-3 and Chang'e-4 Landing Sites Revealed by Ground Penetrating Radar[END_REF][START_REF] Li | The Moon's farside shallow subsurface structure unveiled by Chang'E-4 Lunar Penetrating Radar[END_REF]Zhang et al., 2020a;Zhang et al., 2020b): a topmost finegrained regolith layer (depth, 0 to 12 m), an intermediate layer consisting of coarse materials with embedded rocks (12 to 24 m), and a bottom layer composed of alternating coarse and fine materials (24 to 40 m). As of November 22, 2020, the Yutu-2 rover had already operated for 24 lunar days and traversed about 590 m on the Finsen ejecta unit within the Von Kármán crater. The rover is currently about 30 km away from the Zhinyu crater. The collected low-frequency LPR channel has the maximum detection depth of ~330m [START_REF] Lai | First look by the Yutu-2 rover at the deep subsurface structure at the lunar farside[END_REF], possibly up to ~500 m (Zhang et al., 2020a).

However, due to the low thickness resolution (~10 m) of the LPR low-frequency data and different processing procedures, the deep multi-layered stratigraphic structure revealed by the low-frequency LPR data is not entirely consistent between most of the studies, including [START_REF] Lai | First look by the Yutu-2 rover at the deep subsurface structure at the lunar farside[END_REF], Zhang et al. (2020a), andZhang et al. (2020b).

Nevertheless, the basaltic stratigraphic sequence at the Zhinyu crater proposed in this study roughly corresponds to the D and E layers in the work of [START_REF] Lai | First look by the Yutu-2 rover at the deep subsurface structure at the lunar farside[END_REF]. The rover will continue to traverse westward and it is planned to reach the mare basalt area in the future. The proposed stratigraphic sequence and basalt thickness are expected to be validated by the LPR data, which will provide ground truth and should contribute to a fuller understanding of the basalt thermal history of the Von Kármán crater.

Implications for magma evolution

The ejecta under investigation is composed of materials exhumed by the Zhinyu crater during formation, a process that inverts the target's stratigraphical sequence.

Therefore, the variation of ejecta composition along the radial distances around Zhinyu crater reflects the compositional variation of the basalt source at different stratigraphic depths. A fairly homogenous OL composition (Fo#: 55) of the three layers indicates the sources of these mare basalts (≥320 m in thickness) are almost homogenous in composition, or might even be from a single magma chamber. Generally, for primary magmas, forsteritic (high Fo#) OL is indicative of a primitive source, and more fayalitic (low Fo#) OL is indicative of an evolved source (Basaltic Volcanism Study Project, 1981). The predicted Fo# value points to a moderate magma evolution. The small mineral abundance differences (Fig. 3) among layers may indicate there were at least three mare flooding events occurred inside Von Kármán crater during the Imbrian period, with two dome-like structures being the possible latest volcanic vents of these basalts [START_REF] Qiao | Geological characterization of the Chang'e-4 landing area on the lunar farside[END_REF].

Conclusion

China's Chang'e-4 probe landed inside the Von Kármán crater, whose floor was flooded with mare basalts during the Imbrian period. Zhinyu, a fresh impact crater, is the largest close to the landing site. The materials exhumed by the blast would have distributed the ejecta around the crater, exposing the stratigraphy sequence of the underlying basaltic rock layer(s). The mineral abundances derived from multispectral data, vary radially and discontinuously from the rim craters outwards, forming compositionally distinguished annuli. In particular, variations in the mineral olivine and clinopyroxene abundances support the idea that there could be at least three distinct basalt layers (named Im1, Im2, and Im3 from bottom to top) within the Von Kármán crater. The thickness of each layer is estimated by an empirical relation according to the processes of impact excavation and ejecta deposition, i.e., the shallow material deposits at the further distance whereas the deep material deposits around the crater. As it cannot be determined in this study whether the Zhinyu impact penetrated through the entire Im1 layer, the total thickness of the three basalt layers is estimated to be at least 320 m.

However, the minor differences noted in the mineral abundances and the fairly consistent olivine compositions (Fo#: 55) of these three layers indicate the sources of these mare basalts to be comparable in composition, perhaps even from a single magma chamber where the magma underwent only moderate evolution. The westward traversing Yutu-2 rover may provide ground truth for the findings in this study, which would contribute to the understanding of the thermal history of the Von Kármán crater.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this article.

Fig. 1 .

 1 Fig. 1. Location of the Zhinyu crater inside the Von Kármán crater. The base map

Fig. 2 .

 2 Fig. 2. Three mineral abundance variation rings distributed concentrically around

  distance (Dmax) comes from the surface, and the radial position of the continuous ejecta blanket is inversely related to the excavation depth, i.e., (Hexc -Sce)/Hexc = D/Dmax. Therefore, the source depth of continuous ejecta blanket (Sce) at a given radial position (D) is estimated by equation Sce = Hexc -(D*Hexc/Dmax).

Fig. 3 .

 3 Fig. 3. Mineral abundances of the identified three basalt layers. The error bars

Fig. 4 .

 4 Fig. 4. FeO and TiO2 contents of identified three mare basalt layers.

Acknowledgments

The M 3 data are from PDS Geoscience Node (https://pdsgeosciences.wustl.edu/). The MI data are from Selene Data Archive (http://l2db.selene.darts.isas.jaxa.jp/). The MI-derived mineral maps are from (https://astrogeology.usgs.gov/maps/lunar-kaguya-multiband-imager-mosaics). This work was supported by the Strategic Priority Research Program of Chinese Academy of Sciences (grant No. XDB41000000), National Natural Science Foundation of China (grant Nos. 41702354 and 41941003), Macao Young Scholars Program (grant No. AM201902), and Science and Technology Development Fund of Macau (grant No. 131/2017/A3).

Supplementary material

Supplementary material associated with this article can be found in the online version.

Supplementary Material

Click here to access/download e-Component Mare basalt flooding revealed by Zhinyu crater_SM_1208.docx