domain: spi phys

text:
Keywords: Anisotropy, Geant4, Monte-Carlo simulations, proton, radiation belts, radiation monitor, response function, space environment

The ICARE-NG instrument onboard the Argentinian satellite SAC-D detected much more protons during descending orbits (when latitude decreases) than for ascending orbits (increasing latitudes). In this paper we will investigate on the anisotropy seen by the ICARE-NG Head-A for protons in coincidence mode from Monte-Carlo simulations performed with GEANT4. Our simulations show that the difference in the fluxes observed during ascending and descending orbits comes from the fact that the instrument observed trapped protons or not on each point of the orbits as a result of the instrument and satellite orientations. In addition, we show in this paper that the measurements performed by ICARE-NG can be used in conjunction with our GEANT4 simulations to study the anisotropy of trapped protons, i.e. their distribution relative to their equatorial pitch-angle.

I. INTRODUCTION

INCE their discovery in 1958 by Explorer and Pioneer probes, the Earth's radiation belts are widely studied especially because kinetic energies of particles are very high and may cause damages on human bodies and on spacecraft, sometimes leading to their loss [START_REF] Baker | The occurrence of operational anomalies in spacecraft and their relationship to space weather[END_REF] [START_REF] Townsend | Implications of the space radiation environment for human exploration in deep space[END_REF]. In order to better understand radiation belts, monitors are implemented onboard satellites to give more information about particles trapped in radiation belts. Several studies have been done to investigate anisotropy in the inner belt. In particular, Siegl [START_REF] Siegl | Inner belt anisotropy investigations based on the standard radiation environment monitor (SREM)[END_REF] uses a combination of the AP8 omnidirectional flux model [START_REF] Vette | The NASA/National Space Science Data Center Trapped Radiation Environment Model Program (1964-1991[END_REF] [5] with the Badhwar-Konradi model [START_REF] Badhwar | Conversion of omnidirectional proton fluxes into a pitch angle distribution[END_REF] to investigate anisotropy for SREM [START_REF] Mohammadzadeh | The ESA standard radiation environment monitor program first results from PROBA-I and INTEGRAL[END_REF] coincidence channels. The shape parameter of the Badhwar-Konradi model is calculated by fitting the distribution of the flux given by the model to the AP8 fluxes of 20MeV trapped protons. Based on this model we will calculate the shape parameter out of fluxes recorded by the ICARE-NG Head-A (and not using AP8) for 18.6MeV trapped protons to compare with results obtained by Siegl, and Manuscript submitted December 20, 2018. This work was conducted as part of an ONERA internship and an ONERA / CNES thesis.

M. Ruffenach, S. Bourdarie and D. Falguère are with ONERA The French aerospace lab / Département Physique Instrumentation Environnement et Espace, 2 Av. Edouard Belin, FR-31055 Toulouse Cedex 4, France (telephone: +33 5 62 25 25 68, e-mail: Marine.Ruffenach@onera.fr).

J.R. Vaillé is with the Université de Montpellier, IES -UMR UM/CNRS 5214, 860 Rue de St Priest, Bat. 5, F-34097 Montpellier, France.

J. Mekki is with CNES, 18 av. E. Belin 31401 Toulouse, France.

for 81MeV trapped protons to evaluate the energy dependence of the anisotropy. Other models will be used in this paper to do more comparisons.

In section II, a brief presentation of the ICARE-NG instrument will be given. In section III conditions of GEANT4 [START_REF]Geant4: A toolkit for the simulation of the passage of particles through matter[END_REF] [9] [START_REF] Allison | Geant4 developments and applications[END_REF] simulations are presented with the description of the particles' source and the establishment of the protons response function in coincidence mode of ICARE-NG Head-A. Section IV focuses on GEANT4 simulations for the monitor on-board the SAC-D satellite. In particular, matrices providing the number of particles which can be measured according to their equatorial pitch-angle will be established for several satellite's locations and orbit type (ascending and descending orbit). Finally, in section V we will establish a model describing the anisotropy for 18.6MeV and 81MeV energies protons, respectively the lowest and the highest energies not contaminated by electrons in coincidence mode for the ICARE-NG Head-A.

II. ICARE-NG INSTRUMENT

This paper focuses on the CNES ICARE-NG instrument onboard the low Earth orbiting Argentinian satellite SAC-D which flew between June 10, 2011 and June 8, 2015. The complete description of the satellite and instrument can be found in [START_REF] Boscher | In Flight Measurements of Radiation Environment on Board the French Satellite JASON-2[END_REF] and in [START_REF] Boscher | In-Flight Measurements of Radiation Environment on Board the Argentinean Satellite SAC-D[END_REF]. The satellite was on a sunsynchronous orbit at a mean altitude of 650 km and an orbit inclination of 98.01°. The ICARE-NG instrument consists of three detector heads (A, B, and C). Both sensors A and C have a pair of 500 µm thick silicon diodes such that a coincidence and an anti-coincidence acquisition mode are available. These heads are set to measure protons and electrons. The Head-B is composed of a single 700 µm thick silicon diode to measure electrons. Particle fluxes for electrons from 250keV to 3.2MeV and for protons from 12.8MeV to 190MeV are obtained. SAC-D was a spin-stabilized satellite, thus ICARE-NG always pointed to the same direction which was 45° from the zenith and from the satellite velocity vector (see [START_REF] Boscher | In-Flight Measurements of Radiation Environment on Board the Argentinean Satellite SAC-D[END_REF] and section IV for more details). When the satellite was flying on descending orbits (decreasing latitudes), the ICARE-NG monitor recorded fluxes higher than the ones for ascending orbits (increasing latitudes). An illustration of these measurements is shown in Fig. 1 for protons of 65MeV detected in coincidence mode with the Head-A during January, 2015. This is attributed to the anisotropy of the environment. In this paper we focus on the Head-A of ICARE-NG in coincidence mode for which deposited energy is discretized in 128 channels from 0MeV to 9.51MeV. To improve the signal to noise ratio original channels may be summed up when level 2 data are produced (Table I).

Proton radiation belts anisotropy as seen by ICARE-NG

III. GEANT4 SIMULATIONS FOR THE ICARE-NG HEAD-A

Response function of protons in coincidence mode of the Head-A is computing from a Monte-Carlo simulation based on the GEANT4 toolkit. The geometry of the Head-A was drawn and translated into GEANT-4 compatible mesh using FASTRAD Software [START_REF] Thomas | FASTRAD: A 3D CAD Interface for Radiation Calculation and Shielding[END_REF] [START_REF]FASTRAD Software -The 3D CAD Tool For Radiation Shielding Analysis[END_REF]. In the simulations described in sections III and IV, incident energies are defined from 14MeV to 150MeV by step of 2MeV to cover the full energy range seen by the Head-A. When very good statistical results are required, the Monte-Carlo simulation is performed with 10 million tries otherwise two million tries are assumed to be good enough. In these simulations a particle source with an isotropic distribution is assumed, thus we consider a spherical source with a radius equals to 3.9𝑐𝑐𝑐𝑐 and a cosine-law angular distribution for velocity vector directions [START_REF] Santin | Normalisation modelling sources[END_REF]. With this source, particles come from all around the Head-A.

To have information about privileged detected particles taking into account the shielding of the Head-A surrounding diodes we calculate the geometric factor. The geometric factor 𝐺𝐺𝐺𝐺𝐺𝐺 is obtained from the ratio between the number of particles detected per channel over the total number of particles launched by incident energy, expression given in (1) and calculated using [START_REF] Santin | Normalisation modelling sources[END_REF] and [START_REF] Sullivan | Geometric factor and directional response of single and multi-element particle telescopes[END_REF].

𝐺𝐺𝐺𝐺𝐺𝐺 = 4𝜋𝜋 2 𝑅𝑅 2 𝑁𝑁 𝑑𝑑𝑑𝑑𝑑𝑑 𝑁𝑁 𝑇𝑇𝑇𝑇𝑑𝑑𝑇𝑇𝑇𝑇 ( 1
)
where 𝑅𝑅 is the radius of the source.

Response function to protons in coincidence mode is shown on Fig. 2. In this response function two branches can be seen, the first one is attributed to particles coming from the main entrance of the detector, and the second one is attributed to penetrating particles crossing the shielding. This is the response function for protons in coincidence mode for particles coming from all around the Head-A without taking into account the satellite shielding and any anisotropy of the proton environment. In the following section we will take into consideration the location and orientation of the satellite to calculate which particles' pitch-angle are reaching the detectors.

