Florian Thuillet
email: florian.thuillet@oca.eu

Patrick Michel

Shogo Tachibana

Ronald-L Ballouz

Stephen R Schwartz

Numerical modeling of medium-speed impacts on a granular surface in a low-gravity environment Application to Hayabusa2 sampling mechanism

Keywords: Minor planets, asteroids, individual, (162173) Ryugu -Methods, numerical

Even if craters are very common on Solar System body surfaces, crater formation in granular media such as the ones covering most of visited asteroids still needs to be better understood, above all in low-gravity environments. JAXA's sample return mission Hayabusa2, currently visiting asteroid (162173) Ryugu, is a perfect opportunity for studying medium-speed impacts into granular matter, since its sampling mechanism partly consists of a 300 m s -1 impact. In this paper, we look at medium-speed impacts, from 50 to 300 m s -1 , into a granular material bed, to better understand crater formation and ejecta characteristics. We then consider the sampler horn of Hayabusa2 sampling mechanism and monitor the distribution of particles inside the horn. We find that the cratering process is much longer under low gravity, and that the crater formation mechanism does not seem to depend on the impact speed, in the considered range. The Z-model seems to rightly represent our velocity field for a steady excavation state. From the impact, less than 10% is transmitted into the target, and grains are ejected mostly with angles between 48 • and 54 • . Concerning the sampling mechanism, we find that for most of the simulations, the science goal of 100 mg is fulfilled, and that a second impact increases the number of ejecta but not necessarily the number of collected particles.

INTRODUCTION

Since the discovery of a fine-grained layer on the surface of asteroids (Robinson et al. 2002), termed regolith, understanding the dynamics of granular material in low-gravity environments has become crucial. This knowledge could help us better comprehend the outcomes of both low-speed and high-speed impacts on the surface of small bodies [START_REF] Katsuragi | Physics of Soft Impact and Cratering[END_REF], and their surface evolutions and histories [START_REF] Asphaug | [END_REF][START_REF] Melosh | Planetary Surface Processes[END_REF]. Moreover, the development and launch of two asteroid sample return missions, JAXA's Hayabusa2 (Watanabe et al. 2017) and NASA's OSIRIS-REx (Lauretta et al. 2017), further motivated research in this field. Indeed, understanding granular material serves a scientific purpose, but also an engineering one, particularly in the case of sample collection, where understanding and predicting the interactions with the surface is fundamental.

One of the main sources of information on Solar System bodies surface properties and histories are impact craters. Craters are the most frequently and easily observed surface features by space probes that performed a fly-by or a rendezvous to a small body. A planetary body's crater morphology and dimensions can help infer the physical properties of its surface and sub-surface. Furthermore, the size distribution of craters is a hint for the age of a surface, but craters have to be discriminated into categories to determine if they have been created by an exogenous impactor (primary craters) or by ejecta fallout (secondary craters). The correct interpretation of craters on the surface of asteroids requires a good understanding of crater formation on low-gravity surfaces made of regolith, for both low-and high-speed impacts [START_REF] Hirata | [END_REF]Walsh et al. 2019).

Low-speed impacts on asteroids have also had a human origin, and are one of the chosen solutions for sampler return missions to collect material. JAXA chose for both Hayabusa and Hayabusa2 missions a sampling mechanism consisting of a projectile impacting the asteroid surface at 300 m s -1 (Sawada et al. 2017). Thus, increasing our knowledge on lowspeed impacts on regolith material is also beneficial for the design of sampling mechanisms and for the interpretation of their outcome.

Analytical formulae have been derived from experiments to scale the cratering process in a general way, most often focusing on impacts that are equivalent to an explosion and for which a point source solution applies [START_REF] Melosh | Impact Cratering : A Geologic Process[END_REF]Holsapple 1993). Some were more interested in very low speed impacts (less than 4 m s -1 ), and among them different parts were emphasized, but always in the context of Earth's gravity, such as the restraining force (Katsuragi & Durian 2007;Katsuragi & Blum 2017) or the crater's shape and size (Uehara et al. 2003;Walsh et al. 2003;de Vet & de Bruyn 2007;[START_REF] Nordstrom | APS March Meeting Abstracts[END_REF]. Others conducted experiments with much faster projectiles, closer to the speed used for the Hayabusa2 sampling mechanism: for example [START_REF] Yamamoto | [END_REF] looked at the velocity distribution of ejecta, Yamamoto et al. (2009) at the transient crater growth, and Nakamura et al. (2013) at the penetration depth of the impactor. Experiments have even been possible under low-gravity or microgravity, through the use of an Atwood machine [START_REF] Murdoch | [END_REF] or parabolic flights (Nakamura et al. 2013;Brisset et al. 2018).

However, these experiments require a lot of time to set up, are limited by available material and equipment, and rarely achieve true micro-gravity conditions. In parallel to these experiments, numerical simulations could be developed and then compared to actual data. They enable an exploration of a wide parameter space and make it far easier to measure a system's physical properties that are otherwise difficult to collect. Comparisons between experiments and numerical simulations could then been done, for instance penetration depth measurements (Nakamura et al. 2013) and crater formation (Wada et al. 2006).

Before it was eventually decided that the Hayabusa2 sampling mechanism would adopt the same design as the Hayabusa mechanism with a semi-spherical projectile, experiments were done to study the influence of the projectile shape [START_REF] Makabe | 26th International Conference on Space Technology and Science (ISTS)[END_REF] on the amount of ejected mass. Numerical comparisons by Schwartz et al. (2014) were then done, but limited to experiments with impact speeds of 11 m s -1 in Earth gravity. Our study builds on an improved version of the code used by Schwartz et al. (2014) and considers the low-gravity environment of the Hayabusa2 target (162173) Ryugu (Binzel et al. 2001;Wada et al. 2018). Ryugu's surface is represented as a collection of grains with assumed physical and mechanical properties. We then develop numerical simulations of the Hayabusa2 sampling projectile impact with the surface to provide estimates of the amount of ejecta produced by the impact as well as the excavation and crater formation, as a function of assumed surface properties. We also compare our results to the numerical simulations of Wada et al. (2006) and the experiments of Housen et al. (1983); [START_REF] Yamamoto | Proceedings of the 36th ISAS Lunar and Planetary Symposium[END_REF][START_REF] Yamamoto | [END_REF]Yamamoto et al. ( , 2006)). Finally, we perform numerical simulations that include the exact geometry of Hayabusa2 sampling horn, in addition to the projectile, to monitor the amount of ejecta that may be captured in the different parts of the horn.

In Section 2, we describe the characteristics of our grains and briefly our method, and in Section 3 we present our results, without and with the sampling horn. In Section 4, we provide a discussion of these results and an outlook on future work.

METHOD

In this section, we describe the different choices we made to model Ryugu's surface, and the quantities we investigated and monitored.

Simulation parameters and setup

Our simulations were performed with the N -body gravity tree code pkdgrav (Richardson et al. 2000;Stadel 2001;Richardson et al. 2009Richardson et al. , 2011)). In order to model the interaction between regolith grains, we use the Soft-Sphere Discrete Element Method (SSDEM) version developed by Schwartz et al. (2012), improved with a new rotational resistance model for the grains (Zhang et al. 2017), and the addition of "reactive walls", i.e., inertial walls that react to particles' forces, contrarily to regular walls, [START_REF] Maurel | [END_REF]. The version used here is the same as the one used in [START_REF] Thuillet | [END_REF], and several comparisons with experiments have been run all along the development of the code to check the validity of the results, from silo discharges (Schwartz et al. 2012), to projectile penetrations (Schwartz et al. 2014) and measurements of the angle of repose for different materials (Yu et al. 2014;[START_REF] Maurel | [END_REF].

In order to represent the environment in which the Hayabusa2's sampling will be performed, we considered the low-gravity environment of the asteroid Ryugu, assuming a constant gravitational acceleration. We then investigate the sampling process only in a close area around the sampling system. The gravitational acceleration has been computed from the measurements done by Hayabusa2 (Watanabe et al. 2019) and from the approximate location of the sampling, leading to a value of g = 1.19 • 10 -4 m s -2 .

We first ran simulations of the impact of the projectile without the sampling horn to study the characteristics of the impact, considering a regolith bed contained in a 15 cm radius and 20 cm deep cylinder. The projectile is modeled by a 5 g sphere with a radius of 0.4 cm. The real projectile is not exactly a sphere but the part impacting the regolith bed is spherical, and we assume that the response of the regolith is similar for an entirely spherical projectile. Schwartz et al. (2014) also used a spherical projectile in their comparisons between experiments and simulations with pkdgrav, and found good matches. We use a cylinder radius as a compromise between the computation time (a larger radius implying a much larger number of particles) and the long duration of the cratering process in low gravity. The bed is also defined to be larger than the bottom of the sampling horn (13.8 cm diameter).

Our bed is made of 101, 657 particles, with a particle density of 2.43 g cm -3 as in [START_REF] Thuillet | [END_REF]. The bulk density of the whole medium is about 1.23 g cm -3 , implying a macro-porosity of ∼ 50%. The size distribution of our particles is assumed to be Gaussian, with a mean radius of 0.25 cm, a standard deviation σ of 30% and a cut-off at 1σ. The choice of particle size was motivated by the opening of Hayabusa2 sampler's filter, which prevents particles larger than 1 cm in diameter to be ingested. The compromise was to consider particles small enough to be able to go through the horn and filters (with considerations on the computation time), but not too small to determine whether the filter will be clogged by larger particles.