IV. THE INFLUENCE OF THE SATELLITE'S POSITION ON PRIVILEGED DIRECTIONS OF INCOMING PARTICLES

Since the anisotropy is expected to be more pronounced in coincidence mode than in anticoincidence mode we focus here on the analysis of the response function in coincidence mode of the Head-A. Indeed, to be able to see the impact of the anisotropy we must have a narrow field of view otherwise any anisotropy would be smoothed out because of the angle averaging for a large field of view (in the case of anticoincidence for example).

According to the conservation of the first adiabatic invariant along a drift shell we have:

1 𝐵𝐵 𝑚𝑚 = 𝑠𝑠𝑠𝑠𝑠𝑠²(𝛼𝛼 𝑑𝑑𝑒𝑒 ) 𝐵𝐵 𝑑𝑑𝑒𝑒 = 𝑠𝑠𝑠𝑠𝑠𝑠²(𝛼𝛼 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 ) 𝐵𝐵 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 (2)
with 𝛼𝛼 𝑑𝑑𝑒𝑒 the equatorial pitch-angle and 𝛼𝛼 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 the pitch-angle of the particle at the considered position. 𝐵𝐵 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 , 𝐵𝐵 𝑚𝑚 and 𝐵𝐵 𝑑𝑑𝑒𝑒 are respectively the magnetic field magnitude at the considered point, the magnetic field magnitude at the mirror point, and the magnetic field amplitude at the magnetic equator. We will use this equation in particular to convert local pitch-angles into equatorial pitch-angles.

In order to study the anisotropy it is necessary to know particles velocity vector information and the magnetic field vector for each location along the spacecraft orbit. These values have to be expressed in the same frame to calculate the pitch-angle. Particles coordinates and their velocity vector provided by GEANT4 are expressed in the ICARE-NG coordinate system while the magnetic field vector coordinates are expressed in the geographic system. Knowing the position of ICARE-NG on the satellite [START_REF] Trad | SAC-D CARMEN -Analyse sectorielle[END_REF], and knowing the spacecraft orientation with respect to the Earth, it is then possible to express particle coordinates in the geographic coordinate system to compute the local pitch-angle of incident particle called 𝛼𝛼 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 , impinging the Head-A of ICARE-NG at each orbit location.

The ICARE-NG sensor is part of CARMEN equipment [START_REF] Bezerra | CARMEN/MEX test board for the study of radiation effects on electronic components aboard JASON-2 and SAC-D satellites[END_REF] on-board SAC-D. The location of ICARE-NG on the satellite is illustrated on Fig. 3. The origins of the two coordinate systems are separated by a few dozen centimeters so we can neglect Earth magnetic field variations on such distance. In this respect we assume that the CARMEN coordinate system is equivalent to that of the Head-A. Then, the geometry of the SAC-D satellite and the orientation of the ICARE-NG instrument according to the satellite axis are shown in Fig. 4.

To convert from the detector coordinate system to the satellite coordinate system we need to rotate the coordinates by an angle of -3𝜋𝜋/4 around the Y-axis. So, the relation to express a vector in the satellite coordinate system from the detector coordinate system is given by (3).

� 𝑥𝑥

𝑠𝑠𝑇𝑇𝑑𝑑 𝑦𝑦 𝑠𝑠𝑇𝑇𝑑𝑑 𝑧𝑧 𝑠𝑠𝑇𝑇𝑑𝑑 � = ⎝ ⎜ ⎛ - 1 √2 0 1 √2 0 1 0 - 1 √2 0 - 1 √2 ⎠ ⎟ ⎞ � 𝑥𝑥 𝑑𝑑𝑑𝑑𝑑𝑑 𝑦𝑦 𝑑𝑑𝑑𝑑𝑑𝑑 𝑧𝑧 𝑑𝑑𝑑𝑑𝑑𝑑 � (3)
Now, we have to convert satellite coordinates into geographic coordinates. The velocity of the satellite being along the X-axis of the satellite coordinate system, the Z-axis of the satellite always points toward the Earth and the Y-axis completes this trihedral. For all the satellite locations along its orbit, its geodesic coordinates (longitude, latitude, altitude) are known and can be converted into geographic coordinates. From these geographic coordinates we can calculate the velocity of the satellite using two consecutive points of the orbit and the Earth-Satellite direction in the geographic system. After normalization, it gives us the conversion of the X-axis and the Z-axis from the satellite coordinate system to the geographic system. The Y-axis projected in the geographic coordinate system completes this trihedral. Having the velocity of particles 𝑉𝑉 𝐺𝐺𝐺𝐺𝐺𝐺 ���������⃗ and the magnetic field 𝐵𝐵 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 �����������⃗ in the same coordinate system, it is straightforward to calculate the local pitch-angle using the scalar product given in [START_REF] Vette | The NASA/National Space Science Data Center Trapped Radiation Environment Model Program (1964-1991[END_REF].

𝐵𝐵 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 �����������⃗ . 𝑉𝑉 𝐺𝐺𝐺𝐺𝐺𝐺 ���������⃗ = �𝐵𝐵 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 �����������⃗ �. �𝑉𝑉 𝐺𝐺𝐺𝐺𝐺𝐺 ���������⃗ �. cos (𝛼𝛼 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 ) (4)
Due to the Earth magnetic field properties, eccentric, tilted field with respect to geographic coordinates, SAC-D will encounter a given field line, at different 𝐵𝐵 𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙 𝐵𝐵 𝑒𝑒𝑒𝑒 at different longitudes. This is to say that the maximum equatorial pitchangle ICARE-NG can measure on a given field line depends on the longitude at which it will cross this given field line. The maximum equatorial pitch-angle SAC-D/ICARE-NG can measure as a function of the longitude is illustrated on Fig. 5 for L=1.185 at an altitude of 664km. These data are extracted from the orbitography of the SAC-D satellite.

In order to assess the anisotropy of trapped protons on a given field line, three spacecraft locations are considered (red crosses on Fig. 5). For each position, two configurations are possible: ascending orbit and descending orbit.

The equatorial pitch-angle which defines the loss-cone is such that the particle has a mirror point at an altitude of 100km during its bounce and drift along a given drift shell. Thanks to the equation ( 2) only the minimal magnetic field value obtained at 100km during a drift (limit point for the mirror point beyond which particles precipitate in the upper atmosphere), and the magnetic field value at the magnetic equator are necessary to calculate the loss-cone equatorial pitch-angle. These values are obtained using IRBEM library [START_REF] Bourdarie | The International Radiation Belt Environment Modeling (IRBEM) library[END_REF]. This open-source library allows to calculate magnetic coordinates in the Earth's magnetic field and to perform coordinates conversions. For L=1.185 the equatorial pitchangle of the loss-cone is equal to 66.27°. Fig. 6 represents a magnetic field line with a satellite on it. Fig. 7 gives the number of protons reaching the two diodes (among the 2 million particles per bin of 2MeV being launched) in an incident energy -equatorial pitch-angle map for the three locations represented on Fig. 5 and for ascending and descending orbits. In Fig. 7, the horizontal red line (also represented in red in Fig. 6) indicates the loss cone equatorial pitch-angle (66.27° at L=1.185). Particles with a smaller pitchangle are not trapped and precipitate into the atmosphere. Fig. 6. The satellite on a magnetic field line. The red line represents the limit of the atmosphere where particles precipitate. Hatched area represents particles which are not seen by the satellite since their mirror point is between the satellite and the magnetic equator. Fig. 7. Number of simulated protons in coincidence mode in an incident energy (X-axis) -equatorial pitch-angle (Y-axis) map. This is for a field line with L = 1.185 and an altitude equals to 664km. Ascending orbits (top panel) and descending orbits (bottom panel) are classified according to increasing longitudes from left to right. Hatched areas represent particles that cannot be observed because of satellite's position. For each bin of incident energy, two million of particles are simulated. Fig. 7 does represent what the Head-A has actually observed, but it points out what might be the pitch-angles of observed particles. If no particle is observed by the detector, it means that the instrument has no access to trapped particles because of its orientation. The instrument has no access to particles with an equatorial pitch-angle higher than 55° when the spacecraft is located in the -165° longitude range according to matrices on the left, 73° when the spacecraft is located in the -90° longitude range according to matrices on the center, while all equatorial pitch-angles can be observed when the spacecraft is located in the -63° longitude range according to matrices in the right panels. These equatorial pitch-angle values are the maximum observable equatorial pitch-angle (cf. Fig. 5), beyond which the sensor does not detect particles regardless of the loss cone. Indeed, this is related to the location of the satellite along a given field line: when it has a magnetic latitude high enough then particles with a high equatorial pitch-angle have a mirror point between the magnetic equator and the satellite (hatched area on Fig. 6). Hatched areas on Fig. 7 represent particles that cannot be observed because of the satellite position along a given field line.

Fig. 8. Evolution of the number of simulated trapped protons summed over all the incident energies (from 14MeV to 150MeV) as a function of longitude for ascending orbits (in red) and descending orbits (in blue).

Fig. 8 gives the number of simulated trapped protons summed over all the incident energies (total number of trapped protons summed from 14MeV to 150MeV) as a function of the spacecraft longitude. As a reminder, Fig. 7 and Fig. 8 are obtained thanks to simulations.

Considering maps at the left in Fig. 7, i.e. when the spacecraft is at a longitude lower than -95°, the sensor has only access to protons having an equatorial pitch-angle lower than the one defining the loss-cone i.e. all these protons are precipitating into the upper atmosphere. As a result (see Fig. 8), when the satellite is at a longitude lower than -95°, the Head-A has no access to trapped particles so it cannot detect trapped protons at these longitudes. For longitudes greater than -95° we can see on Fig. 8 that more and more trapped protons can be measured. Also note that more trapped protons are counted in the case of descending orbits. Maps at the center and on the right of Fig. 7 indicate that the sensor has access to particles with an equatorial pitch-angle higher than the loss cone pitch-angle (66.27°). From the middle maps to the right maps of Fig. 7, the number of trapped particle the sensor has access to is increasing. Also, from top to bottom (middle and right maps of Fig. 7), i.e. increasing or decreasing orbits, the number of trapped protons the sensor has access to is increasing. This is consistent with the results shown in Fig. 8 for longitude higher than -95°.

To summarize, Fig. 7 and Fig. 8 are obtained from GEANT4 simulations. First, we have considered an isotropic environment at the outside of the Head-A. Then, for each location of SAC-D along its orbit we can calculate particles' equatorial pitch-angle impinging the diodes in coincidence mode and find out if the sensor has access to trapped particles or not.

In situ measurements from coincidence channels, like the ones obtained thanks to the ICARE-NG Head-A illustrated in Fig. 1, are well understood with these GEANT4 simulations. They explain qualitatively the difference in the number of detected particles between ascending and descending orbits by highlighting the access to trapped protons or not versus spacecraft location. To fully characterize this anisotropy, in the next section we will perform a comparison between different models to describe the evolution of the unidirectional flux according to the equatorial pitch-angle and the McIlwain L parameter.

V. ESTIMATING UNIDIRECTIONAL FLUXES

In this section we will express the unidirectional flux as a function of equatorial pitch-angle according to ICARE-NG measurements during 2013 and local pitch-angle -energy map of this detector. We compare four models of unidirectional fluxes expressed in [𝑐𝑐𝑐𝑐 -2 𝑠𝑠 -1 𝑠𝑠𝑠𝑠 -1 𝑀𝑀𝑀𝑀𝑉𝑉 -1 ]:

1. Badhwar-Konradi [START_REF] Badhwar | Conversion of omnidirectional proton fluxes into a pitch angle distribution[END_REF] using AP8 model 2. Valot [START_REF] Valot | Pitch Angle Distribution of Geomagnetically Trapped Protons for 1[END_REF] using ICARE-NG data 3. Badhwar-Konradi using ICARE-NG data 4. Roederer [START_REF] Roederer | Dynamics of Geomagnetically Trapped Radiation[END_REF] For the second and the third models parameters have to be estimated using in-flight data of ICARE-NG Head-A by comparing counts of particles calculated with these two models from those measured by ICARE-NG. Counts of particles 𝐶𝐶 are calculated thanks to (5):

𝐶𝐶 = � � 𝑓𝑓𝑓𝑓𝑓𝑓𝑥𝑥(𝐺𝐺, 𝛼𝛼 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 ) * 𝑀𝑀𝑀𝑀𝑀𝑀(𝐺𝐺, 𝛼𝛼 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 ) * 𝑑𝑑𝐺𝐺𝑑𝑑𝛼𝛼 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 𝐺𝐺 𝛼𝛼 𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙 (5)
where 𝑀𝑀𝑀𝑀𝑀𝑀(𝐺𝐺, 𝛼𝛼 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 ) is the geometric factor at a given energy and local pitch-angle. 𝐺𝐺𝑓𝑓𝑓𝑓𝑥𝑥(𝐺𝐺, 𝛼𝛼 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 ) is the unidirectional flux as function of incident energy and local pitch-angle. Unidirectional fluxes for these models are calculated for 𝐿𝐿 = 1.25 and for two energies, 18.6𝑀𝑀𝑀𝑀𝑉𝑉 and 81𝑀𝑀𝑀𝑀𝑉𝑉.

A. Model 1: Badhwar-Konradi using AP8 model

This model is the one proposed by NASA to study fluxes at the International Space Station altitude. This distribution is given by ( 6) for equatorial pitch-angles which are between the angle defining the and 𝜋𝜋 2 , otherwise the flux is equal to 0 (in the loss-cone).

𝑗𝑗 𝐵𝐵𝐵𝐵-𝐴𝐴𝐴𝐴8 = 𝐾𝐾𝐾𝐾𝑀𝑀 -𝛽𝛽𝛽𝛽

where 𝐾𝐾 = 𝑠𝑠𝑠𝑠𝑠𝑠�𝛼𝛼 𝑒𝑒𝑒𝑒 � � 𝐵𝐵 𝑒𝑒𝑒𝑒 -1 �𝐵𝐵𝑙𝑙𝑙𝑙 [START_REF] Badhwar | Conversion of omnidirectional proton fluxes into a pitch angle distribution[END_REF] with 𝐵𝐵 𝑇𝑇𝑙𝑙 the magnetic field magnitude for a particle having its mirror point on the Earth atmosphere. Magnetic fields magnitudes are expressed in 𝐺𝐺𝑀𝑀𝑓𝑓𝑠𝑠𝑠𝑠. The 𝐾𝐾 parameter is a normalization factor calculated so that the integral of differential flux with respect to pitchangles is equal to the omnidirectional flux, but in our case we do not need to calculate this parameter since we will compare hereafter normalized unidirectional fluxes. Then, an expression of the 𝛽𝛽 shape parameter is given in [START_REF] Kruglanski | Trapped Proton Anisotropy at Low Altitudes[END_REF] and used by Siegl in [START_REF] Siegl | Standard radiation environment monitor : simulation and inner belt flux anisotropy investigation[END_REF]. The shape parameter is calculated for 20MeV trapped protons using the AP8 model as a function of the McIlwain L parameter as shown in [START_REF] Mohammadzadeh | The ESA standard radiation environment monitor program first results from PROBA-I and INTEGRAL[END_REF].

𝛽𝛽�√𝐺𝐺� = 1 0.13164 -8.8674ln (𝐿𝐿) (7)

B. Model 2: Valot model using ICARE-NG data

The Valot model is given in [START_REF]Geant4: A toolkit for the simulation of the passage of particles through matter[END_REF] where 𝑠𝑠 and 𝑗𝑗 0 have to be estimated, and where 𝛼𝛼 𝑇𝑇𝑙𝑙.𝑑𝑑𝑒𝑒 is the equatorial pitch-angle defining the loss-cone (mirror point at 100km of altitude in the upper atmosphere).

𝑗𝑗 𝑉𝑉𝑇𝑇𝑇𝑇𝑇𝑇𝑑𝑑 = 𝑗𝑗 0 �𝑠𝑠𝑠𝑠𝑠𝑠�𝛼𝛼 𝑑𝑑𝑒𝑒 � 𝑠𝑠 -𝑠𝑠𝑠𝑠𝑠𝑠�𝛼𝛼 𝑙𝑙𝑐𝑐.𝑑𝑑𝑒𝑒 � 𝑠𝑠 � (8)
Count of particles are calculated using local pitch-angles while the Valot model uses equatorial pitch-angles. Moreover, Monte-Carlo simulations gives us the number of particles in an incident energy -local pitch-angle map that we call 𝑠𝑠𝑓𝑓𝑐𝑐𝑛𝑛𝑀𝑀𝑠𝑠(𝐺𝐺, 𝛼𝛼 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 ). Finally count of particles, expressed in particle per second, is calculated with [START_REF] Agostinelli | Geant4-a simulation toolkit[END_REF]. [START_REF] Agostinelli | Geant4-a simulation toolkit[END_REF] where 𝑠𝑠𝑠𝑠𝑠𝑠𝑖𝑖𝑠𝑠𝑀𝑀𝑓𝑓 𝑠𝑠𝑓𝑓𝑐𝑐𝑛𝑛𝑀𝑀𝑠𝑠 is the number of particles launched during the Monte-Carlo simulation for a given incident bin in energy. For these simulations we use an 𝑠𝑠𝑠𝑠𝑠𝑠𝑖𝑖𝑠𝑠𝑀𝑀𝑓𝑓 𝑠𝑠𝑓𝑓𝑐𝑐𝑛𝑛𝑀𝑀𝑠𝑠 equals to ten million of protons per bin of 1MeV from 14 to 150MeV to get high fidelity statistical results.

𝐶𝐶 = � �
In order to have several points of comparison, we consider all positions of the SAC-D satellite corresponding to the field line with L=1.25 at an altitude of 664km during 2013. We focus on 18.6MeV protons, so channels with numbers ranging from 50 to 78 for ICARE-NG Head-A, the lowest proton energy for the coincidence mode. For a given channel, i.e. for a given measured incident energy, the count rate is given by [START_REF] Allison | Geant4 developments and applications[END_REF].

𝐶𝐶 𝑚𝑚𝑑𝑑𝑇𝑇𝑠𝑠𝑚𝑚𝑚𝑚𝑑𝑑𝑑𝑑 = 𝑓𝑓𝑓𝑓𝑓𝑓𝑥𝑥 𝑚𝑚𝑑𝑑𝑇𝑇𝑠𝑠𝑚𝑚𝑚𝑚𝑑𝑑𝑑𝑑 * 𝐺𝐺𝐺𝐺𝐺𝐺 * ∆𝐺𝐺 (10)
The 𝐺𝐺𝐺𝐺𝐺𝐺 * ∆𝐺𝐺 factor for 18.6MeV protons is 0.7𝑐𝑐𝑐𝑐 2 . 𝑠𝑠𝑠𝑠. 𝑀𝑀𝑀𝑀𝑉𝑉. We have to find 𝑠𝑠 and 𝑗𝑗 0 parameters so that counts calculated thanks to [START_REF] Agostinelli | Geant4-a simulation toolkit[END_REF] are equal to those measured by ICARE-NG given by [START_REF] Allison | Geant4 developments and applications[END_REF]. We plot counts measured by ICARE-NG as a function of counts calculated thanks to [START_REF] Agostinelli | Geant4-a simulation toolkit[END_REF], and the objective is to modify 𝑠𝑠 and 𝑗𝑗 0 until we obtain a linear plot with a slope close to 1 (equality between counts measured and calculated). We know that for a scatter plot (x, y), the estimation by least squares of the slope of the linear function that best fits these points is given by 𝑠𝑠𝑓𝑓𝑠𝑠𝑀𝑀𝑀𝑀 =

∑ 𝑥𝑥 𝑖𝑖 𝑦𝑦 𝑖𝑖 𝑖𝑖 ∑ 𝑥𝑥 𝑖𝑖 2 𝑖𝑖
and the estimation of the variance of the linear regression is given by

𝑆𝑆 𝑥𝑥𝑦𝑦 2 = ∑ 𝑦𝑦 𝑖𝑖 2 -𝑠𝑠𝑇𝑇𝑇𝑇𝑐𝑐𝑑𝑑 2 ∑ 𝑥𝑥 𝑖𝑖 2 𝑖𝑖 𝑖𝑖 𝑠𝑠 𝑝𝑝 -1
with 𝑠𝑠 𝑐𝑐 the number of points. We proceed by dichotomy to find the 𝑠𝑠 parameter which minimizes 𝑆𝑆 𝑥𝑥𝑦𝑦 2 and then we adjust the 𝑗𝑗 0 parameter to have a slope equals to 1. Values which give us the best match with measurements of ICARE-NG are 𝑗𝑗 0 = 54400𝑐𝑐𝑐𝑐 -2 𝑠𝑠 -1 𝑠𝑠𝑠𝑠 -1 𝑀𝑀𝑀𝑀𝑉𝑉 -1 and 𝑠𝑠 = 45. The result of the minimization process for E=18.6MeV is shown in Fig. 9 in blue.

For ICARE-NG Head-A, the highest energy for protons for the coincidence mode not contaminated by electrons is 81MeV, corresponding to the channel 12 for which the 𝐺𝐺𝐺𝐺𝐺𝐺 * ∆𝐺𝐺 factor is equal to 4.37𝑐𝑐𝑐𝑐 2 . 𝑠𝑠𝑠𝑠. 𝑀𝑀𝑀𝑀𝑉𝑉. We use the same technic than the one used for the 18.6MeV energy to find the estimated parameters and we obtain 𝑠𝑠 = 24, 𝑗𝑗 0 = 15800𝑐𝑐𝑐𝑐 -2 𝑠𝑠 -1 𝑠𝑠𝑠𝑠 -1 𝑀𝑀𝑀𝑀𝑉𝑉 -1 . The result of the minimization process for E=81MeV is shown in Fig. 10 in blue. Those results indicate that the anisotropy is energy dependent.

C. Model 3: Badhwar-Konradi model using ICARE-NG data

The initial formula is the same than the one used for the Badhwar-Konradi model based on AP8.

𝑗𝑗 𝐵𝐵𝐵𝐵-𝐼𝐼𝐼𝐼𝐴𝐴𝐼𝐼𝐺𝐺 = 𝑘𝑘𝐾𝐾𝑀𝑀 -𝐵𝐵𝛽𝛽 (11)
Instead of using the AP8 model to calculate B we estimate 𝑘𝑘 and 𝐵𝐵 parameters using the same method as used above. We obtain 𝑘𝑘 = 725𝑐𝑐𝑐𝑐 -2 𝑠𝑠 -1 𝑠𝑠𝑠𝑠 -1 𝑀𝑀𝑀𝑀𝑉𝑉 -1 𝐺𝐺 0.5 and 𝐵𝐵 = -10𝐺𝐺 0.5 for an energy of 18.6MeV, and 𝑘𝑘 = 1675𝑐𝑐𝑐𝑐 -2 𝑠𝑠 -1 𝑠𝑠𝑠𝑠 -1 𝑀𝑀𝑀𝑀𝑉𝑉 -1 𝐺𝐺 0.5 and B = -8𝐺𝐺 0.5 for an energy of 81MeV. The results of the minimization process for 18.6MeV and for 81MeV are shown in Fig. 9 and in Fig. 10 in red. Again those results indicate that the anisotropy is energy dependent.

D. Model 4: Roederer using AP8 model

The Roederer model consists in calculating unidirectional flux from omnidirectional flux given by the AP8 model. The unidirectional flux expression perpendicular to the magnetic field extracted from [START_REF] Roederer | Dynamics of Geomagnetically Trapped Radiation[END_REF] is given by [START_REF] Boscher | In-Flight Measurements of Radiation Environment on Board the Argentinean Satellite SAC-D[END_REF].

𝑗𝑗 ⊥ (𝐵𝐵 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 ) = - 𝐵𝐵 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 3/2 2𝜋𝜋 2 � 𝑑𝑑 𝑑𝑑𝐵𝐵 � 𝐽𝐽(𝐵𝐵) 𝐵𝐵 � 𝑑𝑑𝐵𝐵 �𝐵𝐵 -𝐵𝐵 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 𝐵𝐵 𝑙𝑙𝑙𝑙 𝐵𝐵 𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙 ( 12
)
where J is the omnidirectional flux calculated using the AP8min model. The formula given in [START_REF] Boscher | In-Flight Measurements of Radiation Environment on Board the Argentinean Satellite SAC-D[END_REF] allows the calculation of the unidirectional flux as a function of the magnetic field. However, in our case we compute the evolution of the unidirectional flux as a function of the equatorial pitch-angle. The perpendicular unidirectional flux corresponds to a local pitch-angle of 90° which allows us from VI. CONCLUSION Monte-Carlo simulations have been used investigate anisotropy seen by the Head-A of the ICARE-NG radiation monitor. By combining response function calculated by GEANT4 and real positions of the SAC-D satellite along its orbit we were able to calculate pitch-angles of particles which crossed the two diodes by simulations. Then, by comparing these pitch-angles with the one defining the loss-cone, we were able to say if the head-A sensor in coincidence mode has access to trapped particles or not.

Results in this paper highlight the anisotropy seen by the Head-A of the ICARE-NG instrument. Particles' equatorial pitch-angles impinging the sensor are not the same according to the satellite's location and according to ascending or descending orbits so it allows to characterize the anisotropy seen by ICARE-NG Head-A.

Finally, two models have been implemented to determine unidirectional fluxes out of in-flight measurements of ICARE-NG Head-A. The study of various models for the evolution of the unidirectional flux according to the equatorial pitch-angle shows that anisotropy depends on the energy of particles and that it is also present at low energy. This work could be used as a complement of specification models including the anisotropy.

All the methodology described in section IV and in section V is applicable for other missions for which detectors have strong directionality.

Fig. 1 .

 1 Fig. 1. Fluxes recorded by the Head-A during January 2015 for protons of 65MeV. Top: Descending orbit. Bottom: Ascending orbit.

Fig. 2 .

 2 Fig. 2. Geometric factor (cm 2 .sr) in an incident energy (X-axis) -deposited energy (Y-axis) map for ICARE-NG Head-A on-board SAC-D in coincidence mode.

Fig. 3 .

 3 Fig. 3. CARMEN and Head-A coordinate systems.

Fig. 4 .

 4 Fig. 4. Left: Geometry of the SAC-D satellite. Right: Orientation of the ICARE-NG axis in the satellite frame.

Fig. 5 .

 5 Fig. 5. Evolution of the maximum observable equatorial pitch-angle according to the longitude of SAC-D for a given drift shell with L=1.185 and an altitude equals to 664km.

Fig. 9 .Fig. 10 .

 910 Fig. 9. Result of the process for Valot's model and the Badhwar-Konradi model using ICARE-NG Head-A data for L=1.25 and E=18.6MeV with the ideal straight line in black (if calculated and measured counts were equal). Counts are expressed in number of particles detected per second.

TABLE 1 ICARE

 1

	-NG ENERGY DEPOSITION CHANNELS AND DETECTED PROTON
	ENERGY
	Channels	Energy (MeV)
	A coincidence 50-78	18.6
	A coincidence 9	115
	A coincidence 10	100
	A coincidence 11	90
	A coincidence 12	81
	A coincidence 13	74
	A coincidence 14	69
	A coincidence 15	65
	A coincidence 16	63

the local magnetic field 𝐵𝐵 𝑇𝑇𝑇𝑇𝑙𝑙𝑇𝑇𝑇𝑇 , to calculate the corresponding equatorial pitch-angle using [START_REF] Townsend | Implications of the space radiation environment for human exploration in deep space[END_REF].

By comparing results of the minimization process in Fig. 9 and Fig. 10 The model based on the Badhwar-Konradi method using AP8 for protons of 20MeV over-estimates fluxes in comparison to other models. The Badhwar-Konradi method based on ICARE-NG measurements is close to the Roederer model especially for pitch-angles higher than 66°.

Regarding Fig. 12, curves of models based on ICARE-NG Head-A measurements are very close for an energy of 81MeV except for pitch-angles close to the loss-cone pitch-angle. In general, the Valot's model based on ICARE-NG measurements goes down rapidly for angles near to the losscone pitch-angle since this model is simpler than the other ones. Curves for models based on AP8 have the same variation for 81MeV, although the shape parameter was calculated for 20MeV protons. The Badhwar-Konradi model based on ICARE-NG measurements seems to be the most appropriate in comparison with the Valot's model to characterize the evolution of the unidirectional flux with the equatorial pitch-angle for protons.

Some studies established that AP8 model overestimates protons fluxes. As indicated in Siegl et al. [START_REF] Siegl | Inner belt anisotropy investigations based on the standard radiation environment monitor (SREM)[END_REF], AP8 overestimates proton fluxes between L=1.3 and L=1.6. Heynderickx et al. [START_REF] Heynderickx | A low altitude trapped proton model for solar minimum conditions based on SAMPEX/PET data[END_REF] have shown that AP8 overestimates fluxes for low energy protons. These overestimations could have an impact on results obtained thanks to the AP8 model.


text:
P Oudayer

L Monnin

J.-C Matéo-Vélez

S Hess

P Sarrailh

G Murat

J.-F Roussel

Keywords: dust lofting, lunar dust, electrical charge

A key aspect of dust adhesion to space equipment is the accumulation of charge under the space plasma environment. Recent models and experiments show possible negative charging of dust grains under VUV illumination. Macroscopic potential measurements conducted during a test campaign show that both positive and negative average charging can be reached under VUV irradiation pending on vacuum chamber configuration, suggesting that both situations can exist at lunar surface. Simulations of dust charging at microscopic scale are conducted with the SPIS software to evaluate electrical charge and electric field amplifications induced by the granular structure of the lunar regolith. A multi-layer pile of dust is modelled under lunar conditions. Grains from the first two layers tend to microscopically acquire both negative and positive charge patches when illuminated with a 45° VUV incidence angle, this differential charging being less pronounced under normal incidence angle. It is also found that dust deeply buried in the lunar soil may charge more negative due to the collection of environmental electron only. This effect is thought to reinforce the grain supercharging model presented by other authors. We show however that such charge development may be limited by electrical conduction pending on dust electrical properties.

I. INTRODUCTION 1 One of the main fallouts from the Apollo program is the problem of lunar dust adhesion which was initially assumed to be of secondary importance. The dust layer covering the surface of the Moon tended to stick to any human-made equipment, and therefore complicated the full exploitation of the facilities. It also led to human body undesired effects such as eye irritation, respiratory problems, hay fever symptoms, etc [START_REF] Wagner | The Apollo Experience Lessons Learned for COnstellation Lunar Dust Management[END_REF]. The next space missions will be initiated by robots, with no human intervention and no longer confined to the Moon: planets, comets and asteroids are also targeted. It therefore becomes more and more crucial to anticipate and bring under control the overall dust adhesion problem. In order to achieve this goal, the first priority is to understand the physical phenomena at stake throughout experiments. Particle collection on equipment surfaces may be due to human and robotic activity, leading to dust mobilization from the object soil. Particle levitation from the Moon surface may also be due to the space environment interaction with its surface. Indeed, the Moon surface is likely to be negatively charged due to the solar-wind electron collection. On the day-side however, photoemission phenomenon occurs and leads the surface to be positively charged. At the frontier of day and night-side (called terminator), the presence of both negative and positive charges on small distances create an electric field believed to be strong enough to electrostatically loft dust particle at the Moon surface during the so-called horizon glow phenomenon ( [START_REF] Criswell | Horizon-Glow and the Motion of Lunar Dust[END_REF]- [START_REF] Horanyi | A permanent, asymmetric dust cloud around the Moon[END_REF]). This phenomenon is mainly seen at dawn and dusk. Both experimental and numerical investigations have recently been conducted to model the effect of the lunar plasma environment on dust charging and lofting ( [START_REF] Wang | Dust charging and transport on airless planetary bodies: Electrostatic Dust Transport[END_REF], [START_REF] Champlain | Lunar dust simulant charging and transport under UV irradiation in vacuum: Experiments and numerical modeling: ELECTROSTATIC TRANSPORT OF LUNAR DUST[END_REF], [START_REF] Necmi Cihan Örger | Investigation of Electrostatic Transportation of Lunar Dust Grains due to Ambient Plasma Conditions[END_REF]).

Independently of the origin of dust deposit onto surfaces, i.e. mobilization by robotic/rover/lander/human activity or natural lofting induced by the environment, adhesion of dust to exploration unit surfaces is mainly due to two forces: the Van der Waals force, which is a distance-dependent cohesive force, and the electrostatic forces [START_REF] Mittal | Particle Adhesion and Removal: Mittal/Particle[END_REF] which are a direct result of the lunar space environment since the charge depends on the complex interactions between the ambient plasma and dust. The intensity of these forces remains poorly quantified within the regolith.

In day-side and under VUV irradiation, the dust layer is expected to be positively charged [START_REF] Criswell | Intense localized photoelectric charging in the lunar sunset terminator region, 2. Supercharging at the progression of sunset[END_REF]. Experimental evidence of positively charged dust mobilization induced by VUV exposure was shown by [START_REF] Champlain | Lunar dust simulant charging and transport under UV irradiation in vacuum: Experiments and numerical modeling: ELECTROSTATIC TRANSPORT OF LUNAR DUST[END_REF]. Their setup consisted in applying a horizontal electric field using a positive electrode and a negative electrode both located a few millimeters from a dust layer. Dusts were collected by the negatively biased electrode and not on by a positively biased. In other experimental configurations however, with a grid located above the dust layer, dust mobilization occurred after VUV exposure only when a positive voltage was applied to the grid (Prof. M. Horanyi, private communication). This suggests that different experimental setups can produce different charging at large scale, i.e. at a few millimeters to centimeters scale, what is called macroscopic scale in this paper. One possible explanation for this is the presence of macroscopic electric fields that deflect the electrons emitted by VUV impacts on illuminated surfaces. In addition, complex microscopic charging patterns are supposed to be generated due to the interstices between dust grains that make charged particle

Multiscale Modelling of Dust Charging in Simulated

Lunar Environment Conditions P. Oudayer (1) , L. Monnin (1) , J.-C. Matéo-Vélez (1) , S. L. G. Hess (1) , P. Sarrailh (1) , G. Murat (1) , J.-F. Roussel (1) (1) ONERA -The French Aerospace Lab, Toulouse, France collection and emission very different from a smooth flat surface [START_REF] Wang | Dust charging and transport on airless planetary bodies: Electrostatic Dust Transport[END_REF]. Due to the very loose packing at the surface of the Moon, photoemission phenomena can occur on dust buried just beneath the top dust layer, i.e. what can be called the second dust layer. Some of the emitted photoelectrons can be collected at the bottom of particles of the first layer, locally creating a negative potential spot, while the upper side of the layer remains positively charged due to photoemission. This model was confirmed by numerical simulations [START_REF] Zimmerman | Grain-scale supercharging and breakdown on airless regoliths: REGOLITH SUPERCHARGING AND BREAKDOWN[END_REF]. Complementary experimental investigations showed that a typical difference potential of -2.5 V can be achieved between illuminated side and shaded side of two dust grains [START_REF] Wang | Dust charging and transport on airless planetary bodies: Electrostatic Dust Transport[END_REF].

To our knowledge no attempts have been done to validate the numerical modelling of dust charging with dedicated experiments both at the microscopic scale of the dust grains and at the macroscopic scale of the experimental setup. Both scales are indeed involved in the collection of charge carriers (electrons most of the time) by dust. The objective of this paper is to propose a combined study of dust layer charging under VUV exposure both experimentally and numerically.

Section II of this paper is dedicated to simulations of a dust pile under plasma environment using a microscopic grain arrangement similar to [START_REF] Zimmerman | Grain-scale supercharging and breakdown on airless regoliths: REGOLITH SUPERCHARGING AND BREAKDOWN[END_REF]. In Section III we replicate laboratory experiments with simulation at both macroscopic and microscopic scales. Main results are discussed and concluded in Section IV.

II. MICROSCOPIC MODELLING OF DUST LAYERS UNDER SPACE ENVIRONMENT

A. Modelling

Numerical simulations are conducted with the open-source SPIS software version 6.0 (www.spis.org). A superposition of three grains layers is simulated. The simulation domain is bounded by a rectangular domain of dimensions 0.8 mm, 0.6 mm and 1.8 mm along the X, Y and Z axis respectively. A potential of 0 V is applied on top of the simulation domain. The bottom boundary is set to a floating potential. Periodic conditions along the X and Y axis and reflective conditions along the X axis are implemented in order to simulate an infinite surface. Each dust grain is 100 µm in radius and separated by 1/10 th of the particle size with respect to its neighbors. The grains are not in contact with each other. Dust layers lay onto a flat surface, simulating the rest of it. The SPIS material used to characterize the dust refers to "Lunar dust 71501 Mare" and a bulk conductivity of 10 -14 S/m is imposed, as well as a surface resistivity of 10 15 Ω. No conduction path is modelled between each dust grain.

Two sets of simulations have been performed. The first one corresponds to a Sun incidence and to solar wind incidence angle α of 45°. The second one corresponds to a normal incidence. Solar wind plasma particles dynamics, i.e. electrons and protons, are simulated by Particle-In-Cell modeling. Maxwellian distributions are injected on the top boundary. Secondary electron emission by electron impacts (SEE) and by photoelectrons are simulated following PIC modelling and assuming Maxwellian distributions with a temperature of 2.2 eV. Plasma parameters are gathered in Table I. The plasma parameters used during a previous study with α equal to 45° are reported in Table II [START_REF] Zimmerman | Grain-scale supercharging and breakdown on airless regoliths: REGOLITH SUPERCHARGING AND BREAKDOWN[END_REF].

B. Results

Fig. 1 presents the surface potential reached at equilibrium after about 200 s for α = 45° and neglecting SEE.

Simulations show the apparition of both positive and negative patches on top and at the bottom of the first layer, respectively. As in [START_REF] Zimmerman | Grain-scale supercharging and breakdown on airless regoliths: REGOLITH SUPERCHARGING AND BREAKDOWN[END_REF], photoemission on top of the first layer generates positive charge while photoemission on the second Density (#/m 3 ) 10 8 5x10 6 5x10 6

Fig. 1. Surface potential at equilibrium with α = 45° and without SEE. Bottom figure is the cross section along the Y-axis.

layer generates negative charges on the rear side of the top layer. The simulated differential potential obtained in this work is about -7.5 V between the rear and front sides of the top layer. This is comparable to the -2.5 V voltage obtained experimentally on macroscopic objects mimicking dust under VUV illumination [START_REF] Wang | Dust charging and transport on airless planetary bodies: Electrostatic Dust Transport[END_REF]. However, the average potentials on the top layer and on the second layer are +1 and +2.5 V, respectively, suggesting super charging is compensated in average by photoemission on the first layer top surface. This equilibrium state is reached after about 10 seconds. The average potential reached after 200 seconds on the third layer is -13 V. This layer is not illuminated with directional VUV. It collects only solar wind electrons that find their way between grain interstices due to electrostatic deflection by the first and second layers. This collection needs a longer duration to reach equilibrium because electron fluxes are reduced by the presence of the first and second layers. The electric field obtained at 200 seconds is plotted in Fig. 2. Between, the top and second layer, we observe an electric field slightly greater than 2x10 5 V/m which is in good agreement with previous results [START_REF] Zimmerman | Grain-scale supercharging and breakdown on airless regoliths: REGOLITH SUPERCHARGING AND BREAKDOWN[END_REF]. The electric field between the second and third layer is a bit less than one order of magnitude greater than between the first and the second. That suggests that the maximum grain super charging occurs quite deeply in the soil.

Secondary electron emission impact is also studied. The main change when adding SEE is that the third layer is less negatively charged than previously, as seen on Fig. 3, i.e. -10 V instead of -13 V after 200 seconds. Fig. 3 presents the surface potentials for α = 45° and with SEE. No significant change is observed on the top and second layer because their potential is controlled by the photoemission process rather SEE. The third layer is a bit less negative (-10 V) than in the previous case because SEE tends to limit the negative current balance. SEE is important to model mainly for the deeper layers thus.

Results obtained with α = 0° and with SEE are presented in Fig. 4. It appears that the negative and positive patches are less pronounced. The potential differences between the top layer front and rear sides reach about -3.5 V. This is quite well in agreement with the -2.5 V measured in [START_REF] Wang | Dust charging and transport on airless planetary bodies: Electrostatic Dust Transport[END_REF]. For the second layer, the charging behavior is comparable to the previous case.

The main difference between this work and the previous study reported in [START_REF] Zimmerman | Grain-scale supercharging and breakdown on airless regoliths: REGOLITH SUPERCHARGING AND BREAKDOWN[END_REF] is the number of grain layers. In [START_REF] Zimmerman | Grain-scale supercharging and breakdown on airless regoliths: REGOLITH SUPERCHARGING AND BREAKDOWN[END_REF], one layer and a half is simulated instead of three in the present Fig. 2. Electric field at equilibrium in the same configuration as in Fig. 1. work. Same observations are made for the top layers: apparition of positive charges where photoemission occurs (on top of the first and second layer) and negative charges where photoelectrons are collected (at the bottom of the first layer). The same qualitative behaviors are observed for changes in the solar wind electron density and temperature within a factor of 2. The theory of grain super charging seems thus to be corroborated by the present results [START_REF] Wang | Dust charging and transport on airless planetary bodies: Electrostatic Dust Transport[END_REF]. In this work, we showed that grains buried beneath what we call the third layer can charge highly negative due to the infiltration of solar wind electron between the top layers. This may contribute to some change in the estimation of electrostatic forces.

III. MULTISCALE MODELLING IN SIMULATED LUNAR

ENVIRONMENT Microscopic scale modelling is thus possible to assess using numerical tools. Validating these tools with ground experiments is however difficult to achieve because of the difficulties to perform such experimental investigation at grain scale. In this section, we present experiments conducted to understand macroscopic charging and the efforts made to model these experiments with SPIS at this same macroscopic scale and at microscopic scale as well.

A. Experimental setup

Experiments have been conducted in the DROP vacuum chamber described in [START_REF] Champlain | Lunar dust simulant charging and transport under UV irradiation in vacuum: Experiments and numerical modeling: ELECTROSTATIC TRANSPORT OF LUNAR DUST[END_REF]. This facility is a cylindrical tank of dimension 40 centimeters in diameter and 40 centimeters in length covered with a hemispherical cap where a pressure of 10 -6 mbar is obtained with a turbo-molecular pump. The chamber is equipped with a VUV deuterium lamp emitting photons from 120 nm to 400 nm. The illuminated zone is about 10 centimeters in diameter, the lamp being placed 35 centimeters above the sample. Fig. 5 presents the setup used for surface potential measurements after VUV exposure. A total of three samples can be installed at the center of a holder plate made of copper. The samples are electrically connected to a bias voltage supply. A contactless Kelvin probe is mounted on a translating axis to measure the surface potential profile. The photoelectron flux emitted from a copper reference sample is about 40 nA/cm 2 [START_REF] Matéo-Vélez | Development of a predictive discharge numerical model on solar panels -Test Report[END_REF], which is about 20 times the expected photoemission flux PEY at 1 AU from the Sun. The UV lamp irradiance is thus supposed equivalent to 20 times the VUV flux at 1 AU from the Sun. Lunar dust simulant photoemission flux was an unknown parameter of the test setup and has to be determined.

Measurements of the potential are made using a Kelvin probe which is placed 2 mm above the copper plate. The sample thickness is about 200 µm, so the probe is placed 1.8 mm away from the dust layer. Samples used for the experiment are SMA_A lunar dust simulant, which grain size is inferior to 25 µm [START_REF] Freyria | Simulated Moon agglutinates obtained from zeolite precursor by means of a low-cost and scalable synthesis method[END_REF] . The average grain size is 8 µm. The layer is made by depositing dust into a small container of 1.5 cm of diameter using a height gauge to obtain a thickness of about 200 µm. After pumping the vacuum chamber down to 10 -6 mbar, the dust layer is irradiated for a duration of two minutes. After this period, the potential is measured above the centre point of the dust deposit. The first twenty seconds of data are lost by switching off the lamp and moving the probe. However, it is still possible to extrapolate the surface potential at the time the VUV is switched off using an exponential fit.

A first set of experiments have been performed with samples electrically connected to the chamber wall reference ground. A potential of -4 V has been measured just after VUV irradiation. Two hypotheses have been proposed to explain this negative charging. The first hypothesis follows the grain super charging model [START_REF] Wang | Dust charging and transport on airless planetary bodies: Electrostatic Dust Transport[END_REF] that assumes that negative charges on the top layer rear side overpasses the positive charging ion the top layer front side. Even though, the patch model has been corroborated by numerical simulation and by experiments at large scale [START_REF] Zimmerman | Grain-scale supercharging and breakdown on airless regoliths: REGOLITH SUPERCHARGING AND BREAKDOWN[END_REF], this hypothesis is difficult to assess experimentally at grain size. A patch structure does not necessarily means that the average potential would be negative. To clarify this point, a second hypothesis has been proposed. We indeed suspected a strong recollection by the dust layer of photoelectrons emitted by the copper plate. Low energy electrons emitted by the copper plate are indeed efficiently backscattered by the tank walls and reflected back inside the vacuum chamber, i.e. towards the dust layers [START_REF] Cimino | Can Low Energy Electrons Affect High Energy Physics Accelerators?[END_REF]. If that repelled flux locally exceeded the flux emitted by the dust, that would explain the observed negative charging. Under these conditions, the setup would act as if the dust layer was under VUV and very large flux of primary environmental electron. To check this hypothesis, a second experiment has been conducted by applying a bias potential of -45 V to the sample holder. The initial dust potential was thus -45 Volts. This technique is quite similar to the classical ones used for secondary electron emission measurement to reduce the backflow of repelled electrons where a positively biased hemispherical cup is placed above the sample [START_REF] Belhaj | The effects of incident electron current density and temperature on the total electron emission yield of polycrystalline CVD diamond[END_REF]. The potential measured after VUV irradiation is -33 V with respect to chamber ground. This corresponds to a +12 V differential potential with respect to the sample holder. That means that, in that configuration, dusts are positively charged in average. This experimental campaign thus suggests that positive and negative average charging can alternatively be reached experimentally under VUV illumination pending on the electric field at the vicinity of the dust layer. This may be the Fig. 5. Experimental setup inside the DROP chamber reason why both positive and negative dust lofting was observed in the literature under VUV. The next section aims at providing more details on the positive charging experimental configuration under VUV using numerical simulations at macroscopic and microscopic scale.

B. Numerical simulations

Because of the large difference between individual grains size and the chamber dimensions, the electron transport and dust surface potentials are computed in two steps. First, a large scale simulation is performed using the geometry presented in Fig. 6. The DROP chamber is modelled by a 20 cm-radius cupola. The circular target 1.5 cm in diameter surrounded by a rectangle of copper represents the dusts. The hole on the top represents the VUV lamp aperture. The mesh is refined at the dust layer location with mesh grid spacing of 0.5 mm. It progressively increases up to 10 mm on the external part of the copper plate and up to 50 mm on the rest of the vacuum chamber. Similarly to the second experimental campaign, a bias potential of -45 V has been imposed on the copper plate. All surface potentials are fixed except the dust one which is floating. In one hand, this simulation is used to compute the dust average surface potential. On the other hand, it is used to compute the electron fluxes coming in and out the dust layer vicinity to serve as inputs to a microscopic scale simulation.

The second simulation represents a dust pile geometry similar to what is shown in Section II. Grains are 25 µm-radius spheres and there are five layers, which results in a 0.2 mm thick sample. The mesh grid spacing is 2.5 µm on dust grain and it progressively increases up to 125 µm on the top boundary.

C. Results

1) Macroscopic simulation

With the default SPIS material properties for dust, the surface potential is slightly more positive (-30 V) than during the experiments. The default value of the dust photoemission current density at (PEY) is 0.45 nA/cm 2 . The flux of photoelectrons emitted by a surface perpendicular to the VUV source is thus 9 nA/cm2, taking account of the VUV source amplification of 20 with respect to expected lunar conditions. With a PEY of 0.35 nA/cm 2 instead, the average surface potential tends to the experimental value of -33 V. Fig . 7 presents a cross section of the plasma potential at 400 s, i.e. after reaching equilibrium. Fig. 8 shows the photoelectron charge number density. At equilibrium, most of the photoelectrons emitted by the target plate are recollected by the target itself due to the presence if the potential barrier. One could obtain larger positive voltages by using larger dust deposit diameters. In addition, the photoelectron current emitted by the copper plate and recollected by the dust target surface is I recoll = 3.3 nA. Fig. 9 presents the potential along the Z axis. A minimum potential of -34.4 V is obtained at a distance of 3.5 mm. This potential barrier is imposed by the -45 V copper plate and prevents the dust layer from charging more positive (or less negative) than -33 V with respect to the tank walls. Fig. 8. Photoelectron charge number density (m -3 ) at equilibrium.

2) Microscopic simulations