The bed was created by randomly generating particles inside and above the cylinder, and allowing them to free-fall under Earth gravity. Once the bed has begun to settle, we switch the gravity to that on Ryugu, and we let the bed relax until a very low average speed (RMS speed), less than 3 • 10 -4 m s -1 , is achieved. Particle interactions are defined in pkdgrav by several parameters that include various friction and energy dissipation coefficients (Schwartz et al. 2012;[START_REF] Thuillet | [END_REF]. We took the same parameters as for the gravel-like material in [START_REF] Thuillet | [END_REF]. Concerning friction coefficients, the chosen values are µ s = 1.0 for the static friction, µ r = 1.0 for the rolling friction, µ t = 1.0 for the twisting friction, and β = 1.0 for the shape factor, a parameter representing the fact that real particles are not perfectly spherical (Zhang et al. 2017). Concerning the energy dissipation, the normal and tangential coefficients of restitution ε n and ε t are both set to 0.5. Nevertheless, we study the influence of these coefficients in Section 3.3.1. More explanations concerning these coefficients can be found in [START_REF] Thuillet | [END_REF].

In order to increase the range of applications of our modeling work, we consider impact speeds ranging from 50 m s -1 to 300 m s -1 . Cross-section snapshots for 50 m s -1 and 300 m s -1 are shown as examples in Sect. 3.1.

The timestep was defined so that the contact between particles and possible overlaps are well resolved for the considered dynamics of the system. For instance for the 300 m s -1 impact simulations, we chose a timestep ∆ t ≈ 2.5 • 10 -7 s. We look at the influence of these parameters in Section 3.1.

Analysis method

Our objective is to model Hayabusa2's sampling and check its efficiency for the considered regolith properties. Following Wada et al. (2006), we also examine several impact characteristics such as the ejecta volume, the ejecta speed and the evolution of the crater as a function of time. We also look at the energy distribution and the wave propagation during the first impact times, and look for potential links between the impact process and the outcome. Moreover, we compare our results to the Z-model theory [START_REF] Maxwell | Impact and Explosion Cratering: Planetary and Terrestrial Implications[END_REF].

Concerning simulations that include the sampler itself (the horn surrounding the projectile), we measure the quantity of actual sampled material for cases where particles are collected. Also, we measure the volume of material going through each part of the sampler, to determine the effect of the sampler's geometry and of the filter.

RESULTS

In this section, we describe the different results of our simulations. First, we present cases without the sampling horn, covering a range of impact speeds, in order to provide more general results on the impact process and its outcome into a granular material in reduced gravity (here chosen as Ryugu's one). These cases allow us to study more general outcomes, such as the crater formation, the ejecta production and direction, and to compare them with those given by a theoretical model (Sect. 3.1). Then, we present cases where we include the sample horn in order to measure the amount of collected material (Sect. 3.2). Finally, we study the influence of material properties and impact geometry (Sec. 3.3.2).

Impact outcomes for several impact speeds

General results

The impact speeds considered here range from 50 m s -1 to 300 m s -1 , with a particular focus on 300 m s -1 impacts.

Cross-section snapshots of simulations at different times for the 300 m s -1 impact are shown on Fig. 1. The crosssections are 1 cm wide, and represent the slice x = ±0.5cm of the regolith bed, i.e., the directions we see are the yand z-axes. The snapshots are representative of the threedimensional cavity.

When the projectile hits the bed, the wave propagates through the medium and particles are ejected very quickly. Several particles located very close to the impact point, on the top of the surface, and which are not directly hit by the projectile, linger in free fall since particles below them have been pushed downwards or not in their directions. They have a very small speed due to friction with particles they were in contact with, but these speeds are negligible compared to the speeds of the other particles around them. Because the gravity field is weak, they fall very slowly, and we can see one of them in the last snapshot, after 1230 ms (see Fig. 1).

Another way to visualize the cratering process is to consider the whole cylinder instead of a slice, and to take advantage of its axial symmetry. Figure 2 shows the particle density of the regolith bed, for different radii and heights, and for the same times as in Fig. 1. For each cell representing a toroid, the number of particles is summed along the azimuth and then divided by the toroid's volume.

The crater formation is visible in Fig. 2, as well as the particles ejected from the bed. This is also a way to ensure that the crater dimensions observed in cross sections are the same for all directions, and not only along a preferred one.

We notice that the crater formation in our simulations is much longer than in previous studies such as Wada et al. (2006). It is mostly due to the very low value of the gravity acceleration, compared to the terrestrial one used in Wada et al. (2006). Also, we see that radial boundaries certainly play a role in the cratering process, as the radius of the crater is very close to the radius of the cylinder. However, similar processes can be observed in our setup and the one of Wada et al. (2006).

After the impact, particles are pushed away radially, at the same speed, leading to a hemispherical cavity, as can be seen in the top left panel of Fig. 1. This was observed in laboratory experiments [START_REF] Yamamoto | Proceedings of the 36th ISAS Lunar and Planetary Symposium[END_REF]) and simulations (O'Keefe & Ahrens 1999; Wada et al. 2006). However, about 30 ms after the impact, the cavity depth does not increase as fast as the cavity radius anymore. After about 300 ms, the cavity depth has almost reached its final value, meanwhile the cavity's radius keeps growing, as shown in Fig. 1. This was also observed in Earth-gravity laboratory experiments [START_REF] Melosh | Impact Cratering : A Geologic Process[END_REF][START_REF] Yamamoto | Proceedings of the 36th ISAS Lunar and Planetary Symposium[END_REF], as well as in numerical simulations (Wada et al. 2006), and we can see that the phenomenon happens even under very weak gravity.

The cavity radius and depth as function of time are shown in Fig. 3. To calculate the cavity radius and depth, we counted the number of particles in different 1 mm layers, taking into account only the center of the cylinder and horizontal layers for the depth, and the surface and vertical layers for the radius. We then set a threshold of particle density to define the frontier between the cavity and regolith bed and determine the crater dimensions as a function of time.

During the first stage, particles are pushed radially away from the impact, both depth and radius increase at the same speed, and the depth as a function of the radius is close to the y = x line. This stage is however very short and in a second stage, the cavity depth then increases more slowly. The cavity radius increase also slows down, as particles constituting the crater rim have lower speeds than the ones ejected from the center. The cavity depth stops increasing after about 1 s, whereas the lateral growth continues.

We can also compute the cavity depth as a function of the cavity radius, as shown on Fig. 4, for four different speeds. The higher the impact speed, the faster the increase in cavity dimensions. However, we observe that the evolution of the crater depth as a function of the crater radius is very similar regardless of the impact speed. Yamamoto et al. (2006) and Wada et al. (2006) define the transient crater as the time when, the bowl-shaped cavity being formed and the cavity depth being stabilized, the crater rim starts collapsing and an uplift of the crater's bottom occurs. Our simulations do not last long enough to observe such an uplift and a collapse of the crater's walls. Moreover, the cylinder is too small to avoid boundary effects, and enlarging it and extending the duration would significantly increase the computation time. At the end of our simulations, which cannot be run until the transient crater is completely formed, we find a crater depth-to-radius ratio of about 0.6-0.65, which keeps slowly decreasing with time as the crater goes on expanding laterally, to be compared of the value of 0.52-0.54 observed in Yamamoto et al. (2006). The main difference between these experiments and our simulations, besides the fact that the transient crater is not reached yet in our simulations, is the gravity that is much weaker in our simulations. Therefore, particles can travel further away from the impact location, and therefore a crater should form with a lower depth-to-radius ratio. However, this effect is compensated by our boundary conditions. Moreover, the porosity, the internal friction, the grain size and the densities of both the projectile and the target also differ from the experiments, and Yamamoto et al. (2006) showed that these characteristics do not seem to affect very much the transient crater depth-to-radius ratio, even if they certainly affect the transient crater growth.

Our values (0.6-0.65) may seem very high compared to the depth-to-radius ratio of natural craters. For example, the largest craters observed on Ryugu had ratios between 0.28 and 0.4 (Sugita et al. 2019), and on Bennu the ratios were about 0.32 [START_REF] Barnouin | [END_REF]). However, these values correspond to final craters ratio and, as stated previously, our ratios correspond to transient craters. The transient crater ratio decreases with time, when the rim collapses and sometimes the floor lifts up. Yamamoto et al. (2006) found depthto-radius ratios of about 0.52-0.54 for transient craters and 0.22-0.28 for final ones, showing that our simulations would certainly lead to much lower ratios than 0.6 if we could reach the final state of the crater formation. Moreover, for actual asteroids, meteorite impacts can induce seismicity, leading to crater relaxation and lower observed depth-to-radius ratios (Richardson et al. 2004).

We did not only explore the crater dimensions, but also the impact energetics by studying how the energy is distributed during the first instance after impact.

Energy distribution during the first impact instants

We computed the different energies present in the first instants of the simulations, for the different considered impact speeds between 50 m s -1 to 300 m s -1 , with a particular focus on the 300 m s -1 case.

In our simulations, there are four types of energy: the translational kinetic energy K, the rotational energy E rot , the stored elastic energy E el , and the gravitational potential energy. We studied these energies for both the impactor and the target (sum of all the particles forming the granular bed), as well as the total energy in the simulation. We found that the gravitational potential energy was always negligible and therefore is never included in the following figures.

The kinetic energies are computed as following:

K = 1 2 mv 2 , (1)
E rot = 1 2 Iω 2 (2)
where m is the mass, v the translational speed, I the moment of inertia, and ω the angular velocity.

The elastic energies stored in each particle/particle or particle/wall contact come from the repulsive restoring spring force defined by Hooke's law, which is the interaction model used in pkdgrav (Schwartz et al. 2012). The stored normal, tangential, rolling, and twisting elastic energies can be defined as :

E el, norm = 1 2 k n d 2 overlap , (3)
E el, tang = 1 2 k t x 2 tang , (4)
E el, roll = 1 2 k n β 2 R 2 red θ 2 roll , (5)
E el, twist = k t β 2 R 2 red θ 2 twist (6)
where k n is the normal spring constant (see Schwartz et al. (2012)), d overlap the overlap length between the two considered particles, k t = 2 7 k n the tangential spring constant, β the shape parameter, and R red the reduced radius. x tang , θ roll , and θ twist represent the distention/compression of the respective springs, the rolling and twisting springs being angular springs. Then, equations from 1 to 6 should be summed over all considered particles.

We considered that the energies were equally distributed between particles in contact, and therefore Eq. 3, which represents the stored normal elastic energy for a pair of particles in contact, has to be halved after being summed on all particles to consider contacts only once.

We normalized the different energies by dividing them by the impactor's initial translational kinetic energy, the largely predominant energy prior to the impact, in order to facilitate comparisons between various impact speeds involving different orders of magnitude of energies. The different energies as function of time are shown on Fig. 5. The timestep and spring constant, respectively ∆ t and k n , have been scaled depending on the impact speed, i.e., the timestep used for the 300 m s -1 impact is six times shorter than the one used for the 50 m s -1 , and the spring constant is 6 2 = 36 times higher. Indeed, we considered a timestep varying proportionally to the inverse of the impact speed, and according to Schwartz et al. (2012), the timestep is proportional to the inverse of the square root of the spring constant. These choices of timestep and spring constant allow us to keep the same constant value for ∆ t • √ k n , and to have access to longer durations for smaller speeds.

Impacts evolve similarly regardless of the impact speed. The major differences between cases with different impact speeds are the timescale and the absolute amounts of energy. As expected, the higher the impact speed, the faster the energy dissipation and transmission to the target, and the higher the energy at the end of the collision (when the projectile kinetic energy is negligible compared to other energies). However, the normalized energy profiles are very similar, with the same peaks and a little less than 10% energy left at the end (and still decreasing with time). The fact that the total amount of kinetic energy in the ejecta is less than 10% of the initial kinetic energy was also found by Colwell (2003) in microgravity impact experiments and in high-speed experiments [START_REF] Davis | Asteroids III[END_REF].

Looking more closely at the evolution of the energy a a function of time, we can discern several steps. At first, the impactor's kinetic energy is by far the highest energy source (that is also an indicator that our bed is completely settled before the impact, as the sum of the kinetic energies of all particles forming the target is close to zero and negligible compared to the impactor's energy). Quickly, the projectile's kinetic energy is transformed into elastic energy as particle overlaps increase (for both the impactor and the target), and, simultaneously, part of this elastic energy is transformed back into target's kinetic energy. The impact wave propagates through the particles closest to the impact point, and overlaps between target's particles not in contact with the projectile increase, leading to a higher elastic energy for the target. The impact can also make the impactor spin, as it never impacts perfectly on the center of a particle. This spin is represented in Fig. 5 by the rotational energy. The target particles' rotational energies stay very low during the impact; even if particles have a non-zero rotation rate, their mass is much lower compared to the impactor's.