A quasi-1D domain is simulated, dimensioned using the location of the potential barrier created by the copper plate at approximately 3 mm from the dust layer. Above this barrier, all photoelectrons emitted by the dust layer are definitively lost, so that it can be modeled as an open boundary. A difference of potential of -10.6 V is applied between the open boundary and the copper plate corresponding to the difference between the copper plate potential and the minimum barrier potential. For numerical reasons, all potentials are shifted by 34. [START_REF] Wang | Dust charging and transport on airless planetary bodies: Electrostatic Dust Transport[END_REF] V so the open boundary potential is 0 V. This choice has not impact on the results since only potential gradients are important.

At this new scale, two electron populations are simulated: the photoelectrons emitted from the dust considering α = 0° and an isotropic Maxwellian population of electrons standing for the external photoelectron originating from the copper plate. Secondary electron emission is neglected here. The Maxwellian population is injected from the open boundary. Its density is determined in such a way that its thermal current density equals the current density of electrons coming from the copper plate obtained in the large scale simulation: ,

with I recoll being the previously calculated current (in Ampere), S is the target surface (in square meter), q is the electron charge (in Coulomb), k B is the Boltzmann constant (J.K -1 ), T e is the photoelectrons' temperature (in Kelvin) and m e their mass (in kilogram). The thermal electron temperature is fixed at 12 eV, which corresponds to the temperature of the electron emitted by the copper plus the potential difference between the copper plate and the potential barrier. Eventually, a density of N = 1.95 × 10 8 particles/m 3 is found. Fig. 10 presents the time evolution of the average potential over each layer. Fig. 11 shows a cross section of the surface potential after 450 seconds. The potential on the first two layers is controlled by the external electron collection and photoelectron emission and recollection. First, photoelectrons are emitted and the potential goes up. Then, a fraction of them is collected as well as external electrons and the potential goes down until it reaches equilibrium within a few seconds. The differential potential between the top layer 1 and the copper plate is +10 V, in agreement with the experimental results. Grains that are not enlightened at all (three last in depth layers) do not emit any photoelectron, but they can collect electrons from the Maxwellian population. Thus, the potential on their surfaces decreases by a few volts, until equilibrium. The third layer gets a -10 V potential with respect to the copper plate after about 10 minutes. The time scale to reach equilibrium is several hundreds of seconds. That means that experiments should be conducted over more than 10 minutes to get equilibrium state. The fourth and fifth layers collect very small amount of electrons.