When the impactor penetrates into the bed, overlaps with target particles steadily increase, and this stored elastic energy is eventually transformed into kinetic energy. Therefore, the impactor is slowed down while particles are accelerated. Overlaps between the impactor and particles will therefore decrease, leading to a fall of target's elastic energy. During this short moment, particles almost do not stop the motion of the impactor anymore, which is represented by a little "bump" in the impactor's kinetic energy on Fig. 5.

After having reached its peak, the target's kinetic energy decreases with time because of dissipation through the impact wave. Energy dissipation in the medium comes from friction and damping (the coefficients of restitution are smaller than 1). After 10 -4 s for a 300 m s -1 impact (see Fig. 5), most of the energy left in the simulation is contained in the target (the target's kinetic energy is predominant), and the projectile has been almost completely stopped. The energy will then slowly decrease due to friction and energy dissipation in grain interactions in the regolith bed.

These different phases are shown on Fig. 6 for the 300 m s -1 impact, only for particles close to the impact point. The representations are cross-sections along the vertical plane containing both the initial position of the projectile and the center of the first particle impacted, for different times. They aim at a better understanding of the energy distribution, by looking at the first contacts with the impactor and the bed. In Fig. 6a, we show the first projectile-target contact, corresponding to a time of about 0.018 ms in Fig. 5, when the impactor's kinetic energy starts decreasing. The kinetic energy decreases quickly, due to the increasing impactorparticle overlap and overlaps between particles, as shown in Fig. 6b. After about 0.33 ms, the first particle to be impacted has pushed away other particles, among them the particle under the second particle to be in contact with the impactor, which enables particles to find less resistance in their path (and have speeds comparable to the impactor), and therefore the impactor can push them more easily and lose less kinetic energy (Fig. 6c). However, the impactor comes quickly in contact with a new particle (about 9 • 10 -6 s later), and therefore it is slowed down again, as shown in Fig. 6d.

We can also analyze the number of particles the projectile collides with during its penetration. Fig. 7 shows the number of particles that are in contact at least once with the projectile during this short simulation.

In Fig. 7, we see a first overlap between a particle (n • 8860) and the projectile, corresponding to the steep decrease of the target's kinetic energy. Then, the overlap decreases because n • 8860 has gained speed, and this corresponds to the "bump" in Fig. 5. The second drop in the projectile's kinetic energy happens when projectile forms new contacts with target particles (in this case it forms four simultaneously). The projectile cannot transfer as much en- . Zoomed-in cross-section representations of the projectile impacting the regolith bed, in the vertical plane containing both the initial position of the impactor and the first particle in contact, and for different times. The impactor is in green, and in red are represented the velocities of the different particles and impactor in this plane. The impact speed is 300 m s -1 , and the cross-section width is 1 cm. Distances are in centimeters. ergy as before, and therefore the contacts last longer and the projectile's kinetic energy drops although overlaps are small.

Since we ran simulations with lower impact speeds with scaled timestep and spring constant, one can wonder the influence of these parameters, and the changes in the impact process if these are kept the same for every speed. Here we consider a 50 m s -1 impact, with ∆ t 0 = 1.5 • 10 -6 s and k n 0 = 5.47 • 10 9 kg s -2 , i.e., the timestep and spring constant we used for the previous simulations. The considered range of timestep and spring constant is relatively small, and do not change the packing before impact. For example, for the simulations presented in Fig. 8a and8b, there are respectively 65 and 390 iterations before the impact, and the maximum displacement per iteration is 6 • 10 -6 cm, leading to less than 2.5 • 10 -3 cm in total from the beginning of the simulation to the impact.

If we decrease the spring constant from k n 0 to 0.58 k n 0 , the impactor keeps being slowed down by particles, as shown on Fig. 8a, contrarily to what was observed in Fig. 5. Thus, its kinetic energy decreases much faster. Since the spring constant is smaller, the spring repulsive forces are slightly weaker and particles experience longer contacts, leading to higher elastic energies for both the impactor and the target. We particularly notice that the tangential and rolling elastic energies are noticeably higher.

When particles initially pushed by the impactor encounter other particles, they are slowed down, and since in this case the impactor stays in contact, the slowing down increases the impactor-particles overlap and the former feels a higher, upward-directed, spring repulsive force than with a higher spring constant. This force is strong enough to change its trajectory, and the impactor bounces on the surface. However, the energy contained in the ejecta and the total energy are similar to the ones obtained with k n 0 . Since we are more interested in the amount of ejecta than in the fate of the impactor, a lower spring constant can be considered for further simulations.

Figure 9 shows the overlaps between the impactor and the five particles which are in contact with it at one moment of the simulation. It confirms that overlaps are larger and last longer in this case, i.e., when the spring constant is smaller. The difference is even clearer for the contact with the second particle (particle n • 37046).

If we look more closely at the moment that was just before the phase with less resistance and compare what happens for the two different spring constants, we notice that the contact between the second particle and the impactor is slightly different in both cases (see Fig. 10).

Indeed, the velocity of the second particle (the one below the impactor) does not have the same orientation, due to the deeper penetration of the impactor into the first particle. Since granular media are usually anisotropic, the bed does not answer in the same way if pushed in one direction or another, and this answer depends on the structure or preparation of the medium. For each contact, the target particle does not have the same mobility freedom in every direction, and this has to be taken into account for each particle along the created force chain. In Fig. 10b, the particle encounters more resistance along the direction of its velocity, and cannot move as well as in Fig. 10a, leading to the impactor bouncing away from the bed.

We explained how we scaled the timestep and the spring constant depending on the impact speed, and found very similar profiles. However, we also tried to compare different speeds with the same parameters. For this purpose, for simulations with impact speeds from 50 to 300 m s -1 , we harmonized timesteps and spring constants to respectively the lowest and highest ones, i.e., the ones of the 300 m s -1 impact. For example, for the 50 m s -1 , this leads to ∆ t = ∆ t 0 6 and k n = 36 k n 0 . The energy distribution for these parameters is shown on Fig. 8b. In this case, with a higher spring constant, the first impactor-particles contact pushes particles away quicker, and elastic energy increases very fast. Also, the repulsive force from particles is higher and the impactor's kinetic energy falls as fast as the elastic energy increases.

However, since the impactor and particles are pushed away from one another efficiently, the impactor-particles overlaps also shorten soon. The same phenomenon happens as it does with a spring constant of ∆ t 0 , where the impactorparticles spring forces become weaker as they move away, except this time the higher spring constant implies a longer time for the impactor to encounter low to no resistance from the bed. Where it does again, we can see once again the two phases "quick slow down" and "free penetration", softened because of a lower speed and a medium becoming more dense. At the end, the kinetic energy of the target and the total energy are almost exactly the same as for other simulations. We therefore expect a similar behavior for the target in a larger point of view, for example the cavity dimensions as we have seen in Section 3.1.1, or the amount of ejecta (see Section 3.1.5).

Wave propagation at the first impact instants

We stated in 3.1.1 that, just after the impact, the impact wave propagates through the medium. This is shown for the 300 m s -1 impact on Fig. 11, where we represent the speed of particles as a function of their distance from the impact point. Only particles located under the impact point, within a 60 • half-apex-angle downward cone, are considered in this figure, to avoid taking into account ejecta and to capture only the propagation of the wave inside the granular medium.

In Fig. 11, we see that particles far from the impact point have smaller speeds and are more numerous than close particles. This is due to the wave propagating through the medium, and being attenuated by it. In every contact, energy is lost due to friction and dissipation, leading to the decrease of particle speed with distance. Also, we can see the propagation of the wave in Fig. 11b, since particles further from 0.1 m have been significantly accelerated, and fastest particles, close to the impact point, have shared their energy with their neighbours and slowed down.

Shock propagation theories [START_REF] Melosh | Impact Cratering : A Geologic Process[END_REF]) state that particle speed decreases with distance following a power law, and that the power-law exponent is contained between -1 (for the energy conservation mode) and -2 (for the momentum conservation mode). The particle speed decrease with distance is best seen in Fig. 11b for t = 0.08 ms, and the power-law exponent can be estimated from a fit on the fastest particles for a range of distances. We consider a range of distances with enough particles to be representative of the wave propagation, but not too far from the impact point to consider only particles directly set in motion by the wave. We find a slope of about -1.91; therefore, in our simulations, the impact wave propagation mostly exists within the momentum conservation regime. Wada et al. (2006), in their SSDEM simulations with similar impact speeds, found that their impacts belonged to the same regime.

We find similar results with smaller impact speeds, although the particle speeds are obviously slower with smaller impact speeds. Power-law exponents are also found to be very close to -2, and indicate a pressure wave propagation in the momentum conservation mode.

Comparison with Z-model

We also compared the first instants of our simulations to the Z-model, which is an analytical model that represents the ejecta velocity field with streamlines after a vertical impact [START_REF] Maxwell | Impact and Explosion Cratering: Planetary and Terrestrial Implications[END_REF]). If we assume that the material is incompressible, the geometry of the streamline can be expressed in two dimensions (taking advantage of the axial symmetry) by a simple equation, in polar coordinates (r, θ), where the center of the coordinates system is the impact point, and θ = 0 corresponds to the downward vertical direction. The equation of a given streamline is:

r = R(1 -cos θ) 1 Z -2 , ( 7
)
where R is the horizontal distance from the impact point to the intersection between the considered streamline and the surface, and Z the main parameter of the Z-model, a decay parameter that defines the curvature of the streamlines. An example of streamlines is shown in Fig. 12, with the ejection angles.

From Eq. 7, the ejection angle of particles θ e can be determined. This ejection angle is supposed to be the same whatever the horizontal distance R, and equal to:

θ e = arctan (Z -2) , (8)
With the same method used for Fig. 2, we considered toroids with 0.25 cm-side square sections and looked at the average particle velocity in each of these toroids. We were more interested by the velocity direction than the magnitude, and for a clear figure, we divided the velocity vectors by their norm to look only at their orientation.

In order to find the corresponding Z for the beginning of our simulation, we consider the ejection angles assumed to not be influenced by the boundary conditions. We also do not account for the vicinity of the center where there are very few particles and therefore where variations in the ejection angles are much larger. Thus we consider ejected particles in a region that is at a distance from the center from 4 to 10 cm and in the upper layer (height from -0.25 to 0 cm). This gives us an average ejection angle of about 44 • , and Z = 2.966. However, with such a value for Z, data do not seem to match the streamlines defined by the Z-model in deeper layers.

To improve our determination of Z, we do not consider only the uppermost layer but also deeper ones (from about -6.5 to 0 cm). From Eq. 7, the velocity angle θ v (with the horizontal) for any point in the bed along a streamline can be determined:

θ v = arctan (Z -2) sin θ -(Z -1) sin θ cos θ 1 + (Z -2) cos θ -(Z -1) cos 2 θ . ( 9
)
If we consider the ejection angle, i.e., θ = π 2 , we find that Eq. 9 is compatible with Eq. 8, namely θ v (θ = π 2 ) = θ e . From Eq. 9, we can compare each average velocity angle per section in the considered area to the theoretical ones given by the Z-model. We determine the corresponding R and θ for each center of the areas from Eq. 10 and 11, and compare the velocity angle in our simulation with the theoretical one for a varying Z. We choose the Z value on the basis of the minimum root-mean-square of the differences Particle speed as function of the distance from the impact point for a 300 m s -1 impact, and for two different times, 0.012 ms after impact on the left, and 0.062 ms after impact on the right, corresponding respectively to t = 0.03 and t = 0.08 ms after the start of the simulation. The dotted lines represent the slope of the fastest particles for different distances, in orange for t = 0.03 ms, and in red for t = 0.08 ms. The orange slope has a value of -1.792 and the red one -1.91. We only consider particles within the vertical 60 • half-apex-angle cone under the impact point.

between theoretical velocity angles and numerical ones.

θ = arctan x area -z area , (10)
R = x 2 area + z 2 area (1 -cos θ) 1 Z -2 (11)
where (x area , z area ) are the coordinates of the center of each 0.25 cm square area, in the plane (x area corresponding to the distance from the cylinder central axis, and z area to the height). This leads to Z = 2.675, with an average angle error of 8.9 • and a standard deviation of 6.8 • . The streamlines and the velocity angles are shown in Fig. 13 for the computed Z.

We see that the deeper and the further from the impact point, the larger the angle differences. Since the Z-model is supposed to be for a stationary state, we looked at other instants much further in time, 1.25 • 10 -3 s and 1.25 • 10 -2 s after impact, and try to associate a Z to those instants. The further in time we consider, the greater the influence of the boundary conditions on the velocity distribution, as the cavity expands.

For 1.25 • 10 -3 s after the impact, the Z-model works better than just after the impact, and for the area considered, we find that the best match is for Z = 2.836, with an Figure 13. View of the average normalized velocities in the bed as function of their distance from the center of the cylinder and the height, at time t = 4.52 • 10 -5 s after the impact. Only a part of the cylinder is shown here. Each red arrow corresponds to the average normalized velocity for particles inside the corresponding 0.25 cm-side square section toroid. Blue lines correspond to the streamlines associated with Z = 2.675 (see text) for different R separated from each other by 5 cm, and grey arrows correspond to the normalized velocities associated with these streamlines, in the center of the toroids.

average angle difference between the model and our data of 3.6 ± 3.2 • (see Fig. 14).

When we try to fit our data to the Z-model further in time, i.e., for 1.25 • 10 -2 s, we see that the discrepancies increase. The best fit is for Z = 3.388, with an average angle difference of 9.8±8.7 • , which is much larger than the previous one (see Fig. 15).

The change in Z as a function of time can be explained by the fact that an almost stationary state as assumed in the Z-model is only reached near 1.25 • 10 -3 s after impact. The angle errors are indeed much higher (more than twice) if we consider a time after impact that is too short (stationary state not reached yet) or too long (effects of boundary conditions). This implies that the Z parameter associated with our impact is close to 2.836. According to Housen et al. (1983), momentum conservation is associated to a Z of 3, whereas energy conservation corresponds to Z = 2. In our simulations, Z is much closer to 3, and therefore to the momentum conservation regime, than to 2. This is consistent with the regime we deduced previously from the application of the shock propagation theory to the particle speed in the bed.

We also notice that the best fit for Z increases with time. An increasing Z with time for a normal impact was previously observed in other numerical models [START_REF] Austin | Multi-ring basins: Formation and Evolution[END_REF]. Studies from Cintala et al. (1999) and [START_REF] Anderson | Lunar and Planetary Science Conference[END_REF] also suggested, from the change of ejection angles and positions with time, that either the value of Z or the depth of the flow-field center (starting point of all streamlines) should change over time during the cratering process.

A migrating flow-field center was observed in numerical models [START_REF] Thomsen | Lunar and Planetary Science Conference Proceedings[END_REF][START_REF] Austin | Lunar and Planetary Science Conference Proceedings[END_REF]) and generalized in Croft (1980). With a flow-field center not at the impact point, the ejection velocity depends on the horizontal distance R. Eq. 10 and 11 need to be changed into Eq. 12 and 13:

θ = arctan x area H 0 -z area , (12)
R = x 2 area + (H 0 -z area ) 2 (1 -cos θ ) 1 Z -2 (13)
where θ and R correspond to previously defined θ and R but from the vertically translated flow-field center, as shown in Fig. 16. Therefore, R corresponds to the horizontal distance for the intersection between the streamline and y = H 0 . By replacing R and θ into R and θ in Eq. 9, we can determine the velocity angle in ((x area , z area ), and use the same method as before to find the best matches for H 0 and Z.

With the new addition of a migrating flow-field center, we look at 1.25 • 10 -3 s after impact, when we had the best fit of the three considered instants. We find that the best fit is Z = 2.857 and H 0 = 0.05 cm. H 0 is very low and we can consider that the best fit was in fact obtained for a flow-field center at the impact point. When we try to fit (Z, H 0 ) for 1.25 • 10 -2 s after impact, we find that Z should be higher than 5, which, according to Kurosawa & Takada (2019) should be the upper limit for Z. This means that the Z-model is not adapted to the geometry of the streamlines that we find for this impact instant. As previously stated, this may be due to the crater not being in the stationary excavation flow phase anymore, and this is potentially due to boundary condition effects. This may also be due to the complexities of a granular medium.

Ejecta speed, volume, and direction

After having described the cavity formation and characterized the impact by its first instants, we look at the ejecta, i.e., particles ejected from the regolith bed. In this study, we consider as ejecta target's particles whose heights are higher than the maximum height of surface particles before the impact. Concerning the ejecta speed, we analyze their dependency on the ejection position, i.e., the distance between the impact point and the position from where the particle is ejected. According to Housen et al. (1983), from dimensional analysis and laboratory experiments, Equation 14 links the normalized ejection speed and the normalized ejection position:

v e √ gR = 0.62 d R -2.55 (14)
where v e is the ejection speed, d the ejection position, R the crater radius, and g the gravitational acceleration. Wada et al. (2006) also derived such a law from numerical simulations, but using the transient crater radius and not the final one:

v e √ gR = 0.923 d R tr -2.50 (15)
where R tr is the transient crater radius.

In our simulations, we reach neither the transient crater nor the final one. To normalize the ejection speed and distance, we use instead a characteristic length of our simulation, i.e., the cylinder's radius R cyl = 15 cm, which leads to Fig. 17. We used the data at time about 1.25 s after impact because it corresponds to the end of the main increase of the cavity depth. Later than 1.25 s after impact, the cavity bottom is more or less stabilized, and the depth is roughly constant.

We see in Fig. 17 that there are several particles with relatively small speed and very close to the impact point. These are particles near the surface, not directly hit by the impactor but barely in contact with other ejected particles, which felt a very small upward acceleration and are therefore considered as ejecta. We also see that the further from the impact point, the slower the ejecta, as expected. We also see that in our simulations, the ejecta do not behave as predicted by Housen et al. (1983) and Wada et al. (2006). There are much more low-speed particles far from the impact point than predicted with a power law like Eq. 14 and 15. This may be due to the earliness of the crater formation, to the boundary conditions, or to the material properties. Furthermore, although the gravitational acceleration is taken into account in the normalization of the ejection speed, since power laws from Housen et al. (1983) and Wada et al. (2006) come both from Earth's gravity environment, the significant difference in the value of g may also be an explanation for these dissimilarities.

However, despite these variations, data from this study are not so far from previous data. We can try to derive a final crater radius, or a transient crater radius, from our data. As shown in Fig. 17, if the power laws from Wada et al. (2006) and Housen et al. (1983) are applied to our simulations, we expect a transient crater of about 3 -3.5R cyl and a final crater of about 4 -4.5R cyl . This confirms that boundary conditions certainly play a role in our simulations for the crater dimensions, at least when the cavity radius gets close to the cylinder radius.

We observe similar behaviors for lower speeds, leading to smaller crater radius equivalences. The slower the impact, the smaller the final crater, at least from our equivalences explained previously. We find that our data for 50, 100, and 200 m s -1 correspond to final crater radii of respectively 2.5 -3R cyl , 3 -3.5R cyl , and 4R cyl .

We also study the amount of ejecta as a function of time, still without considering the presence of the sampling horn. The ejecta volume as a function of time is shown on Fig. 18a for four different impact speeds. We ran these simulations for scaled timesteps and spring constants proportional to impact speeds, as well as with the same timesteps and spring constants for all of them to be sure that the choice has no big influence on the outcome. We comfirm the absence of significant influence (in the short range we covered), as expected from Sec. 3.1.2. The ejecta volume clearly depends on the impact speed: a higher speed leads to more ejecta, if we look at a certain time after impact, or same amounts of ejecta are reached sooner for a higher impact speed. This is what was expected, as higher speeds mean higher incoming energies to be distributed to the target.

By scaling the time with the impact speed (see Fig. 18b), we notice that the relation previously indicated is now directly a proportionality. If the impact speed is twice higher, it takes half the time to reach the same amount of ejecta. Indeed, if we remove this time difference, we have very similar amounts of ejecta. This could be expected when we scale the timesteps and the spring constants with the impact speed, since the main force on the bed is the impactor's penetration force, but we also have this result when we set the same timesteps and spring constants for all speeds. This is particularly due to the very low gravity, allowing particles with even small speeds to escape from the bed. Indeed, we expect particles to have higher ejection speeds with a 300 m s -1 impact than a 50 m s-1, but a great amount of particles in both cases is ejected in this low-gravity environment.

The ejected volume can also be expressed as a function of the ejection speed. Here we consider the total volume of all particles faster than a given speed, normalized by a characteristic length, as a function of the normalized ejection speed. Due to the transient crater not being reached and the boundary condition effects (a too small cylinder), we are missing ejecta. Therefore, we do not expect to have as much ejecta volume as in Housen et al. (1983) and Wada et al. (2006), even more since the ejected volume in Housen et al. (1983) comes from measurements of ejecta blankets, and therefore corresponds to a porous global volume and not a volume computed from the sum of grains. However, by using the equivalent final crater radii given previously in this section, we can at least check the agreement of the slope, shown in Fig. 19. Since we saw that it takes twice the time to have the same amount of ejecta for an impact with half the speed, we considered each simulation at different times, scaled according to the impact speed in order to represent the same stage in the crater formation.

We observe in Fig. 19 similarities with Wada et al. (2006), i.e., a depletion of low-volume particles, certainly due to the too small number of grains in our simulations, and a decrease of volume with an increasing ejection speed. We also notice that Wada et al. (2006) had a higher volume of slow ejecta but a much lower amount of fast ejecta. This may be due to the difference in the material properties, as the energy is not necessarily transmitted with the same efficiency through the bed in both simulations.

According to Housen et al. (1983), the power law that should apply to the normalized volume of ejecta faster than a given speed is:

V(> v e ) R 3 = 0.32 v e √ gR -1.22 ( 16
)
where V is the volume, v e is the ejection speed, R the crater radius, and g the gravitational acceleration. This power law is shown in Fig. 19 and is as expected higher than for our data. Also, the section we obtain with a clear slope is for particles faster than the ones used for deriving the power law in Eq. 16. Nevertheless, our slope is very close to the theoretical one.

Finally, to conclude this study on the ejecta, we analyze the ejecta velocity directions. The prediction of the ejection direction of particles can be useful for a spacecraft shooting on a low-gravity surface, both for risk minimization or sampling efficiency.

For different times, we compute the direction of the velocity of each ejecta. We can then deduce the preferred direction of ejection. Since the setup shows axial symmetry, directions are represented by angles ranging from 0 • (horizontal plane) to 90 • . The ejecta velocity directions for a 300 m s -1 impact are shown on Fig. 20, for a duration up to 1 s. Since there is gravity, even if low, the velocity considered for computing the angle is the original one when leaving the bed, and not the current one. This means that variations with time are due to new particles being ejected, and not to gravity affecting the ejecta speed directions. The ejecta velocity angle therefore corresponds to an ejection angle.

During the whole simulation, most particles are ejected with an angle of about 50 • , and in majority between 48 • and 54 • . This means that a spacecraft located above for example 70 • should be safe from any ejecta, or that, in the opposite case when the spacecraft wants to capture particles, it needs to cover angles higher than about 45 • to increase its chances of sampling. As a comparison, Wada et al. (2006) also found a majority of ejecta having an ejection angle around 45 and 50 • , with a same impact speed and a projectile slightly smaller (radius of 3 mm instead of 4 mm).

We also studied the ejection angles for lower impact speeds down to 50 m s -1 and found that the preferred direc- Volume V( > ve)/R 3 50 ms at 6 t0 100 ms at 3 t0 200 ms at 2 t0 300 ms at t0 Housen et al. 1983 Figure 19. Normalized total ejecta volume of particles whose speed is larger than a given speed as a function of normalized ejection speed, where R is the radius of the equivalent final crater (see text for details). We considered R 50m/s = 2.75R cyl , R 100m/s = 3.25R cyl , R 200m/s = 4R cyl , and R 300m/s = 4.25R cyl . The data correspond to different times depending on the impact speed to reach the same crater phase, as shown in the legend, in which t 0 = 0.6 s. The line labeled Housen et al. 1983 comes from Eq. 16. tions seem to be similar whatever the impact speed. However, by looking closely at the results, the preferred ejection angle slightly increases with the impact speed. For example, for a 50 m s -1 impact, there are much more particles ejected with angles between 42 • and 48 • , even if the majority is still in the previous angular section. This is due to the impactor penetrating deeper into the bed for higher impact speeds, and therefore ejecting particles with higher angles. Now that we have characterized the impact, from the cavity's dimensions to the ejecta volume, we add the sampling horn in our simulations to apply these results to Hayabusa2 sampling.

Application to Hayabusa2 sampling

The sampler of Hayabusa2 is composed of a long horn with a non uniform radius, a filter located near the top of the horn, and a narrow path to the sample return capsule. Its geometry is described in Sawada et al. (2017). We modeled the horn and a part of the sampler storage and transfer mechanism (that we incorporate to the horn for more simplicity in this study) thanks to five circular frustums, i.e., cylinders with tapered radii, or truncated cones. The four pieces of the filter were also modeled faithfully to the original one, and prevent particles with a diameter larger than 1 cm to go through, in order to avoid any jamming upstream from the sample return capsule. The bottom of the horn is supposed to come in contact with the surface, and the projectile is fired immediately after the contact with the surface. However, to model the contact between the horn and the surface would require the knowledge of all moments of inertia, including the spacecraft, and to implement the reaction of the whole system to the surface. In order to simplify, we consider the spacecraft and the sampling mechanism as hovering motionless 0.5 cm above the surface, and we are more interested in particles ejected from the projectile's impact than due to the contact between the horn and the surface.

In the actual sampling mechanism, particles in the top cylinder (later called Cylinder 5) are transferred to the storage area, located at the same height as the cylinder but not represented here for more clarity. We did not seal the top of the sampler in order to avoid jamming inside Cylinder 5, where there would not be any in reality. However, when we consider the mass collected, we take into account particles that went through the top of Cylinder 5 as well as particles still inside it.

General results

First, we consider a simulation of the 300 m s -1 impact, with a nominal timestep of dt norm = and a nominal spring constant of k n , i.e., the same ones we took for our 300 m s -1 impacts in Section 3.1. A snapshot of the simulation is shown in Fig. 21a, where each color represents a different cylinder, i.e., a different part of the sampler. Since the sampler is supposed to stay about 1 s at the surface before activating thrusters and leaving the asteroid's surface, we represent here the state of our simulation 1 s after the impact.

In order to check what could happen if we stay longer on the surface or what could be seen from the spacecraft, we ran some simulations with a larger timestep, for which we adapted the spring constant with the rule proposed by Schwartz et al. (2012) as explained in Section 3.1.2. This allows us to simulate the sampling up to 4.5 s, as shown on 21b. We observe that, even if we have a bed limited in size, we already have a lot of ejecta that can also be seen from the spacecraft as they are flying outside the sampling horn. That was indeed observed by Hayabusa2 on-board small monitor camera (CAM-H) right after the sampling. We also see that particles all along the simulations tend to gather together inside what we called Cylinder 4, i.e., the frustum that contains the filter. It is due to a combined effect of the filter and the bottleneck-like geometry of this part of the sampling horn. This is confirmed by a representation of the packing fraction inside the sampling horn, as shown in Fig. 22, also for 1 and 4.5 s after the impact. The packing fraction is defined here, for each 0.1 cm horizontal slice of the sampling horn and the bed, as the volume occupied by particles divided by the volume of the slice. We see that the packing fraction increases with time near the filter, as particles accumulate. However, we see that the packing fraction does not exceed 0.3, which implies that the part upstream the filter is still not jammed after 4.5 s. Indeed, we see peaks of packing fraction inside the upper cylinder on both panels of Fig. 22, as well as particles in Cylinder 5 in Fig. 21. We also see in Fig. 22 that the packing fraction in the bed drops off as particles are ejected and the crater forms.

There is stochasticity in our simulations. Even with exact same initial conditions, a very slight change in a calculation at the moment of impact can have consequences on the direction of the ejecta, and thus on the direction they impact the walls of the sampling horn. These errors can accumulate and lead to different particles density in the cylinders. Naturally, in the highest cylinders, the particle density is smaller, and a slight difference at impact can imply a particle not going through the bottleneck or the filter. Therefore, the number of particles collected can vary from one simulation to the other. The relatively low number of particles and associated large particular size in our simulations consequently also leads to a low number of particles that go through the filter, which emphasizes this effect. We did not consider very small particles, which would be ejected faster by the impact, and would have more chances to go through the filter. Therefore, if smaller grains are present on Ryugu, our simulations could represent a worse case scenario since from them, we should expect to collect less particles than the actual sampling. What Nevertheless, our simulations can help to determine the expected ratio of particles between the different cylinders.

We used data from two different simulations with different timesteps and spring constants for Fig. 21 and Fig. 22, and we saw in Sec. 3.1.2 that the energy distribution differs depending on these parameters. However, there is no noticeable influence on the ejecta volume (see Sec. 3.1.5). In order to confirm the lack of influence on the sampling, we ran several simulations with different ∆t/k n pairs, and we show the distribution of particles in each cylinder in Fig. 23. This is also a way to reduce the stochasticity by running more simulations.

First, we see that the ∆t/k n pairs have no major influence on the outcomes of the simulation, as long as they stay in an appropriate domain where the simulations do not crash and that the first contact between the projectile and the surface lasts enough iterations. There are dissimilarities between the simulations, but there are as large between different ∆t/k n pairs as they are for same pairs, and are therefore due to stochasticity. Nevertheless, even in Cylinder 4, the number of particles in all simulations is almost equivalent. The largest discrepancies appear for Cylinder 5, the one where particles are collected, where the number of particles is the smallest, and therefore where every particle matters.

We see that the peak of particle number moves with time among the cylinders as particles go up in the sampling horn, and that the bed does not replenish the lower cylinders with particles during the whole simulation. However, we see that after 4.5 s, Cylinder 2 still contains a huge amount of particles, almost a hundred, and this is a reminder of the lowgravity environment, where gravity has a very slow influence on ejected particles.

For almost all of the simulations, at least one particle would be collected after 1 s, i.e., has penetrated into the final cylinder. The collected mass after 1 s can go up to 0.4 g, which is higher than 100 mg, the scientific objective of Hayabusa2 (Sawada et al. 2017). The density of particles in our simulations being 2.43 g cm -3 , one particle with a radius larger than 0.215 cm or two particles with any sizes are enough to fulfill the scientific requirements. The mean radius in our particle distribution is 0.25 cm, and since we get more than one particle in most of the simulations, the scientific goal is usually fulfilled.

After 1 s, the density is still increasing in Cylinder 4, upstream to the filter, meaning that a longer time eventually leads to more collected particles, if a jamming does not occur. As we previously said, the number of particles collected in our simulations is not necessarily representative of the actual one, as smaller particles would be ejected with a higher speed. Also, the contact between the surface and the toothed bottom of the horn could also increase the volume of ejected material. Despite these dissimilarities, we still have collected particles for most of the simulations, which is auspicious for the actual sampling.

In order to continue the comparison with lower speeds as it was done in Sec. 3.1, we did simulations with 50 m s -1 impact speed. As it can be expected, the ejecta take much more time to reach the highest cylinders, and by 1.5 there are less than 50 particles in Cylinder 3 (and none has reached Cylinder 4), whereas for 300 m s -1 impacts, we have more than 300 particles. We saw in Sec. 3.1.5 that a six times lower speed implies a six times larger duration to reach the same amount of ejecta. We interestingly find the same correlation in each cylinder (at least before Cylinder 5, in which the number of particles is too stochastic). Indeed, if we look at the amount of ejecta for the 50 m s -1 impact in Cylinder 4 after 1.2 s for example, it is close to the amount of particles in the same cylinder for the 300 m s -1 impact after 0.2 s. Peaks of population also follow this rule in lower cylinders. The volume of particles inside the cylinders seems to scale linearly with the impact speed.

t = t nom , k n = k n, nom t = 2 t nom , k n = 1 4 k n, nom t = 4 t nom , k n = 1 16 k n, nom t = t nom , k n = 0.56k n, nom t = 1 2 t nom , k n = 0.56k n,

Second projectile

Since Hayabusa2 spacecraft has three projectiles at its disposal (Sawada et al. 2017), it is interesting to check if firing the second projectile a short moment after the first one during the same sampling could significantly increase the number of ejecta and collected particles, or if the changes would stay marginal. In order to investigate this, we model the impact of a second impactor 0.2 s after the first one. The projectile is fired from the bottom of Cylinder 4.

The number of particles per cylinder for simulations with and without a second shot is shown in Fig. 24.

As expected, firing a second impactor increases even more stochasticity. In effect, discrepancies between simulations increase as the second projectile does not hit the same surface for each simulation, and can even hit ejecta before impacting the surface. In all simulations, firing a second projectile significantly increases the ejecta volume in lower cylinders. However, for upper cylinders, beginning with Cylinder 3, a second shot can have the effect of reducing the number of particles or at least delaying its growth. Also, it is noticeable that increasing the particle volume in Cylinder 4 does not necessarily increase the number of particles in Cylinder 5, and therefore the number of particles collected. Increasing the ejecta number also increases collisions between particles, and therefore can reduce the speeds of particles or change their directions. Now that we have characterized the impact, and looked at the outcomes for such an impact in the context of Hayabusa2 sampling, we try to expand the parameter space to see if some of the parameters we initially chose have a significant influence on results of Sec. 3.1 and 3.2. In this section, we investigate the influence of the target's normal and tangential coefficients of restitution ε n and ε t on the impact and the sampling. We begin with the influence of ε n on the energy distribution during the first instants. The influence of this coefficient has already been studied in [START_REF] Thuillet | [END_REF] concerning the interaction of the Hayabusa2 French-German (CNES-DLR)landing package MASCOT with Ryugu's modeled surface. It was found that a lower ε n leads to a harder surface and therefore to a higher outgoing-to-incoming speed ratio for MASCOT.

t = t nom , k n = k n, nom t = t nom , k n = k n, nom (2 nd impact) t = 2 t nom , k n = 1 4 k n, nom t = 2 t nom , k n = 1 4 k n, nom (2 nd impact) t = 4 t nom , k n = 1 16 k n, nom t = 4 t nom , k n = 1 16 k n,
First, we noticed that the spring constant had to be lowered for simulations with a very low normal coefficient of restitution. Otherwise, the simulation would crash because particles went too far from each other in a too short duration, i.e., they would reject each other too strongly. However, one could think that decreasing the normal coefficient of restitution would be not equivalent but at least similar to decreasing the spring constant. Indeed, according to Schwartz et al. (2012), the oscillation half-period of an isolated twoparticle collision τ overlap in pkdgrav can be computed with Eq. 17.

τ overlap = π ω 0n 1 -ξ 2 (17)
where ω 0n = k n µ is the undamped harmonic frequency, µ is the reduced mass of the colliding pair, and ξ is the damping coefficient that can be computed from the normal coefficient of restitution ε n thanks to Eq. 18.

ξ = -ln ε n π 2 + (ln ε n ) 2 (18)
Equations 17 and 18 lead to the expression of τ overlap as a function of k n and ε n shown in Fig. 19.

τ overlap = µ π 2 + (ln ε n ) 2 k n (19)
In Eq. 19, if ε n decreases, since 0 < ε n < 1, (ln ε n ) 2 increases, and we could expect overlap durations to increase, and to find a similar behaviour as if k n would decrease. However, this is only true for a single contact, and it is much more complex for a bed composed of thousands of particles. For example, we do not see a 121% increase of overlap duration as predicted when changing ε n from 0.5 to 0.1, even if we see slight differences in the overlaps between particles. If we look at the energy distribution at the first impact instants as we did in Sec. 3.1.2, when ε n reaches very small values like 0.1, we see a very high amount of energy, much higher than for any other parameters. The amount of energy in the projectile stays more or less the same but the target's rotational energy can reach very high values, and the energy after 1.2 • 10 -4 s is about 3.75 times the one for ε n 0.2. If 0.2 ε n 0.5, or 0.5 ε n 0.9, we do not see any significant difference in the energy distribution. We find that this energy surge effect for very low ε n disappears for smaller timesteps, for example with a timestep two times smaller. This phenomenon is nonetheless interesting because, from Eq. 19, a decrease in ε n should increase the overlap durations, and therefore the timestep could be increased. On the contrary, it has to be decreased to faithfully represent the impact mechanism. Once again, we find an unexpected behavior for the normal coefficient of restitution.

Once the timestep is small enough, we find that the distributions are very similar whatever ε n . The main differences are that the higher ε n , the higher the peak of target's kinetic energy and the higher this kinetic energy after 1.2 • 10 -4 s. If we look at the target's kinetic energy after 1.2 • 10 -4 s for ε n = 0.1 and ε n = 0.9, these variations account for about 6% of the initial total energy, which can seem low, but represents a change from 4.3% to 10.4% of the total initial energy, and therefore this energy more than doubles from a very low ε n to a high one.

The amount of ejecta does not seem to depend on ε n . By looking at the number of particles located at higher levels than the initial surface like we did in Sec. 3.1.5, we find the same amounts for ε n equal to 0.1, 0.5, and 0.9. Therefore, even if a higher normal coefficient of restitution leads to slightly larger target's kinetic energy, this does not affect the ejecta amount. This is true only if the timestep is small enough, otherwise we find that the amount of ejecta is much larger for ε = 0.1.

Concerning the direction of the velocities of these ejecta, the higher ε n , the smaller the ejection angle (between the horizontal and the velocity vector). Most particles are Volume V( > ve)/R 3 cyl ( n, t) = (0.5, 0.1) ( n, t) = (0.5, 0.5) ( n, t) = (0.5, 0.9) ejected with an angle between 48 and 54 • for ε n = 0.1, whereas for ε n = 0.9 most ejecta velocity angles are in the range between 42 and 48 • .

With a higher ε n , particles are also ejected with higher speeds, as it is shown in Fig. 25. We see that the volume of ejecta for very low ejection speeds are very similar in all cases, but with high ε n , we have a slightly higher number of high-speed particles. Therefore, even if we saw no direct correlation between the amount of ejecta and the energy distribution during the first instants, we see that the higher target's kinetic energy comes from the fastest ejecta.

When ε n is higher, the impact wave feels less damping and spread more easily inside the regolith bed. This leads to particles ejected with a higher speed. However, not only ejected particles get a higher speed, but also particles inside the bed, and this implies a visible difference in the cratering process. Indeed, both the lateral and vertical growths of the crater are faster with high ε n , as shown in Fig. 26. If we look at the depth-to-radius ratio, it is almost the same for ε n = 0.1 and ε n = 0.5, i.e., about 0.68 after 0.25 s, but is much higher for ε n = 0.9 and is about 0.78 at the same time. Yamamoto et al. (2006) found that cratering depends on, among others, internal friction, and our results also provide information on how cratering depends on the physical properties of the target medium.

Concerning the tangential coefficient of restitution ε t , we see higher target's kinetic energies as well with higher ε t , but these do not have the same outcomes as for ε n . Indeed, ejecta volumes and ejecta speeds are very similar whatever ε t . Also, as shown in Fig. 26, the relation between ε t and cavity radius or depth is not as simple as it is with ε n . Even if the cavity does grow faster with higher ε t at the beginning of the impact (which was expected from mapping the energy distribution at the first instants), it changes with time and for example after about 0.15 s, a low ε t can lead to a deeper cavity than found with a high one. Therefore, the relation between ε t , crater dimensions, ejecta volume, and ejecta speed is a more complex one than for ε n , and no clear trend can be established.

We also checked the influence of the normal coefficient of restitution ε n on the sampling (see Fig. 27). We saw that ( n , t ) = (0.1, 0.5) ( n , t ) = (0.5, 0.5) ( n , t ) = (0.9, 0.5) ( n , t ) = (0.5, 0.1) ( n , t ) = (0.5, 0.9) ( n, t) = (0.1, 0.5) ( n, t) = (0.5, 0.5) ( n, t) = (0.9, 0.5) ( n, t) = (0.5, 0.1) ( n, t) = (0.5, 0.9) Figure 26. Cavity dimensions (depth and radius) in cm as a function of time for various simulations, differing from each other in their coefficients of restitution ε n and ε t . All represent a 300 m s -1 impact. a lower ε n implies higher ejection angles but particles are ejected with smaller speeds. In Fig. 27, we observe that ε n = 0.9 leads to high particle densities in Cylinders 1 to 4 for the first tenths of seconds, due to higher ejecta speeds and to a high coefficient of restitution on the walls constituting the sampler horn. With a very low coefficient of restitution (ε n = 0.1), particles are ejected more slowly and therefore the particle density inside the cylinders takes more time to increase. However, since the ejection angle is higher, the particle density increases more than for higher ε n in lower cylinders. Because of the low coefficient of restitution on the walls, it is more difficult for particles to reach the highest cylinders. With a small ε n , particles that succeed in reaching Cylinders 4 and 5 are ejected with higher angles and bounce less against the sampler's walls than with a high ε n . We find the same results with a 50 m s -1 impact, which confirms that particles are ejected faster with a high coefficient of restitution but that ejecta inside the horn are less numerous in total.

Since the amount of ejecta seemed to be the same for all ε n , the low number of ejecta with high ε n inside the horn is certainly due to the crater forming faster, and therefore the cavity's boundaries quickly go over the mouth of the sampler horn. Hence, particles are still ejected but outside of the horn.

Influence of the impact geometry

Since the regolith bed cannot be perfectly flat, being composed of grains, it is interesting to look at the influence of the impact position with respect to the surface grains on the energy distribution, the crater formation, and the sampling.

First, we look at the energy distribution in the case of a slight shift along the x-axis. Since the average particle radius is 0.25 cm, we chose this value for the shifts, in order to be sure to have a different layout at the impact location than the one already tried. If we call x 0 the regular impact location along the x-axis (the cylinder's axis is located at x = x 0 and y = y 0 ), the energy distribution and impactor-particles overlaps for shifts of x = x 0 ± 0.25 cm are shown on Fig. 29 and28.

We see that the energy distribution at the first instants of the impact depends on the impact location. On Fig. 28, we see that the projectile comes in contact with three particles at the same time, i.e., it goes into a "hole" in the middle of these particles. That is why the impact happens later than for the other simulations, because the impact location height is lower. Also, since it is stopped by three particles at the same time, we see that the overlaps are shallower and the projectile's kinetic energy decreases very fast. In this case, we do not have two phases as in the regular impact location. Nonetheless, the target's kinetic energy after 1.2 • 10 -4 s is very similar in both cases.

With a -0.25 cm shift (see Fig. 29), the projectile hits the same top particle as with the regular location but deeper, and this time there are no other very close particles. Because of that, the projectile's kinetic energy does not decrease as fast as with previous impact locations. The projectile eventually hits other particles and slows down. The elastic energies and the overlaps show that the first contact is longer than for example the +0.25 cm case. We see that, even if the energy distribution depends on the geometry of the surface where the impact happens, there is always a fast decrease of target's energy, peaks in elastic energies and target's kinetic energy. Moreover, the kinetic energy contained in the grains after 1.2 • 10 -4 s is of the same order of magnitude.

Concerning the whole crater formation, the cavity depth and radius do not seem to depend on the impact location, and the small disparities are more due to stochasticity than to actual differences in the formation mechanisms (see Fig. 30). When we look at the ejecta volume as a function of the distance from the impact point, the shift in the impact point seems to have no influence. We can see some dissimilarities when looking at the ejecta volume as a function of the ejection speed, but these faint variations do not correspond to faster cavity growths, and are therefore not as noticeable as they were for ε n in Sec. 3.3.1.

Because the impact may not be perfectly vertical to the surface (and will certainly not be as the projectors are not parallel to the horn central axis), we also study the influence of a tilted impact, with an angle of -25 • or +25 • from the vertical in the x-z plane. Whereas the evolution of the cavity radius is very similar whatever the geometry of the impact, the cavity is significantly shallower in the case of a +25 • tilted impact (see Fig. 30). When we also look at the amount of ejected particles, we see that a +25 • tilted impact generates much fewer ejecta than the other geometries. A -25 • tilted impact also produces less ejecta than a normal impact during the first instants but the volume increases faster than for other geometries, and the volume of ejecta becomes equal for all geometries except the +25 • case at about 0.8 s after the impact. This shows that the regolith bed is anisotropic, and this is certainly due to the layout of the grains and to the force chains in the bed itself. In the vertical direction, particles are subject to very different constraints depending on the particle layout, and it can be easier or harder for them to move down, whereas in the horizontal direction, the constraints are weaker due to the surface, and therefore these variations of constraints are fainter. We also find that the directions of ejection do not depend on the impact angle within the considered range, and that most of the particles are ejected with an angle between 48 • and 54 • whatever the impact angle.

Looking at the high-speed ejecta, we see that they are faster and more numerous from a tilted impact than from a vertical one (see Fig. 31). Indeed, for both simulations with a tilted impact, we find a larger population of ejecta with an ejection speed higher than 10 m s -1 . This is in agreement with [START_REF] Yamamoto | [END_REF], who also find an increase in high-velocity ejecta with an increasing impact angle (from the vertical).

Simulations including the sampling horn show a shortage of ejecta observed for the +25 • tilted impact, which has an influence on the particle density in Cylinders 1 to 4, as shown in Fig. 32. Surprisingly, more particles go through the filter and into Cylinder 5. This is certainly due to the high stochasticity of the particle number in Cylinder 5; the fact that a particle penetrates or not into Cylinder 5 depends on very slight differences in the particle direction. It may also be due to the lower particle density in lower cylinders, and therefore fewer collisions between particles. More particles can then reach Cylinder 5 without experiencing any collision with other ejecta. We also observe a temporary shortage of ejecta at the beginning of simulations of the -25 • tilted impact. In summary, tilted impacts produce less ejecta during the first instants, which can be seen in the lowest cylinders. However, that does not mean in our case that the number of collected particles is smaller. The medium appears to be anisotropic and the direction of impact can lead to a different number of ejecta.

We notice that the impact location has a small influence on the sampling. For a -0.25 cm shifted impact, particles ejected with a high speed are less numerous, leading to more particles in the lowest cylinders and less in the highest ones soon after the impact, but the particle density in these cylinders then increases with time as particles have time to go up. The impact location has therefore more influence on the ejecta speed than on the ejecta volume.

CONCLUSIONS

In this paper, we presented the results of our simulations of 300 m s -1 impacts into a regolith bed, under the low-gravity environment of Ryugu. We also considered for comparison impacts speeds of 50, 100, and 200 m s -1 . We first presented our simulations without the sampler horn, to characterize the impact and the reaction of our bed, without adding the complexity and influence of the sampling mechanism.

We found that the cratering process is much longer under Ryugu's low gravity than for 1 g simulations and experiments. At first, we see a hemispherical growth, then the depth increase slows down, leading to a mainly lateral growth phase. Higher impact speeds lead to larger cavities, but the evolutions of the cavity depth as a function of the cavity radius are very similar. We also studied the energy distribution during the first instants after impact, and we find that less than 10% of the initial kinetic energy is transmitted into the ejecta kinetic energy, which is consistent with what is typically assumed in collisional evolution models (e.g. [START_REF] Davis | Asteroids III[END_REF].

Confronting our simulations to the shock propagation theory, we find that the impact shows characteristics of the momentum conservation mode. Moreover, by applying the Z-model to our simulations, we find that it matches very well only for a certain time, when we have a stationary excavation state. When it matches, we find again that we are in the momentum conservation regime. From the three considered instants, Z seems to increase with time, but the flow-field center appears to be located at the impact point.

We also studied the particles ejected from the bed. Based on ejecta charcteristics, we can expect a transient crater between 3 and 5 cylinder radii, decreasing with the impact speed. Most of particles are ejected with an angle between 48 and 54 • .

When we consider the sampling horn of the Hayabusa2 sampling mechanisms, we observe some stochasticity as similar impacts can give different results, due to the high number of interactions between the projectile, the particles, and the horn. However, for most of the simulations, the scientific goal of 100 mg is almost always fulfilled after 1 s, even without taking into account the interaction between the horn and the surface, and the scoop-up part. Generally, we see that a higher number of ejecta in the lower parts of the horn does not mean a higher density in the upper part. As an example of this observation, a second shot 0.2 s after the first one would increase the initial amount of ejecta but not necessarily the portion that goes in the upper parts of the horn, due to a rising probability of collisions between grains.

Finally, we looked at the influence of the physical parameter ε n and of the impact geometry on our results. The influence of ε n is not as straightforward as basic equations could let us think, as we had to decrease the timestep for low ε n . The amount of ejecta does not seem to depend on this coefficient, but the ejection angle and speed do. With the sampling horn, we find that there are less particles in upper cylinders for a high ε n but they reach them more quickly. Our study on the impact geometry, with slightly translated and tilted impacts, shows that the regolith bed is anisotropic, and that the outcome of the impact depends more on the impact angle than on the exact location of the impact. With tilted impacts, high-speed ejecta are more numerous and faster, but in the first instants the total amount of ejected particles is usually smaller.

To conclude, this study provides many results concerning the outcome of impacts in a low-gravity environment, under conditions that are hard to reproduce on Earth. Our results are in agreement with many previous studies, even if we sometimes find discrepancies due to our relatively restraining boundary conditions.

This study is based on several assumptions (e.g., the size of the cylinder containing the bed, the grain sizes, the friction coefficients etc.) and in order to better understand the cratering process under low gravity, future numerical simulations should be done, for instance, with a larger cylinder to avoid boundary effects and enable the cavity to grow freely. To do so, solutions have to be found to decrease the computation time, or a large amount of time should be dedicated to run these simulations, because increasing the cylinder radius has a significant effect on computation time. Moreover, to better represent the actual rocky surface of Ryugu, grains could be represented as aggregates of smaller particles, that break when submitted to strong forces. Finally, the bottom teeth could be modeled, as well as the interaction of the horn with the surface, taking into account the inertia of the whole spacecraft, and making much more complex the simulations. This is something we will consider in future studies so that we full understanding of such a sampling mechanism efficiency can be assessed over a wide range of parameters and gravity conditions than can be achieved with experiments.

Figure 1 .Figure 2 .

 12 Figure1. Cross-section snapshots of the projectile impacting the regolith bed, in the yz-plane. In red are represented the velocities of the different particles in this plane. The cross-section width is 1 cm and the impactor is not represented in this section.

Figure 3 .

 3 Figure 3. Cavity dimensions (depth and radius) as a function of time, for a 300 m s -1 impact.

Figure 4 .

 4 Figure 4. Cavity depth as a function of cavity radius for four different impact speeds.

  Figure6. Zoomed-in cross-section representations of the projectile impacting the regolith bed, in the vertical plane containing both the initial position of the impactor and the first particle in contact, and for different times. The impactor is in green, and in red are represented the velocities of the different particles and impactor in this plane. The impact speed is 300 m s -1 , and the cross-section width is 1 cm. Distances are in centimeters.

Figure 7 .

 7 Figure7. Overlaps between impactor and different particles in contact with it during the simulation, for a 300 m s -1 impact.

Figure 9 .

 9 Figure9. Overlaps between impactor and different particles in contact with it during the simulation, for a 50 m s -1 impact, for k n = 0.58 k n 0 .

  Figure10. Zoomed-in cross-section representations of the projectile impacting the regolith bed, in the vertical plane containing both the initial position of the impactor and the first particle in contact, for different spring constants and a same time after impact. The impactor is in green, and in red are represented the velocities of the different particles and impactor in this plane. The cross-section width is 1 cm. Distances are in centimeters.

Figure 12 .

 12 Figure 12. Examples of streamlines (for R = 5, R = 10, and R = 15 cm) defined by the Z-model, with several values for the parameter Z. The ejection angle θ e is also shown at the surface for the considered Z.

Figure 14 .

 14 Figure14. View of the average normalized velocities in the bed as function of their distance from the center of the cylinder and the height, at time t = 1.25 • 10 -3 s after the impact. Only a part of the cylinder is shown here. Each red arrow corresponds to the average normalized velocity for particles inside the corresponding 0.25 cm-side square section toroid. Blue lines correspond to the streamlines associated with Z = 2.836 (see text for explanation) for different R separated from each other by 5 cm, and grey arrows correspond to the normalized velocities associated with these streamlines, in the center of the toroids.

Figure 15 .

 15 Figure15. View of the average normalized velocities in the bed as function of their distance from the center of the cylinder and the height, at time t = 1.25 • 10 -2 s after the impact. Only a part of the cylinder is shown here. Each red arrow corresponds to the average normalized velocity for particles inside the corresponding 0.25 cm-side square section toroid. Blue lines correspond to the streamlines associated with Z = 3.388 (see text for explanation) for different R separated from each other by 5 cm, and grey arrows correspond to the normalized velocities associated with these streamlines, in the center of the toroids.

Figure 16 .

 16 Figure 16. Examples of streamlines going through (x cur , z cur ) defined by the Z-model, with several values for the height of the flow-field center H 0 . Regular θ and R are shown, as well as θ and R , the equivalent values for a flow-field center not at the surface.

Figure 17 .

 17 Figure17. Normalized ejection speed as a function of normalized ejection distance, i.e., distance between impact point and position from where the particle is ejected, divided by the cylinder's radius, after a time of about 1.25 s, corresponding to the end of the cavity depth main increase. Our simulation data are compared toHousen et al. (1983) andWada et al. (2006) by respectively considering different final crater and transient crater radii.

Figure 18 .

 18 Figure 18. Ejecta volume as a function of time for four different impact speeds, with a regular timescale and an adapted one.

Figure 20 .

 20 Figure 20. Ejecta velocity angles from the horizontal as function of time, for a 300 m s -1 impact. Each ring represents an instant, and the quarter circle is divided into several 6 • angular sectors.

Figure 22 .

 22 Figure 21. Snapshots of simulations of the sampler and the regolith bed about 1 s and 4.5 s after the projectile has been fired. The colors correspond to the different tapered cylinders composing the sampling horn, where the smallest value corresponds to the lowest cylinder. The filter can be seen near the top of Cylinder 4.

Figure 23 .

 23 Figure 23. Number of particles in different cylinders of the sampler. Cylinder 1 is the lowest tapered cylinder and Cylinder 5 the highest one. Each color represents a simulation. Different ∆t/k n pairs were chosen to run these simulations, and are indicated in the legend.

Figure 24 .

 24 Figure24. Number of particles in different cylinders of the sampler. Cylinder 1 is the lowest tapered cylinder and Cylinder 5 the highest one. Each color represents a simulation. Different ∆t/k n pairs were chosen to run these simulations, and are indicated in the legend. In the legend is also indicated simulations with a second fired projectile after 0.2 s. The black dotted line represents the moment of second firing for the relevant simulations.

Figure 25 .

 25 Figure25. Normalized total ejecta volume of particles faster than a given speed as a function of normalized ejection speed, where R cyl is the radius of the cylinder. The data correspond to different ε n , 0.25 s after the impact.

Figure 27 .

 27 Figure 27. Number of particles in different cylinders of the sampler. Cylinder 1 is the lowest tapered cylinder and Cylinder 5 the highest one. Each color represents a simulation. Different normal coefficients of restitution were chosen to run these simulations, and are indicated in the legend.

Figure 29 .

 29 Figure29. Energy distribution (left) and impactor-particles overlaps (right) for a 300 m s -1 impact as functions of time, with an impact shift of x = x 0 -0.25 cm.

Figure 30 .

 30 Figure30. Cavity radius in cm as a function of time for various simulations, differing from each other in the impact geometry. All represent a 300 m s -1 impact. In the legend, degrees are from the vertical, in the x-z plane.

Figure 31 .

 31 Figure 31. Number of ejected particles whose ejection speed is larger than a given speed as a function of the ejection speed, for different impact angles. On this figure is only considered a speed range from about 10 m s -1 to 70 m s -1 .

Figure 32 .

 32 Figure32. Number of particles in different cylinders of the sampler. Cylinder 1 is the lowest tapered cylinder and Cylinder 5 the highest one. Each color represents a simulation. Different impact angles and locations were chosen to run these simulations, and are indicated in the legend. In the legend, degrees are from the vertical, in the x-z plane.

  Figure8. Different energies at stake for 50 m s -1 impacts with different timescales and spring constants. Energies are normalized with the impactor's initial kinetic energy and are represented for the first instants of the impact.

	Normalized energy	0.2 0.4 0.6 0.8 1.0	total Kproj Rproj Eel,norm,proj Eel,tang,proj Eel,roll,proj Eel,twist,proj Ktarget Rtarget Eel,norm,target Eel,tang,target Eel,roll,target Eel,twist,target	Normalized energy	0.2 0.4 0.6 0.8 1.0	total Kproj Rproj Eel,norm,proj Eel,tang,proj Eel,roll,proj Eel,twist,proj Ktarget Rtarget Eel,norm,target Eel,tang,target Eel,roll,target Eel,twist,target
		0.0			0.0
			0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 Time (s)		0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 Time (s)
						(b) ∆ t =	∆t 0
	Distance between impactor and other particle [cm]	0.0 0.2 0.4 0.6 0.8	Part n 8860 Part n 37046 Part n 62505 Part n 72841 Part n 8405
				0.0000 0.0001 0.0002 0.0003 0.0004 0.0005 0.0006 0.0007 Time [s]

(a) ∆ t = ∆ t 0 , k n = 0.58 k n 0 6 , k n = 36 k n 0

  Figure28. Energy distribution (left) and impactor-particles overlaps (right) for a 300 m s -1 impact as function of time, with an impact shift of x = x 0 + 0.25 cm.

	Normalized energy	0.0 0.2 0.4 0.6 0.8 1.0	total Kproj Rproj Eel,norm,proj Eel,tang,proj Eel,roll,proj Eel,twist,proj Ktarget Rtarget Eel,norm,target Eel,tang,target Eel,roll,target Eel,twist,target	Distance between impactor and other particle [cm]	0.0 0.2 0.4 0.6 0.8	Part n 62505 Part n 8860 Part n 37046 Part n 8405 Part n 72841
		0.00000 0.00002 0.00004 0.00006 0.00008 0.00010 0.00012 Time (s)			0.00000 0.00002 0.00004 0.00006 0.00008 0.00010 0.00012 Time [s]
	Normalized energy	0.0 0.2 0.4 0.6 0.8 1.0	total Kproj Rproj Eel,norm,proj Eel,tang,proj Eel,roll,proj Eel,twist,proj Ktarget Rtarget Eel,norm,target Eel,tang,target Eel,roll,target Eel,twist,target	Distance between impactor and other particle [cm]	0.0 0.2 0.4 0.6 0.8	Part n 8860 Part n 72841 Part n 12063 Part n 37046 Part n 8405
		0.00000 0.00002 0.00004 0.00006 0.00008 0.00010 0.00012 Time (s)			0.00000 0.00002 0.00004 0.00006 0.00008 0.00010 0.00012 Time [s]

© 2019 The Authors

MNRAS 000, 1-27 (2019)

This paper has been typeset from a T E X/L A T E X file prepared by the author.MNRAS 000, 1-27 (2019)

ACKNOWLEDGEMENTS F.T. is supported by a PhD fellowship of the University of Nice-Sophia. P.M. and F.T. acknowledge funding sup-port from the French space Agency CNES, as well as from Academies of Excellence: Complex systems and Space, environment, risk, and resilience, part of the IDEX JEDI of the Université Côte d'Azur in connection with its Center for Planetary Origins (C4PO). Simulations were performed at Mésocentre SIGAMM hosted at the Observatoire de la Côte d'Azur.