The results show the same trends as in Section II suggesting that the experimental conditions produce results comparable to the expected lunar conditions. The chamber may be used to both simulate photoemission by dust and collection of ambient electrons. As for the simulation presented in Section II with α = 0°, the charge patch patterns are less pronounced. It would be interesting to test the effect of grazing incidence angle.

However, no charge conduction between layers is taken into account, which is not physically realistic.

3) Potential relaxation modelling

We report here a study of charge relaxation after VUV illumination both experimentally and numerically. The time evolution of the surface potential have been measured at the center of the dust sample after 2 minutes of VUV irradiation is presented. In order to model numerically the potential relaxation, a current leakage is introduced between the layers by imposing an arbitrary resistance between two consecutive layers of R = 10 17 Ω. The conducted current is computed by dividing the average potential difference between two layers by R. The value of 10 17 Ω is about ten times greater than the bulk internal resistance of each grain. The simulation duration is composed of 10 minutes initial charging under VUV followed by a relaxation of 110 minutes. Fig. 12 compares the potential evolution obtained both experimentally and numerically. Experimentally, surface potentials range from + 12 V to 0 V within 90 min. Numerically, during the charging period, i.e. during the first 10 minutes, a behavior similar to the previous 3-layers simulations is observed. The first two layers reach an average potential of approximately +8 and +5 V. The third layer potential is decreasing due to electron collection. It reaches a minimum value at -6 V, which is quite superior to the -25 V previously obtained for the same charging duration. The difference is due to the contact resistance simulated between layers which limit charge leakage from one layer to another. When the simulation reaches 10 min (VUV lamp is switch off), the two first layers show a potential decrease (the first decreasing more rapidly than the second one). Final simulated potentials reached are +4 and +2 V after 120 min. This difference potential is less important than the one reached from experimental results, indicating the simulated resistance contact value is too high. New estimation of the resistance contact is done using slope ratio between simulated and experimental potentials. In order to superpose both potential, it is found the simulated resistance contact must be about 10 16 Ω.

Third layer shows a potential increase, seemingly providing charge to the fourth and fifth layer, as these two last ones show a potential decrease. Complex charge redistribution occurs and would need to be taken into during both charge and relaxation phases. However, supposing a relaxation time long enough, all potentials are expected to tend towards the copper plate potential.

IV. CONCLUSIONS

This paper presented experimental and numerical evidence of complex charging patterns of dust submitted to lunar like environments. Average surface charging in sunlight may range from a few volts positive to a few volts negative pending on the soil geometrical configuration and electrical properties and on plasma conditions as well, including solar wind and sunlight incidence angle. We have shown that potential barriers are of prime importance to understand macroscopic charging behavior during ground testing. Extrapolating this observation to the lunar surface, one would expect for instance sunlit areas located at the border of negatively charged shaded surfaces to charge negatively due to the recollection of photoelectrons induced by potential barrier effects. This paper also corroborates literature models where dusts located on top of the lunar surface may exhibit positive and negative charge patches. We complement this model by showing the possibility of some dust buried quite deeply beneath the surface can get highly negative due to ambient electron collection. The charging time scale for these buried dusts is about a few hundreds of seconds while charging of the top layers dust is about a few seconds. That makes them more sensitive to electrical conduction mechanisms. The equilibrium potential of buried dust depends a lot on the soil bulk and surface conductivity as well as on the contact resistance between grains. The potential of dust located at the top surface is less sensitive but probably not totally insensitive to conduction processes, especially for smaller grains.

Electrostatic dust lofting has been reported in [START_REF] Wang | Dust charging and transport on airless planetary bodies: Electrostatic Dust Transport[END_REF] on laboratory experiments even though no electrostatic force estimation could theoretically overcome adhesion and gravity. The observations made by the authors indicated aggregate lofting rather than single grain. The present paper, by showing that dust buried well beneath the top surface can get highly negatively charged, tends to corroborate what is suspected in [START_REF] Wang | Dust charging and transport on airless planetary bodies: Electrostatic Dust Transport[END_REF]: electrical charge carriers can reach deeply inside aggregates and provide them sufficient charge level to get them mobilized by electrostatic lofting by an external electric field, that depends on large scale geometry, and/or by electrostatic repulsion with other dust and/or agglomerates, that is more related to small scale interaction.

Future work should include more conduction current measurements in order to assess charging attenuation by current leakage through dust and from grain to grain, pending on the lunar soil to be simulated. This is indeed a key parameter to understand charging patterns and their effect on dust mobilization and adhesion both at macroscopic and microscopic scales.

Fig. 3 .Fig. 4 .

 34 Fig.3. Surface potential at equilibrium, with α = 45° and with SEE

Fig. 6 .Fig. 7 .

 67 Fig. 6. 3D modelling of the DROP chamber simulated at macroscopic scale

Fig. 9 .Fig. 10 .Fig. 11 .

 91011 Fig. 9. Potential profile along the Z axis above the dust target.

Fig. 113

 113 Fig. 113 presents a cross section of the five layers simulation at 120 min. No more charge patches are observable, as resistance and conductivity have homogenized the potentials along the grains.

Fig. 12 .

 12 Fig. 12. Evolution of numerical surface potentials of layers 1 to 5 versus time (in continuous lines) and envelope of experimental potential relaxation measurements (in grey).

Fig. 13 .

 13 Fig. 13. Cross section along the Y axis showing surface potential at 120 min during potential relaxation.

V. ACKNOWLEDGEMENTS

Authors would like to thank Serena Esposito from the Università degli Studi di Cassino and Francesca S. Freyria from the Instituto Politecnico di Torino for providing the SMA_A lunar dust simulant used in this paper.

VI.