Henri Gouin
email: henri.gouin@univ-cezanne.fr

Laurent Espanet

Bubble Number in a Caviting Flow

Keywords: Gas/liquid flows, Bubble formation, Bubble dynamics, Cavitation PACS: 47.55.Ca, 47.55.db, 47.55.dd, 47.55.dp

Cavitation is a general phenomenon of the fluid flows with obstacles. It appears in the cooling conduits of the fast nuclear engines. A model of this phenomenon using the theory of Laplace and a common non-convex energy for the liquid and vapour bulks is proposed. This model makes it possible to determine a higher limit of the density of bubbles (a number of bubbles per unit of volume in the flow). The maximum intensity of cavitation is associated with the mechanical and thermal characteristics of the fluid flow.

Introduction

The Laplace theory applied to a closed tank allows us to obtain an equilibrium equation of a liquid in presence of vapour bubbles. We analytically prove that the total number of bubbles cannot exceed a maximum value. We extend the model to the mixture of fluid and gas in the case of steady flows without viscosity: the number of bubbles transported by the current admits an upper limit. It was observed in experiments that the number of bubbles created in a cavitating flow depends on the quality of the carrying fluid. It is well known that an injection of gas bubbles does not affect this upper limit. The number of particles and microscopic gas bubbles is a significant cause of the intensity of cavitation [START_REF] Briançon-Marjollet | La prévision de la cavitation en fonction de l'état de nucléation de l'eau et des conditions hydrodynamiques[END_REF][START_REF] Lecoffre | Aspects pratiques du contrôle des germes de cavitation en moyens d'essais[END_REF][START_REF] Gindroz | Lois de similitude dans les essais de cavitation des turbines Francis[END_REF]. Unfortunately, this parameter is generally inaccessible to measurement. The practical case is associated with the flows existing in coolant circuits of fast nuclear engines; the used cooling agent consists of a mixture of sodium and argon. Sodium in liquid form or vapour is called "fluid". Argon in form of gas, possibly dissolved, is a neutral constituent at a temperature much higher than its critical value. It is called "gas". It is not possible to analyze the quality of the mixture. To estimate a maximum number of bubbles produced by a possible phenomenon of cavitation is very significant. It is an essential data of the nuclear security.

Simplified model of cavitating flow

We consider an isothermal permanent two-phase flow made up of fluid (for example sodium at the temperature of 400 • C) and neutral gas far from its critical point (for example argon). We consider the approximation of a flow without viscosity. With an aim of simplifying the diagram of the system of cooling of fast nuclear engines, we separate the coolant circuit in two parts: -a first part of constant section s and volume v, in which the flow is carried along with cavitation. Cooling agent consists of sodium either in liquid form or in form of vapour bubbles; argon is dissolved in the liquid sodium or mixed with its vapour in the bubbles. The mixture contains n bubbles each one of volume v 1 . The bubbles are supposed in equilibrium in a reference frame convected by the flow and the edge effects of the pipe are not taken into account. The bubbles supposed without interactions are identical and in spherical form.

-a second part of constant section s 0 and volume v 0 without cavitation in which sodium is in liquid form and saturated by dissolved argon.

The motion of the mixture in each compartment is supposed to be uniform. One denotes respectively by u and u 0 the velocity values in each compartment. The loop is closed and contains a total mass M of fluid (sodium). A pump discharges the fluid according to an imposed flow d. We denote by ρ m the average density of the fluid in the cavitating part. The theorem of Bernoulli is supposed to be applicable to each of the two components of the fluid-gas mixture. This simplifying assumption associated with a barotropic permanent motion makes it possible to connect the velocity to the density of each component. The theorem of Bernoulli extends to dynamics the equality of the chemical potentials between phases and components [START_REF] Gouin | Variational theory of mixtures in continuum mechanics[END_REF]. These relations can be replaced by more advancing equations of the motion representing, for a permanent flow, the connections between velocities and densities.

The volumic free energy of the fluid-gas mixture is supposed to be the sum of partial volumic free energies. The pressure of mixture is the sum of the partial pressures. The set of the assumptions leads to the following relations:

ρ m v = ρ l v + nv 1 (ρ v -ρ l ) ,
which gives the value of the fluid mass in the zone of cavitation; ρ l and ρ v stand for the densities of liquid and vapour. It follows:

n = v v 1 ρ l -ρ m ρ l -ρ v . (1)
The total fluid mass in the coolant circuit is:

M = ρ m v + ρ 0 v 0 , (2)
where ρ 0 denotes the density of the fluid in the part without cavitation. A pump gives the fluid flow d in the circuit,

s ρ m u = s 0 ρ 0 u 0 = d. (3)
Let us denote by F the volumic free energy of the fluid. The theorem of Bernoulli applied to each phase of the fluid component gives:

1 2 u 2 + F ′ (ρ l ) = 1 2 u 2 0 + F ′ (ρ 0 ) (for the liquid phase), (4)
1 2 u 2 + F ′ (ρ v ) = 1 2 u 2 0 + F ′ (ρ 0 )
(for the vapour phase).

(

) 5
The approximations suppose a moderate cavitation. Although very synopsis, this schematization allows to get a first estimate of the maximum number of bubbles in the zone with cavitation. We take a state law of van der Waals type for expressing the fluid partial pressure. The free volumic energy F is given by a single expression in the liquid and the vapour phases [START_REF] Rocard | Thermodynamique[END_REF]. The relations (4) and ( 5) imply:

F ′ (ρ l ) = F ′ (ρ v ) (6)
which is the equality of the chemical potentials of the two phases of the fluid component. The free energy of the fluid being a non-convex function of ρ, the values of ρ l and ρ v are different. We denote by G the volumic free energy of the gas. We obtain G ′ (r v ) = G ′ (r l ) where r v is the density of argon in the bubble and r l the density of saturated gas dissolved in the liquid. Argon being a gas at a temperature far from its critical value, the function G is a convex function of the density r.

Consequently, r v = r l .
The partial pressures of the fluid and the gas are respectively:

p s (ρ) = ρF ′ (ρ) -F (ρ) et p a (r) = r G ′ (r) -G (r) .
Within the framework of the Laplace model, the equilibrium of a bubble in the zone of cavitation is:

p s (ρ v ) + p a (r v ) -p s (ρ l ) -p a (r l ) = 2σ R , (7)
where σ is the constant surface tension at the temperature of the flow and R is the radius of the bubble [START_REF] Bruhat | [END_REF]. Owing to the fact that r v = r l and with K denoting (36π) 1/3 σ, the relation ( 7) is written in the equivalent form:

F (ρ l ) -F (ρ v ) -F ′ (ρ l )(ρ l -ρ v ) = 2 3 Kv -1/3 1 . ( 8
)
This relation remains unchanged if adding to F an unspecified linear function of ρ. Thus, F ′ (ρ) is chosen null for the densities of saturation ρ ls and ρ vs of the liquid and the vapour.

The volumic free energy is supposed to be a regular function near the densities of saturation. The Taylor expansions to the second order in the vicinity of ρ ls and ρ vs yield:

F (ρ v ) = 1 2 (ρ v -ρ vs ) 2 F ′′ (ρ vs ), F (ρ l ) = 1 2 (ρ l -ρ ls ) 2 F ′′ (ρ ls ).
With taking into account the relation ( 6), we have

F ′′ (ρ vs )(ρ v -ρ vs ) = F ′′ (ρ ls )(ρ l -ρ ls ),
and the relation (8) yields:

F ′′ (ρ ls )(ρ ls -ρ l )(ρ l -ρ v ) 1 + 1 2 ρ ls -ρ l ρ l -ρ v 1 - F ′′ (ρ ls ) F ′′ (ρ vs ) = 2 3 Kv -1/3 1 . (9)
The value C l of the sound velocity in the liquid is

C 2 l = ∂p s ∂ρ l ρ ls = (ρF ′ (ρ) -F (ρ)) ′ ρ ls = ρ ls F ′′ (ρ ls ) .
The value C v of the sound velocity in the saturating vapour is

C 2 v = ρ vs F ′′ (ρ vs ) .
So, we obtain

F ′′ (ρ ls ) F ′′ (ρ vs ) = C 2 l C 2 v ρ vs ρ ls .
Far from the critical point of the fluid, we have

F ′′ (ρ ls ) F ′′ (ρ vs ) ≪ 1.
Such is the case of liquid sodium at the temperature of the nuclear engine [START_REF]Handbook of Chemistry and Physics[END_REF].

The range of values of ρ v and ρ l allows to write:

ρ ls -ρ l ρ l -ρ v ≪ 1.
So, from the respective values of the various terms of relation (9), one deduces the approximated relation:

v 1 = 2K 3C 2 l 3 1 (ρ ls -ρ l ) 3 .
(10)

An estimation of an upper bound value of the number of bubbles

With the approximation (10), the relation ( 1) is written in the form:

n = v v 1 (ρ l -ρ m ) ρ ls .
By taking into account relation (10), one obtains the number of bubbles in the flow:

n = 3C 6 l 32πσ 3 ρ ls v (ρ ls -ρ l ) 3 (ρ l -ρ m ) . (11)
Let us notice that ρ l and ρ m are not independent but bounded by the relations (2), ( 3) and (4) expressing ρ m as function of ρ l . Let us regard n as a continuous quantity; from these relations one deduces the differential system:

                   v dρ m + v 0 dρ 0 = 0, u dρ m + ρ m du = 0, u 0 dρ 0 + ρ 0 du 0 = 0, u du + F ′′ (ρ l )dρ l -u 0 du 0 -F ′′ (ρ 0 )dρ 0 = 0.
By eliminating du, du 0 et dρ 0 between the four relations one obtains:

u 2 ρ m + u 2 0 v ρ 0 v 0 -F ′′ (ρ 0 ) v v 0 dρ m = F ′′ (ρ l )dρ l .
The respective orders of magnitude of the different terms of the preceding relation and the fact that ρ 0 and ρ l are close to ρ ls yield:

-F ′′ (ρ 0 ) v v 0 dρ m = F ′′ (ρ l )dρ l .
One obtains the derivative of ρ m as a function of ρ l in the simplified form:

dρ m dρ l = - v 0 v .
Relation (11) where n is considered as a function only of ρ m , implies:

n ′ (ρ l ) = 3C 6 l 32πσ 3 ρ ls v -3 (ρ ls -ρ l ) 2 (ρ l -ρ m ) + (ρ ls -ρ l ) 3 1 + v 0 v .
For ρ l different of ρ ls , n ′ (ρ l ) is null for:

ρ l -ρ m = 1 3 (ρ ls -ρ l ) 1 + v 0 v . ( 12
)
Let us notice that n = 0 for ρ l = ρ ls and ρ l = ρ m . The value of ρ l associated with relation (12) corresponds to a maximum of bubbles in the flow. The value of the number of bubbles per unit of volume is N = n v such as:

N = C 6 l 32πσ 3 ρ ls (ρ ls -ρ l ) 4 1 + v 0 v . (13)
Instead of the density of liquid, it is more realistic to consider the flow pressure in the part of the circuit subjected to cavitation. Indeed, this physical quantity is more accessible to measurements.

Let us denote by P = p s (ρ l ) + p a (r v ) and P sat = p s (ρ ls ) + p a (r v ), where p s (ρ ls ) is the saturated vapour pressure and p a (r v ) is the partial pressure of the gas. One obtains P sat -P = C 2 l (ρ ls -ρ l ). One deduces the value of the maximum number of bubbles per unit of volume in the form:

N = (P sat -P) 4 32πσ 3 C 2 l ρ ls 1 + v 0 v . (14)
For example, we consider the physical values of liquid sodium at the temperature of 400 • C: surface tension σ = 0, 166 N/m ; sound velocity liquid sodium C l = 2350 m/s ; density of saturating sodium ρ ls = 856 Kg/m 3 [START_REF]Handbook of Chemistry and Physics[END_REF].

With choosing for the coefficient 1 + v 0 v the value 10 and for difference of pressures P sat -P = 4000 P a, one obtains in the part prone to cavitation a maximum of bubbles N = 1200 per liter .

Conclusion

We presented an approach allowing to consider a higher limit number of the density of bubbles present in a closed loop of cavitating flow. This number is given by the relation ( 14). The effects of uncondensable gases (argon in our case) do not appear in the relation owing to the fact that the fluid (sodium) is saturated with gas. The bubbles are supposed in local equilibrium, which requires that their times of creation and collapse are weak with respect to the running time in the cavitating part of the circuit (a fine study could be carried out starting from the equation of Rayleigh-Plesset [START_REF] Franc | La cavitation[END_REF]). The relation binding the average densities of the fluid ρ m and of the liquid ρ l partially results from the relation ( 4). This relation corresponds to a permanent flow. It may be replaced by all other relations binding u, u 0 , ρ l et ρ 0 . The method corresponds to a very simplified model of a flow of complex nature. Nevertheless, this model offers the advantage of showing that in a closed pipe, cavitation is a limited phenomenon. It should be noted that relation ( 14) takes into account the ratio of dimensions between the zones of cavitation and without cavitation. To decrease cavitation, one may reduce this ratio by creating additional zones where cavitation is of no importance on the maintenance of the coolant circuit.

Abridged French version

Résumé

La cavitation est un phénomène général des écoulements fluides avec obstacles. Elle apparait dans les conduites de refroidissement des réacteurs nucléaires rapides. Un modèle de ce phénomène utilisant la théorie de Laplace et une énergie non convexe commune aux phases liquide et vapeur est proposé. Il permet de déterminer une borne supérieure de la densité de bulles (nombre de bulles par unité de volume dans l'écoulement). L'intensité maximum de la cavitation est associée aux caractéristiques mécaniques et thermiques de l'écoulement fluide.

Version française abrégée

Il a été constaté expérimentalement que le nombre de bulles dans un écoulement cavitant dépend de la qualité du fluide. En particulier, le nombre de particules et de bulles microscopiques de gaz est un paramètre qu'il est important de conna ître [START_REF] Briançon-Marjollet | La prévision de la cavitation en fonction de l'état de nucléation de l'eau et des conditions hydrodynamiques[END_REF][START_REF] Lecoffre | Aspects pratiques du contrôle des germes de cavitation en moyens d'essais[END_REF]. Malheureusement, ce nombre est difficilement mesurable dans la plupart des écoulements. C'est notamment le cas pour l'écoulement de sodium liquide en réacteur nucléaire de type RNR. D'autres expériences ont révélé l'existence, pour un débit fixé, d'une valeur maximale pour le nombre de bulles; au delà de ce seuil, l'injection de microbulles de gaz ne semble plus modifier l'intensité de la cavitation [START_REF] Gindroz | Lois de similitude dans les essais de cavitation des turbines Francis[END_REF]. L'objectif de ce travail est de quantifier l'intensité de la cavitation à partir de l'estimation analytique du nombre de bulles dans un écoulement. Nous considérons un écoulement diphasique permanent isotherme constitué d'un "fluide" (par exemple du sodium à la température de 400 • C) et d'un "gaz" neutre très au-dessus de son point critique (par exemple de l'argon). Nous envisageons l'approximation d'un écoulement sans viscosité. Dans le but de simplifier le schéma du système de refroidissement des réacteurs nucléaires rapides, nous séparons le circuit de refroidissement en deux parties: -Une première partie de section constante s et de volume v, dans laquelle l'écoulement s'effectue avec cavitation. Le réfrigérant est constitué de sodium soit sous forme liquide soit sous forme de bulles de vapeur; l'argon est dissous dans le sodium liquide ou mélangé à sa vapeur dans les bulles. Le mélange contient n bulles, chacune de volume v 1 . Les bulles sont supposées en équilibre dans un repère entra îné par l'écoulement et les effets de bord de la conduite ne sont pas pris en compte. Les bulles supposées sans interactions entre elles sont alors identiques et d'une forme sphérique. La zone de cavitation correspond à la région d'échange de chaleur entre les barres radioactives et le fluide réfrigérant.

-Une deuxième partie supposée sans cavitation, de section constante s 0 et de volume v 0 dans laquelle le sodium est sous la forme liquide et saturé par de l'argon dissous. Le mouvement du mélange dans chaque compartiment est pris uniforme. On note respectivement par u et u 0 les valeurs de la vitesse dans chaque compartiment. Le circuit est fermé et contient une masse totale donnée M de fluide. Une pompe entra îne le fluide suivant un débit imposé d. Nous noterons ρ m la masse volumique moyenne du fluide dans la partie cavitante. Le théorème de Bernoulli est supposé applicable à chacun des deux constituants du mélange fluide-gaz. Cette hypothèse simplificatrice associée à un mouvement permanent et barotrope permet de relier la vitesse à la masse volumique de chaque constituant. Ces relations peuvent être remplacées par des équations du mouvement plus évoluées représentant, lorsque l'écoulement est permanent, les liaisons entre la vitesse et la masse volumique de chaque constituant. Considérons pour simplifier des bulles en équilibre dans un repère lié à cet écoulement. En utilisant la théorie de Laplace, l'énergie libre totale du milieu s'écrit comme la somme des énergies de chaque phase et de l'énergie superficielle des bulles. On obtient alors l'égalité des potentiels chimiques pour le constituant fluide et pour le gaz ainsi que l'expression du saut de pression à l'interface. Ces relations déterminent les caractéristiques géométriques des bulles de vapeur. On représente l'énergie libre du fluide supposée non convexe par deux paraboles dont les rayons de courbure aux sommets sont reliés à la vitesse du son dans chaque phase. L'ensemble des hypothèses effectuées conduit à l'expression (1) du nombre de bulles dans la zone de cavitation en fonction des densités du liquide et vapeur et de leur dimension, à la conservation de la masse totale de fluide (relation ( 2)), à l'expression du débit imposé (relation ( 3)). Le théorème de Bernoulli appliqué à chacune des phases du fluide conduit aux relations (4) et (5) d'où on déduit l'égalité (6) des potentiels chimiques [START_REF] Gouin | Variational theory of mixtures in continuum mechanics[END_REF]. De la définition de la pression à partir de l'énergie libre et de la loi de Laplace on obtient l'expression (8) reliant le rayon des bulles aux densités du liquide et de sa vapeur. En se limitant aux termes principaux, on montre que les deux relations ( 6) et [START_REF] Franc | La cavitation[END_REF] se réduisent alors à une seule expression: l'équation (10) reliant le volume des bulles à la densité du liquide dans la zone cavitante ρ l . La présence de gaz qui sature le fluide ne modifie pas la forme de l'équation de Laplace. Du fait des relations (2), (3) et (4), le nombre de bulles n peut être considéré comme un fonction uniquement de ρ l . On constate alors que le nombre n de bulles admet une valeur maximale. Cette valeur prend une forme explicitable en fonction de l'écart entre la densité du liquide et sa densité à la pression de saturation. C'est l'expression (13) donnant le nombre maximum de bulles par unité de volume. Elle constitue donc une première majoration analytique du nombre de bulles dans un écoulement cavitant. La densité du sodium liquide n'étant pas une grandeur physiquement mesurable, il convient de transformer la relation (13). En utilisant une loi d'état linéaire, on exprime le nombre maximum de bulles non plus en fonction de la densité du sodium liquide dans la zone de cavitation, mais en fonction de la pression totale. La relation obtenue (14) montre que le nombre maximum de bulles par unité de volume est fonction de l'écart entre la pression totale et la pression totale à la saturation. Un calcul numérique est effectué pour du sodium à la température des circuits de refroidissement des réacteurs nucléaires. Il constitue un ordre de grandeur réaliste du nombre de bulles dès lors que l'on est capable de mesurer la pression dans la zone cavitante.


Harrison F Agrusa
email: hagrusa@astro.umd.edu

Derek C Richardson

Alex B Davis

Eugene Fahnestock

Masatoshi Hirabayashi

Nancy L Chabot

Andrew F Cheng

Andrew S Rivkin

Patrick Michel

A Benchmarking and Sensitivity Study of the Full Two-Body Gravitational Dynamics of the DART Mission Target, Binary Asteroid 65803 Didymos

Keywords: Asteroids, dynamics, Rotational dynamics, Celestial mechanics, Near-Earth objects

NASA's Double Asteroid Redirection Test (DART) is designed to be the first demonstration of a kinetic impactor for planetary defense against a small-body impact hazard. The target is the smaller component of the binary asteroid 65803 Didymos. We have conducted high-fidelity rigid full two-body simulations of the mutual dynamics of this system in a broad benchmarking exercise to find the best simulation methodologies, and to understand the sensitivity of the system to initial conditions. Due to the non-spherical shapes of the components and their close proximity, the components cannot be treated as point masses and so the dynamics differ significantly from a simple Keplerian orbit, necessitating the use of numerical simulations to fully capture the system's dynamics. We find that the orbit phase (angular position or true anomaly) of the secondary is highly sensitive to the initial rotation phase of the primary, making prediction of the secondary's location from numerical simulation challenging. Finally, we show that the DART impact should induce significant free and forced librations on the secondary. If this libration can be measured by ESA's recently approved follow-up spacecraft, Hera, it may be possible to constrain properties of the secondary's interior structure.

Introduction

The Asteroid Impact & Deflection Assessment (AIDA) collaboration is a NASA-and ESAsupported effort to test the capability of a kinetic impactor for hazardous asteroid mitigation.

NASA will lead the Double Asteroid Redirection Test (DART) mission, which will achieve a kinetic impact on the secondary (Didymos B) of 65803 Didymos, a near-Earth binary asteroid, in the fall of 2022 [START_REF] Cheng | AIDA DART asteroid deflection test: Planetary defense and science objectives[END_REF]. ESA will lead Hera, a follow-up mission to rendezvous with Didymos to characterize the system and visible effects of the impact [START_REF] Michel | European component of the AIDA mission to a binary asteroid: Characterization and interpretation of the impact of the DART mission[END_REF].

The main goal of the DART mission is to demonstrate the kinetic impactor technique by intercepting the secondary, causing a change in the binary orbital period that can be measured with ground-based observations.

In this work, we present results from a suite of benchmarking simulations conducted by the DART Investigation Team's Dynamics Working Group to better understand the complex mutual dynamics, to constrain the sensitivity of the simulated Didymos system to initial conditions, and to identify the appropriate numerical methods to fully capture the dynamics. In Section 1.1 we describe the physical and dynamical properties of the Didymos system. Section 2 introduces the four dynamics models used in this study and the initial conditions of the simulation cases. The results are presented in Section 3. Finally, we discuss future work to better understand the Didymos binary in Section 4.

The Didymos System

Discovered in 1996, 65803 Didymos is classified as a near-Earth object and potentially hazardous asteroid with a heliocentric semi-major axis of 1.644 au [START_REF] Alday | Minor Planet Electronic Circulars[END_REF]. In 2003, it was discovered that Didymos is a binary system (Pravec et al., 2003). The binary has a mutual orbit period of P orb ∼ 11.9217 h and a semi-major axis of a orb ∼ 1.19 km [START_REF] Pravec | Photometric survey of binary near-Earth asteroids[END_REF]Naidu et al., 2020a). According to the binary mean separation and orbital period, Kepler's third law for assumed point masses implies a system mass of M sys ∼ 5.37 × 10 11 kg.

The primary (Didymos A) is ∼780 m across and has an oblate shape and equatorial ridge, and the secondary (Didymos B) is approximately 164 m across and we assume an ellipsoidal shape similar to that of other NEO binary secondaries. Didymos A is a fast rotator, with a spin period of 2.26 h. A polyhedral shape model with 1996 facets was derived by Naidu et al. (2020a) from combined radar and light curve data. In the simulations presented here, Didymos B is assumed to be synchronous (i.e., tidally locked) with its long axis initially aligned with the line of centers. It is also assumed that both bodies are in principal axis rotation and that their spin poles are initially aligned with the binary orbit normal.

Table 1 summarizes the relevant physical and dynamical parameters of the Didymos system. These are the nominal system parameters adopted by the DART investigation team at the current time and will be updated throughout the DART mission as new measurements become available.

Methodology

The Didymos binary is an example of the full two-body problem (F2BP), where the rotational and translational dynamics are fully coupled, due to the objects' irregular shapes and the close proximity of the components. As a result, the system's dynamical evolution is especially sensitive to the shapes and initial positions and orientations of each component, thus F2BP simulation codes are necessary to fully capture the system's dynamics.

The Simulation Codes

Four different codes were tested in this study, each developed by team members at NASA JPL, University of Colorado Boulder (UCB), Auburn University, and the University of Maryland (UMd), respectively. Only some of these codes have official names, so we refer to each code by the institution that developed it for simplicity. Brief descriptions of the codes are given below.

NASA JPL. The JPL code is based on the formulation of the mutual gravitational potential between two polyhedral bodies developed by [START_REF] Werner | Mutual Potential of Homogeneous Polyhedra[END_REF]. It calculates the mutual gravitational potential and its gradients through a Legendre polynomial series expansion, truncated to a desired order, and integrates the discrete-time Hamiltonian equations of motion using the Lie-Group Variational Integrator (LGVI) developed by [START_REF] Lee | Lie group variational integrators for the full body problem in orbital mechanics[END_REF]. This code was written in C++ Assumed.

Table 1: Physical and dynamical parameters of the Didymos System. These are the current nominal values adopted by the DART investigation team. Because these parameters are constantly being refined by ongoing observations, these are not exactly the same parameters used in this study. The initial conditions of the simulations presented here differ slightly, but remain within the uncertainty bounds given here. (See Table 3 for the simulation initial conditions.)

and parallelized to run on a cluster computer environment, due to the high computational cost of the potential and gradients evaluation at each timestep.

University of Colorado Boulder (UCB). The recently developed UCB code utilizes inertia integrals to expand the mutual gravitational potential according to the formalism derived by [START_REF] Hou | Mutual potential between two rigid bodies with arbitrary shapes and mass distributions[END_REF]. This tool, known as the General Use Binary Asteroid Simulator (GUBAS), is now publicly available1 and can easily be run on a single desktop computer. Despite the different mathematical formulations for the mutual gravitational potential and its gradients between the JPL and UCB codes, they agree to near-machine precision for the same given expansion order of the mutual potential, since they used the same numerical integrator (LGVI) for the simulations run herein.

However, the inertia integral formulation allows for the attitude and mass distribution to be decoupled and computed separately, which allows for a more computationally efficient implementation and thus faster runtimes. The present study served as a convenient test to confirm that the UCB code does in fact achieve the same result as the JPL code.

Auburn University. The Auburn code is a simplified version of the UCB code. It expands the inertia integrals only to second order according to the formulation given by [START_REF] Hirabayashi | Recursive computation of mutual potential between two polyhedra[END_REF]. The equations of motion are solved with an 8th-order Runge-Kutta scheme. This code can be thought of as evaluating the mutual gravitational potential of the system as if the Didymos A shape model were replaced with a best-fit ellipsoid. Although the Auburn code does not fully capture perturbations due to the asymmetric shape of the primary, it is extremely fast and is a useful reference point to understand the effect of higher-order perturbations due to Didymos A's shape.

University of Maryland (UMd). Unlike the other three codes that represent the primary and secondary as monolithic and homogeneous bodies of some arbitrary polyhedral or ellipsoidal shape, the UMd code treats each body as a rigid aggregate of many spherical particles. The code, called pkdgrav, is a parallel N-body tree code [START_REF] Richardson | Direct Large-Scale N-Body Simulations of Planetesimal Dynamics[END_REF][START_REF] Stadel | Cosmological N-body simulations and their analysis[END_REF]. The UMd code uses a primary consisting of ∼3500 particles in order for the average particle diameter (∼42 m) to be within the spatial resolution of the radar shape model (∼50 m). Details of UMd's representation of each body are shown in Table 2. The translational motion is integrated with a fixed-step second-order leapfrog integrator, while the rotational motion is integrated with a time-adaptive fourth-order Runge-Kutta scheme within each leapfrog step [START_REF] Richardson | Numerical simulations of asteroids modelled as gravitational aggregates with cohesion[END_REF]. Note that pkdgrav's k-d tree is not used, so the forces and torques are computed by summing over every particle at every timestep to ensure the highest possible accuracy at the expense of speed.

Snapshots of the simulations of Didymos are shown in Fig. 1. The key differences are that JPL, UCB, and Auburn simulate the full radar-derived Didymos A shape model with an ellipsoidal Didymos B, with the mutual gravitational potential expanded to some desired order of accuracy.

The UMd method fills a volume with randomly packed spherical particles, then carves each body to match the desired shape and computes the mutual gravitational potential of the packed spheres exactly without truncation to any order. Because of the different mass representations of each method, there will be inherent variations between these approaches.

Initial Conditions

The 11 simulations presented here comprise a small subset of cases that the Dynamics Working Group has studied thus far. These simulations were selected to compare code performance and better understand the system's sensitivity to uncertainty in its initial state. Using the nominal values for the mass of each body and their separation, we computed the initial conditions to approximately put the system on a circular Keplerian orbit. This is considered the "nominal" simulation case. More details on the initial conditions for the nominal case are shown in Table 3. We then give these initial conditions slight perturbations to test the system's sensitivity to the initial relative velocity of the secondary and the initial rotation phase of the primary. See Table 4 for details of the 11 test cases, along with a schematic in Fig. 2. Each group selected a timestep for their respective code, based on numerical convergence and runtime constraints, with each code conserving energy to one part in a million or better over the entire simulation. Each of the 11 test cases was simulated for a total of 150 days of simulation time. All codes modeled the system as two rigid bodies interacting purely through their mutual gravity, with all additional forces or torques such as solar tides, BYORP, or internal dissipation turned off.

Each group output previously agreed-upon state variables at a 1-minute cadence (except for UMd which had hourly outputs due to data storage constraints). With such a high output cadence, we were able to compare both short-and long-term evolution of the binary system with each simulation code. The JPL and UCB codes were run with the mutual gravity expansion set to 4th order; this choice is accurate enough to capture the dynamics with high fidelity while keeping the computation time manageable. UCB also repeated the nominal case with the gravity expansion order set to 8th order to confirm that the choice of 4th order was indeed sufficient to accurately model the system. Again, the Auburn code is limited to 2nd order, while the UMd code has no order truncation.

Results

Code Performance

Since each code was run on different machines (see acknowledgments) with different timesteps and numerical routines, normalized performance comparisons can be troublesome. In Table 5, we simply show the runtimes for the nominal case along with the number of processors and timesteps used by each code. The Auburn code is orders of magnitude faster than the other codes, given its 2nd-order approximation of the mutual potential. This makes it a useful tool for quick tests, however it does not capture higher-order perturbations due to the asymmetric shape of the primary.

It should also be noted that the UCB or JPL codes would have similar performance if the mutual potential approximation were set to 2nd order. The UMd code had the longest runtime, due to a combination of its small timestep and requirement to compute the gravitational potential on Equivalent to a spin period of 11.9216 h, in order to match the Keplerian orbit period. Aligned with orbit pole.

Table 3: Initial conditions of the nominal simulation.

Name Description nominal

The nominal, unperturbed initial state. The initial conditions are calculated using Newton's version of Kepler's Third Law to give a circular orbit based on the total system mass and mean separation.

posR, negR

The initial velocity of the secondary's barycenter is perturbed by ±0.0005 m/s in the instantaneous orbital radial direction.

posT, negT

The initial velocity of the secondary's barycenter is perturbed by ±0.0005 m/s in the instantaneous orbital tangential direction (along-track).

posN, negN

The initial velocity of the secondary's barycenter is perturbed by ±0.0005 m/s in the instantaneous orbital normal direction (out-of-plane).

ph+1, ph+3, ph-1, ph-3 the initial rotation phase of primary shape model is adjusted by rotating ±1 or ±3 degrees from nominal, around the primary spin pole. a particle-by-particle basis. The 4th-order UCB code offers the best combination of speed and accuracy, as we will see in the following section. Therefore, the Dynamics Working Group has recommended that the UCB code be adopted for future rigid-body dynamics studies related to DART.

The Nominal Case

The orbit period, semi-major axis, and eccentricity for each code's nominal case is shown in Table 6 along with what those values would be if the system were Keplerian. As expected, the NASA JPL and UCB (4th Order) results are nearly identical. They also match closely to the 8th-order result, indicating that the 4th-order approximation is capturing the mutual gravitational potential with high fidelity. The deviations in the orbit period and semi-major axis are driven by each code's representation of the mass distribution, and thus the mutual potential, of the two bodies.

Figure 3 shows the evolution of the system over 2 d (∼4 orbital periods) as determined by each simulation code for the nominal case. Since the bodies do not follow a precise Keplerian orbit due to their irregular shapes, the orbital eccentricity and inclination are osculating-they are instantaneous values evaluated based on the position and motion of the body centers for each simulation output.

Since each group used a numerically converged timestep, the differences they show in the system's orbital evolution are attributable to how each code represents the mass distribution of each body. The oscillations in the various orbital elements are small and are driven by the shape perturbations of the primary. The bottom two plots in Fig. 3 show the obliquities of each body, defined as the angle between the body's spin axis and the mutual orbit pole. Due to its 2nd-order gravity approximation, the Auburn code shows negligible out-of-plane motion, indicating that the small changes in the inclination and obliquities with the other three codes are mainly driven by asymmetries in the primary shape.

The NASA JPL and UCB codes at 4th order are indistinguishable in Fig. 3 and the 8th-order version shows almost no appreciable difference. All codes show qualitative agreement, given their known differences. Because the UMd code had an output frequency of one hour, the inclination and obliquity plots look artificially jagged. The NASA JPL and UCB codes show perfect agreement at 4th order, and the UCB 8th-order version matches closely as well. The Auburn code is evaluating the mutual gravitational potential to 2nd order, so it doesn't capture higher-order effects of the primary's asymmetric shape on the mutual orbit. The choppy noise in UMd's plot of Didymos B's Obliquity is a result of its coarser output cadence.

In order to distinguish each curve, this figure is best viewed in the online color version. To study the orbital elements over longer time scales, we take a running average to remove short-term behavior. This is shown for the nominal case over the full 5-month simulation in Fig. 4. The long-term evolution is quite stable, with the orbital elements remaining constant or drifting by a very small amount. One apparent effect is that each code has a different average semimajor axis and eccentricity. This is again a result of different mass representations, which sets the initial mutual potential, and thus the eccentricity and equilibrium separation. This highlights the influence of the primary's shape on the orbital properties of the system. In general, the codes agree well, with differences attributable to their respective representations of each body. Because we expect the Didymos binary to be stable over long periods, our ability to capture long-term stability in its mutual orbit is reassuring.

Primary Rotation Phase

The only constraints we place on the orientations of the primary and secondary are that their spin axes be initially aligned with the mutual orbit pole. In our nominal case, the primary's long axis is aligned with the line of centers at t = 0, but this choice is arbitrary. A precise measurement of the primary's orientation relative to the secondary at a given epoch with ground-based observations prior to the DART impact will be very challenging, so we treat the primary's initial rotation phase as a free parameter. Therefore, understanding the system's sensitivity to the initial primary rotation phase is essential to developing methods for predicting the position of the secondary at later times.

To test this sensitivity, we varied the initial primary rotation phase with respect to the nominal case by ±3 • and ±1 • . All codes showed a non-negligible dependence on this slight change. Fig. 5 shows the orbital phase (angular position of secondary) relative to each code's respective nominal case. The Auburn code is comparatively insensitive to the initial primary phase since it is only approximating the mass distribution to second order. NASA JPL and UCB have identical results, and UMd is slightly more sensitive to the initial primary phase.

If the initial primary rotation phase is altered, the initial mass distribution will be slightly different, resulting in a different mutual potential and thus a different orbital period. Due to the asymmetry of the primary and the binary's small separation, this is a non-negligible effect, especially if we want to accurately predict the position of the secondary.

After determining that the orbit may be sensitive to the initial primary rotation phase, we performed another set of simulations over a wider range of initial rotation phases with the UMd code only. The results in Fig. 6 show that the initial primary rotation phase has a significant influence on the orbital evolution of the secondary. After an integration time of only 30 d, there is a spread of ∼15 • in the relative positions of the secondaries. The DART Investigation Team has a requirement to predict the orbit phase at the impact epoch to within ±45 • 60 days prior to launch and to within ±15 • 55 days prior to impact (3σ errors). The high sensitivity to Didymos A's initial rotation phase means that it will likely be impossible to meet this requirement with dynamical simulations, especially considering the uncertainties in the other initial conditions. However, the Observing Working Group should be able to meet this orbital phase prediction requirement through This represents the spread in the possible locations of Didymos B after a fixed interval of time given some random initial primary rotation phase. In order to distinguish each curve, this figure is best viewed in the online color version.

fitting a weighted least-squares model to observed timing of mutual events (Naidu et al., 2020b).

Didymos B Libration

Through tidal dissipation, we expect that the mutual orbit has circularized, the secondary is tidally locked, and any libration of the secondary's spin state has damped to a minimum. So the Didymos system should be in or close to a dynamically relaxed state prior to the DART impact.

The impact will nearly instantaneously reduce the instantaneous orbital velocity of the second-ary, decreasing the orbit period and increasing the eccentricity. Therefore, significant libration of the secondary should be induced. A libration angle is a measure of the orientation of a satellite's long axis relative to the line of centers between the two components' centers of mass. In the following analysis, we show only results from the UCB 4th order code for brevity, although we note that all 4 codes show good agreement, given the known differences among the codes.

In the coupled spin-orbit problem in which a synchronous, ellipsoidal secondary orbits a spherical or point-mass primary with its spin axis aligned with the mutual orbit pole, there are two modes of libration: excited and relaxed. In the decoupled spin-orbit problem, the excited and relaxed modes are analogous to free and forced librations, respectively. See [START_REF] Naidu | Near-earth asteroid satellite spins under spin-orbit coupling[END_REF] for a detailed discussion on these two libration modes in both the coupled and uncoupled scenarios.

The frequency of free libration for a synchronous satellite on a circular orbit is given by,

ω 0 = n 3 B -A C 1/2 , ( 1
)
where n is the mean motion, and A, B, and C are the three principal moments of inertia of the secondary, where A < B < C [START_REF] Murray | Solar System Dynamics[END_REF]. This mode is analogous to a pendulum's natural frequency, depending on its length and the gravitational acceleration. This libration mode is thought to be damped away due to tidal friction, especially if the secondary has a rubble-pile structure [START_REF] Murray | Solar System Dynamics[END_REF][START_REF] Goldreich | Tidal evolution of rubble piles[END_REF]. However, the forced (relaxed) mode necessarily exists for a synchronous secondary on an eccentric orbit. The secondary will feel a periodic restoring torque, with a frequency equal to the mean motion, due to the misalignment of the long axis with the line of centers, resulting from the orbital angular velocity varying over the course of a single orbit.

This picture is complicated when we consider libration in the full two-body problem. The theory on spin-orbit coupling discussed above makes two critical assumptions: 1) that the orbit is fixed (no apsis precession) and 2) that the ellipsoidal secondary is orbiting a spherically symmetric primary [START_REF] Wisdom | Urey prize lecture: Chaotic dynamics in the solar system[END_REF]. As a result, we will see some differences between the classic theory and our simulation results.

The nominal DART spacecraft impact is designed to hit the secondary's center-of-figure, in a direction nearly opposite its orbital motion at a 15-to-25-degree angle with respect to the orbital plane (depending on DART's trajectory), imparting a near-instantaneous change to its orbital velocity without significantly altering its spin state [START_REF] Cheng | AIDA DART asteroid deflection test: Planetary defense and science objectives[END_REF]. This will induce both longitudinal (in-plane) and latitudinal (out-of-plane) librations that will have both relaxed and excited components. By changing the mean motion without a matching change in the secondary's spin rate, we introduce excited (free) libration modes on top of those that exist already (if any). Further, the impact will increase the eccentricity of the system, exciting a stronger relaxed (forced) libration mode. Therefore, studying the resulting libration for the benchmarking cases where we perturb the orbital motion of the secondary reveals the extent to which DART may effect a libration in the secondary. Further, understanding the behavior of induced librations may be an important tool for interpreting the results of the DART impact, if the libration amplitude or frequency is observable.

The longitudinal libration for the nominal case is shown in Fig. 7a for the first 10 days of the simulation. The distinct beating pattern in the libration is a signature of both excited-mode and relaxed-mode librations [START_REF] Naidu | Near-earth asteroid satellite spins under spin-orbit coupling[END_REF]. A Fourier transform of the libration pattern shows the two distinct libration modes (Fig. 7b).

The frequency of the relaxed (forced) mode is the frequency at which the orbital angular velocity oscillates, which in this case is the epicyclic or radial frequency. A key assumption in the classic spin-orbit problem is that the orbit is fixed (i.e., no apsidal precession), in which case the epicyclic frequency would match the mean motion. However, the oblate shape of the primary and the close orbit of the secondary results in an extremely fast precession of the periapse. In the nominal case, the mean motion differs from the epicyclic frequency by ∼1%, which corresponds to a precession rate of ∼3.5 • per orbit.

The excited (free) libration mode has a frequency close to the theoretical prediction given by Eq. 1. These frequencies don't match perfectly because Eq. 1 assumes a spherically symmetric primary on a fixed orbit. Because this excited libration frequency will depend on the secondary's moments of inertia in a fashion similar to Eq. 1, it may be possible to infer something about the mass distribution and interior structure from a careful measurement of the libration frequency with Hera.

Fig. 8 shows the libration for the nominal case and the 2 cases where the secondary was given an along-track velocity perturbation (posT/negT). The libration amplitude is driven by the initial difference between the orbital angular velocity and the secondary's spin rate. The posT case is where the secondary is given a slightly larger initial tangential velocity, so its orbit expands, increasing the orbital period to ∼11.911 h, closely matching Didymos B's initial spin period and mode is driven by the epicyclic (radial) frequency, while the excited (free) mode is controlled by the moments of inertia of the secondary. The theoretical free libration frequency doesn't match the excited mode perfectly because its derivation assumes a spherically symmetric primary.

thus decreasing the libration amplitude. The negT case is the opposite: a smaller initial tangential velocity shrinks the orbit and shortens the orbital period (∼11.703 h), producing a larger discrepancy between the secondary's initial spin and orbital periods, thus a libration amplitude reaching ∼8 • at its maximum.

The perturbation to the secondary's linear momentum in the posT and negT cases is approximately one-half of the momentum carried by the DART spacecraft, so these perturbations are of the same order of magnitude that DART may produce. Since the posT case has a relatively small libration amplitude, we can think of it as being close to the "true" relaxed state of the system (in which the excited libration mode has nearly damped away but relaxed librations persist). Then, the jump from the posT to the negT cases corresponds to a rough conservative estimate of the effect of the DART impact on the libration, when the momentum perturbation to the secondary is approximately equal to the momentum carried by the DART spacecraft. In reality, we would expect the momentum perturbation to the secondary to be considerably larger, due to the contributions of ejecta to the net momentum transfer. Further, if the DART spacecraft impacts several meters off of center-of-figure, which terminal guidance simulations at JHU/APL suggest is likely, the torque applied to the secondary will also alter its spin state, nearly instantaneously. The DART terminal guidance system will likely result in an impact location biased toward the illuminated portion of the secondary, which will be the side opposite the primary based on DART's viewing geometry at the impact epoch. Therefore, such an off-center impact is likely to reduce the secondary's angular velocity, further increasing the maximum possible libration amplitude. For these reasons, the simulations presented here are a conservative estimate of the possible post-impact libration state of Didymos B, given our current knowledge of the state of the system. It is also important to note, that the libration amplitude and frequency are dependent on Didymos B's moments of inertia which are computed based on our assumptions of constant density and its ellipsoidal shape.

Studying the dependence of Didymos B's libration on its mass distribution is planned for a future study. [START_REF] Naidu | Near-earth asteroid satellite spins under spin-orbit coupling[END_REF] show that the libration of a synchronous satellite may be detectable with radar, if the secondary is large enough compared to the primary (D S /D P 0.2). This threshold is barely satisfied by the Didymos system and therefore the libration may be measurable with The longitudinal libration is insensitive to the normal perturbation as expected, while the small induced latitudinal libration is caused by initial out-of-plane motion of the secondary. In order to distinguish each curve, this figure is best viewed in the online color version.

Conclusions

In this work, we found that the simulation package provided by UC Boulder is well-suited to studying the orbital dynamics of the Didymos system, due to its accuracy and speed. The Dynamics Working Group has recommended the adoption of this code for future dynamics studies in support of DART.

The results of this benchmarking study show that: 1) shape perturbations cause a non-negligible deviation from a Keplerian orbit; 2) the orbit phase of the secondary is highly dependent on the initial orientation of the primary; and 3) the system will be highly susceptible to induced librations resulting from the DART impact, which may be measurable from ground-based radar or with Hera.

If measurable, Didymos B's libration may be a useful probe of its internal structure.

The first two results indicate that predicting the orbital phase of the secondary may not be feasible with numerical simulations, given the uncertainties in the initial conditions and body shapes. However, the Observing Working Group will be able to meet this orbital phase prediction requirement through fitting an analytic model to observed timing of mutual events.

We have begun a comprehensive study with the UCB code on the strength and frequency of post-impact librations as a function of the mass distribution of Didymos B and momentum transferred by DART. This study will be used to constrain the range of possible impact outcomes in order to better infer the result of the actual DART experiment.

In reality, Didymos A is likely a rubble pile given its shape and fast rotation. A rubble-pile structure may play an important role for the binary dynamics due to processes such as landslides [START_REF] Hirabayashi | Assessing possible mutual orbit period change by shape deformation of Didymos after a kinetic impact in the NASA-led Double Asteroid Redirection Test[END_REF] and tidal dissipation. Didymos B may also be a rubble pile in which case its free libration modes will dampen via internal friction. Therefore, we also plan to use pkdgrav with an implementation of a soft-sphere discrete element method to numerically investigate whether this damping may be a noticeable effect over the timescales between the DART and Hera missions.

  (a) Radar shape model of Didymos A and the assumed ellipsoidal shape of Didymos B. JPL, UCB, and Auburn simulate the full radar-derived primary shape model with the mutual potential expanded to various orders of accuracy.(b) The UMd representation of the binary, where randomly packed spherical particles fill the shapes of each body and the potential is computed explicitly over every particle.

Figure 1 :

 1 Figure 1: Representations of the Didymos binary among the different codes. Both of these images are a top-down view (i.e. from the mutual orbit north pole) at the start of the simulation.

  Figure 2: Schematic of the perturbations to Didymos B. The X, Y, and Z axes are the three principal axes of Didymos A. The nominal case has the long axes of both bodies aligned. The ph±3 and ph±1 cases have Didymos A rotated about its spin axis such that its X axis is pointed ±3 or ±1 degrees away from the direction to Didymos B.

Figure 3 :

 3 Figure3: The nominal case over 2 d (∼4 orbital periods). The NASA JPL and UCB codes show perfect agreement at 4th order, and the UCB 8th-order version matches closely as well. The Auburn code is evaluating the mutual gravitational potential to 2nd order, so it doesn't capture higher-order effects of the primary's asymmetric shape on the mutual orbit. The choppy noise in UMd's plot of Didymos B's Obliquity is a result of its coarser output cadence. In order to distinguish each curve, this figure is best viewed in the online color version.

Figure 4 :

 4 Figure 4: A running average of orbit parameters for the nominal case over 5 months. The first and last 10 d of data are cut off due to the running average using ∼ 10 days worth of data (20 orbital periods). In order to distinguish each curve, this figure is best viewed in the online color version.

  vs. Initial Primary Rotation PhaseAll quantitites are relative to their respective nominal run

Figure 5 :Figure 6 :

 56 Figure5: Orbit phase relative to each code's nominal run (the horizontal line). Small differences in the initial primary rotation phase give a slightly different orbital period, making it difficult to predict the position of the secondary at later times. Note the differences in vertical scale between codes.

  Longitudinal libration of UCB nominal case over 10 days. The libration pattern is consistent over the full 150-day simulation. transform of longitudinal libration.

Figure 7 :

 7 Figure7: The Fourier transform of the longitudinal libration reveals the two libration modes. The relaxed (forced) mode is driven by the epicyclic (radial) frequency, while the excited (free) mode is controlled by the moments of inertia of the secondary. The theoretical free libration frequency doesn't match the excited mode perfectly because its derivation assumes a spherically symmetric primary.

  Longitudinal libration for out-of-plane perturbations.

  Latitudinal libration for out-of-plane perturbations.

Figure 9 :

 9 Figure9: Longitudinal and latitudinal libration vs. time for normal (out-of-plane) perturbations, using the UCB code. The longitudinal libration is insensitive to the normal perturbation as expected, while the small induced latitudinal libration is caused by initial out-of-plane motion of the secondary. In order to distinguish each curve, this figure is best viewed in the online color version.

  Derived based on D P and M sys .

	Symbol	Parameter	Value	Comments/References
	a orb	Semi-major Axis	1.19 ± 0.03 km	(Naidu et al., 2020a)
	c : b : a Secondary Axis Ratios	1:1.2:1.56	Assumed, based on other binary systems.
	D P	Diameter of Primary	780 ± 30 m	(Naidu et al., 2020a)
	D S	Diameter of Secondary	164 ± 18 m	Derived from D P and D S /D P
	D S /D P	Size Ratio	0.21 ± 0.01	(Scheirich and Pravec, 2009)
	e orb	Binary Orbit	e orb < 0.03	Upper limit, assumed zero. (Scheirich and
		Eccentricity		Pravec, 2009)
	i orb	Binary Orbit	0.0	Assumed.
		Inclination
	(λ, β) (310 P orb Mutual Orbit Pole Binary Orbit Period 11.9217 ± 0.0002 h	One possible orbit solution. (Scheirich and
				Pravec, 2009)
	P P	Primary Spin Period	2.2600 ± 0.0001 h	(Pravec et al., 2006)
	P S	Secondary Spin Period	11.9217 h	Assumed.
	ρ P	Primary Bulk Density	2170 ± 350 kg m -3
	ρ S	Secondary Bulk	2170 ± 350 kg m -3
		Density

• , -84 • ) ± 10 • Ecliptic coordinates, (Scheirich and Pravec, 2009; Naidu et al., 2020a) M sys Total System Mass (5.37 ± 0.44) × 10 11 kg Derived via Kepler's 3rd Law with P orb and a orb .

Table 2 :

 2 Physical parameters for pkdgrav particles. All particles in each respective body are given a uniform particle density in order to achieve the desired bulk density of 2.104 g cm -3 . r avg is the mean particle radius. The particle size distributions for the primary and secondary are sampled from a normal distribution with mean µ = 21 m, standard deviation σ = 4.2 m and µ = 4.4 m, σ = 0.88 m, respectively. The size distributions both have ±1σ cutoffs.

		Primary Secondary
	N	3355	3546
	ρ bulk [g cm -3 ]	2.104	2.104
	ρ particle [g cm -3 ] 3.8947	3.7329
	r avg [m]	21	4.4

Table 4 :

 4 Description of the 11 simulation cases.

Table 6 :

 6 Simulated time-averaged orbital period, semi-major axis, and eccentricity for the nominal case. All codes use exactly the same initial conditions and body masses, so the deviations from a Keplerian orbit and among the codes themselves are due to different mass representations of the primary and secondary. Each quantity is rounded to enough decimal points to show deviation between nearly identical numbers.

		Orbit Period [h] Semi-Major Axis [km]	Eccentricity
	Auburn	11.8138797	1.1747401	0.0055961390
	NASA JPL	11.8062721	1.1743513	0.0059765762
	UCB (4th Order)	11.8062824	1.1743520	0.0059765772
	UCB (8th Order)	11.8053794	1.1743062	0.0060210983
	UMd	11.8211246	1.1750370	0.0051450228
	Kepler Orbit	11.9216030	1.1800000	0.0000000000

https://github.com/alex-b-davis/gubas

Acknowledgments

This study was supported in part by the DART mission, NASA Contract # NNN06AA01C to JHU/APL. Some of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). H.F.A. would like to thank Doug Hamilton for useful discussions. A.B.D. was supported by the National Science Foundation Graduate Research Fellowship Program under Grant No. DGE 1650115. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. M.H. acknowledges the Auburn University intramural grant program support. P.M. acknowledges the European Space Agency (ESA) and the French Space Agency (CNES) for support. Some of the simulations herein were carried out on: 1) The University of Maryland Astronomy Department's YORP cluster, administered by the Center for Theory and Computation; 2) the Stam-pede2 cluster at the Texas Advanced Computing Center (TACC), University of Texas at Austin, funded by NSF award ACI-1134872; and 3) the nexus cluster of the Mission Design and Navigation Section at JPL. Raytracing for Fig. 1b was performed using the Persistence of Vision Raytracer (http://povray.org/).

radar, given adequate observing conditions, a favorably shaped secondary, and a sufficiently large momentum transfer.

When the secondary is given a normal (out-of-plane) perturbation, we are inducing out-ofplane motion in the secondary and therefore a noticeable, but small, latitudinal libration in the secondary, while the longitudinal libration is effectively unchanged (see Fig. 9). The librations resulting from the benchmarking cases where we apply radial velocity perturbations (posR/negR) and different primary rotation phases (ph±3 and ph±1) show almost no sensitivity, so we exclude showing them here.


M Arakawa

T Saiki

K Wada

K Ogawa

T Kadono

K Shirai

H Sawada

K Ishibashi

R Honda

N Sakatani

K Ogawa 21

Y Iijima

C Okamoto

H Yano

Y Takagi

M Hayakawa

P Michel

M Jutzi

Y Shimaki

S Kimura

Y Mimasu

T Toda

H Imamura

S Nakazawa

H Hayakawa

S Sugita

T Morota

S Kameda

E Tatsumi

Y Cho

K Yoshioka

Y Yokota

M Matsuoka

M Yamada

T Kouyama

C Honda

Y Tsuda

S Watanabe

M Yoshikawa

S Tanaka

F Terui

S Kikuchi

T Yamaguchi

N Ogawa

G Ono

K Yoshikawa

T Takahashi

Y Takei

A Fujii

H Takeuchi

Y Yamamoto

T Okada

C Hirose

S Hosoda

O Mori

T Shimada

S Soldini

R Tsukizaki

T Iwata

M Ozaki

M Abe

N Namiki

K Kitazato

S Tachibana

H Ikeda

N Hirata

R Noguchi

A Miura

An artificial impact on the asteroid (162173) Ryugu formed a

come   L'archive ouverte pluridisciplinaire

Abstract: The Hayabusa2 spacecraft investigated the small asteroid Ryugu, which has a rubble pile structure. We describe an impact experiment on Ryugu using Hayabusa2's Small Carry-on Impactor (SCI). The impact produced an artificial crater with a diameter >10 meters, which has a semicircular shape, an elevated rim and a central pit. Images of the impact and resulting ejecta were recorded by the Deployable CAMera 3 (DCAM3) for >8 minutes, showing the growth of an ejecta curtain and deposition of ejecta onto the surface. The ejecta curtain, outer edge of ejecta, was asymmetric, heterogeneous, and never fully detached from the surface. The artificial crater was formed in the gravity-dominated regime: the gravity limited the crater growth; we discuss the implication of this observation for Ryugu's surface age.

Main Text: The Hayabusa2 spacecraft reached the near-Earth asteroid (162173) Ryugu in June 2018, showing that this small asteroid is covered with boulders with sizes up to ~160 m [START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF][START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu-A spinning top-shaped rubble pile[END_REF]. The surface is covered with a regolith layer composed of boulders and granules under microgravity conditions (gravitational acceleration about 1×10 -4 m s -2 ), which was expected to have strength originating from cohesion forces between regolith grains, with the maximum theoretical strength of the surface layer predicted to be 1 kPa [START_REF] Wada | Asteroid Ryugu before the Hayabusa2 encounter[END_REF]. Impact craters formed under those conditions should be limited in size by this surface strength, which reduces the expected crater diameter compared to a strengthless surface [START_REF] Holsapple | Asteroid impacts: Laboratory experiments and scaling laws[END_REF]. Crater sizes produced by hypervelocity impacts on small bodies can be predicted using a variety of crater scaling laws. The number and sizes of craters on asteroids (including Ryugu) can be used to date their surfaces, but different crater scaling laws give ages that differ by more than an order of magnitude (1).

Hayabusa2 spacecraft was equipped with a Small Carry-on Impactor (SCI) designed to launch a 2kg copper projectile at a velocity of 2 km s -1 to form an artificial impact crater (referred to hereafter as the SCI crater) on the surface of Ryugu. The SCI crater was intended to expose the asteroid subsurface for remote sensing and sample collection [START_REF] Arakawa | Scientific objectives of Small Carry-on Impactor (SCI) and Deployable Camera 3 Digital (DCAM3-D): observation of an ejecta curtain and a crater formed on the surface of Ryugu by an artificial high-velocity impact[END_REF][START_REF] Tachibana | Hayabusa2: Scientific importance of samples returned from C-type near-Earth asteroid (162173) 1999 JU3[END_REF]. The Deep Impact spacecraft performed an impact experiment on a small body and formed an artificial impact crater on the nucleus of the comet 9P/Tempel 1 (7)[START_REF] Schultz | Expectations for crater size and photometric evolution from the Deep Impact collision[END_REF][START_REF] A'hearn | Deep impact: excavating comet Tempel 1[END_REF]. The SCI impact can be used to test crater scaling laws, for the crater size and ejecta velocity distribution, which were constructed by laboratory experiments [START_REF] Matsue | Measurements of seismic waves induced by high-velocity impacts: Implications for seismic shaking surrounding impact craters on asteroids[END_REF][START_REF] Tatsumi | Cratering efficiency on coarse-grain targets: Implications for the dynamical evolution of asteroid 25143 Itokawa[END_REF][START_REF] Cintala | Ejection-velocity distributions from impacts into coarsegrained sand[END_REF], and be compared with numerical simulations of the impact cratering processes in a microgravity environment [START_REF] Jutzi | Modeling asteroid collisions and impact processes[END_REF].

The SCI impact operation was carried out on 5 April 2019 (14). The impact and evolution of impact ejecta were simultaneously observed by Hayabusa2's Deployable CAMera 3 (DCAM3) (14). About 3 weeks after the SCI impact, Hayabusa2 searched for evidence of the impact crater using its Optical Navigation Camera-Telescopic (ONC-T) from an altitude of 1.7 km with a resolution of ~18 cm/pixel. The fresh crater was identified at latitude 7.9° N and longitude 301.3° E (Figure 1A-B), ~20 m from the aiming point (which was latitude 6.00° N, longitude 303.00° E). This location is in a northern part of Ryugu's equatorial ridge (2). The impact angle was estimated to be ~60° from the local horizontal surface (14).

Comparison of ONC-T images taken from the same altitude before and after the impact (Fig. 1A-B) shows the formation of an artificial impact crater. A dozen boulders with sizes of at least tens of cm were driven away from the surface, and a large block with a size of 5 m (which we refer to as the Mobile Block, MB) was excavated and moved ~3 m northwest and >1 m above the initial surface by the impact. Before the impact, MB was adjacent to another large block (Stable Block, SB) (Fig. 1A), but SB was hardly moved by the SCI impact. SB might be part of a larger block buried below the crater floor. We speculate that SB and its subsurface extension halted the crater growth toward the southeast. Alternatively, SCI may have impacted closer to MB than SB, allowing the higher pressure to move MB but not SB.

To investigate the crater morphology profile, we generated a Digital Elevation Map (DEM) (Fig. 2A) from ONC-T images (14). This shows an elevated rim with a height of ~40 cm surrounding the SCI crater (Fig. 2B). The elevated rim could be composed of ejecta deposits and/or structural uplifts. We used the DEM to produce an ortho-corrected (a geometrically accurate) image on which we traced the top of the rim to determine the crater shape (Fig. 1C).

The SCI crater is semi-circular with a rim to rim diameter 𝐷 rim of 17.6 ± 0.7 m. This shape is not due to the oblique impact but the effect of a large block (possibly SB) southeast of the impact point (Fig. 1B) (14). Impact craters are usually symmetrical at this impact angle, so we assume that the SCI impact point is the center of the semicircular crater rim. A crater's apparent diameter 𝐷 is conventionally measured at an initial surface elevation. Using the DEM, we determine the SCI crater to have D = 14.5 ± 0.8 m (Fig. 2A).

There is a small pit close to the impact point on the crater floor. The pit entrance is at a depth of 1.7 m from the initial surface (Fig. 2B); the diameter of the pit is ~3 m and its depth is 0.6 m. The pit has a conical shape similar to the simple craters formed in laboratory experiments [START_REF] Matsue | Measurements of seismic waves induced by high-velocity impacts: Implications for seismic shaking surrounding impact craters on asteroids[END_REF][START_REF] Melosh | Impact cratering: A geologic process[END_REF]. Central pits are commonly observed in craters smaller than a few hundred meters on the Moon because the lunar subsurface has a layered structure consisting of a relatively hard layer covered with softer cohesionless regolith [START_REF] Yue | Lunar regolith thickness deduced from concentric craters in the CE-5 landing area[END_REF]. Therefore, the pit could be explained if the SCI impact occurred above a subsurface layer with a cohesional strength. Some small natural craters on Ryugu also have pit features (fig. S4). The cohesional strength of this putative subsurface layer would be between 140 and 670 Pa (14). Alternatively, a central pit could be formed on compressive material if the projectile density is larger than that of the targets (7,[START_REF] Schultz | Expectations for crater size and photometric evolution from the Deep Impact collision[END_REF]. Because the pit depth is 2.7 m from the rim top, the depth to rim diameter ratio of the SCI crater is 0.15 ± 0.01. This value is consistent with the natural craters found on Ryugu, which range from 0.14 to 0.20 [START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF].

Several boulders with sizes of 0.6 to 1.9 m are present within the crater. The crater wall appears to be smooth (Fig. 1D) because the boulders are buried inside the wall. They could have been exposed by the excavation flow of subsurface materials during the crater growth. The sizefrequency distribution (fig. S5) of boulders embedded in the wall is about one third of that on a nearby area of undisturbed surface (the C01-C region at 10°N, 300°E, about 20 m northwest of the impact point), where the size-frequency distribution in the undisturbed area could be similar to the impact point before the SCI impact in submeter range. The SCI crater wall must be composed primarily of regolith with rock sizes smaller than 20 cm. Assuming that the crater wall represents the subsurface structure, we conclude that the subsurface layer is dominated by regolith with rock sizes smaller than 20 cm. DCAM3 observed the impact ejecta from the moment of impact until more than 8 minutes later. Figure 3A-F shows the images taken by DCAM3 from 5 s to 498 s after the impact. The image at 5 s shows the impact ejecta curtain growing northwards from the surface at an angle of 40° (Fig. 3A, fig. S6A). Subsequent images show the impact ejecta curtain was asymmetric and contained multiple rays in the early stages. Ejecta growth toward the south is not observed, but low speed ejecta at the base of the ejecta curtain are visible in faint scattered light.

At 36 s, the northward ejecta ray becomes clearer and more extensive. Two more ejecta rays appear at this time (Fig. 3B, fig. S6B). Multiple rays of ejecta were previously observed in the Deep Impact experiments (7). A continuous portion of ejecta curtain, which appears as a shadow at the base of the ejecta curtain, begins to grow vertically and horizontally. This growth could originate from the expanding crater due to ongoing excavation of the subsurface material.

These three ejecta rays continue to grow in both size and thickness, becoming more visible at 100 s and 192 s (Fig. 3C-D; fig. S6C-D), indicating that the crater growth is continuing at these times. An additional ejecta ray appears as a dark shadow in these images at the center of the ejected region, so we interpret it as moving toward the line of sight of DCAM3, but there is no visible ejecta curtain in the southeast area. The surface of Ryugu is covered with many large boulders with a power law size distribution (1), which we presume also applies to boulders buried in the subsurface. The SCI impact area was covered with boulders >0.6 m, some of which had most of their volume buried in the ground (Fig. 1A). These boulders could cause a heterogeneous ejecta growth, which would then result in ejecta rays, i.e. the boulders on or under the surface might stop the excavation flow, separating the ejecta curtain into several rays. The ejecta curtain formed by the Deep Impact collision exhibited an up-range zone of avoidance (ZoA) (7). We see a similar ZoA for the SCI impact, which may be caused by the oblique nature of the impact, as in Deep Impact. However, the resulting semicircular shape of the SCI crater demonstrates the influence of large boulders.

We analyzed the distribution of ejecta rays around the SCI crater using the image at 192 s (Fig. 4A) together with ONC-T images of the SCI crater area obtained before the impact. The direction and covered surface area of each ejecta ray were identified on the images taken by ONC-T with reference to the DCAM3 line of sight. Fig. 4B shows the distribution of the main ejecta rays mapped on an ONC-T surface image; each numbered sector on the map corresponds to the numbered separated ejecta ray in Fig. 4A, while the continuous ejecta curtain corresponds with a semicircle at the center. The whole ejecta, including the rays as well as the continuous curtain, was only distributed on the northwest-side. This distribution map can be compared with an ONC-T v-band (0.55 µm) reflectance difference map shown in Fig. 4C (14): we interpret regions with larger decrease of the reflectance factor as being covered with a thicker ejecta deposit. These regions are arranged asymmetrically around the crater: there are 4 discrete extended areas and none (or only a very faint one) to the southeast of the crater. These features in the reflectance difference map are consistent with the ejecta distribution map, so we are convinced that these maps represent the spatial distribution of the ejecta deposit and the ejected subsurface material has a lower reflectance than the average Ryugu surface. At 396 s and 498 s, the surface of Ryugu is visible through the ejecta curtain (Fig. 3, E andF; fig. S6, E andF), indicating that the optical thickness of the ejecta curtain was lower than at 192 s. A deep gray elliptical area appears, surrounded by the ejecta rays. We interpret this as the intrinsic color of the surface, and suggest that it represents the excavated area and/or the ejecta deposit area. Because crater growth is expected to have halted after 250 s (14), the ejecta curtain observed in these images must have been falling towards the surface on a ballistic trajectory. The major axis of the ellipse at 498 s was determined to be 36 m (14). Detachment of the whole ejecta curtain from the ground is generally expected during crater formation in the strength-dominated regime but not in the gravity-dominated regime [START_REF] Arakawa | Scientific objectives of Small Carry-on Impactor (SCI) and Deployable Camera 3 Digital (DCAM3-D): observation of an ejecta curtain and a crater formed on the surface of Ryugu by an artificial high-velocity impact[END_REF]. No detachment of the ejecta curtain was observed for the SCI crater formation, including in the final observation at 498 s (Fig. 3A-F). In Fig. 3 and Movie S2, the ejecta curtain appears to be deposited at the crater rim. This is supported by the ONC-T reflectance difference map in Fig. 4C, so the elevated rim surrounding the SCI crater could be (at least partially) composed of the ejecta deposits. This is expected for crater formation in the gravity-dominated regime [START_REF] Melosh | Impact cratering: A geologic process[END_REF][START_REF] Housen | Crater ejecta scaling laws: Fundamental forms based on dimensional analysis[END_REF]. We conclude that the structures of the SCI crater and the ejecta curtain indicate that they formed in the gravity-dominated regime.

For the crater formation to be dominated by gravity, the effective strength of the target must be less than 1.3 Pa (14). Cohesionless materials such as sand can satisfy that small strength at the depth shallower than 10 m, because the sub-surface shear strength originates from the frictional stress exerted on cohesionless materials and this stress is proportional to the surface gravity. We calculate the expected SCI crater radius using the conventional scaling law in the gravitydominated regime: 𝜋 $ = 𝑘 % • 𝜋 & '(/(&+() 𝜋 - (&+('./)/(0(&+()) [START_REF] Housen | Ejecta from impact craters[END_REF], where 𝜋 $ = 𝑅(𝜌 𝑚 ⁄ ) %/0 , 𝜋 & = (𝑔𝑎 𝑢 & ⁄ ), 𝜋 -= (𝜌 𝛿 ⁄ ), and 𝑅, 𝜌, 𝑚, 𝑎, 𝑢, 𝛿 are the apparent crater radius, the target density, the projectile mass, the projectile radius, the impact velocity, and the projectile density, respectively, and 𝑘 % , 𝜇, and 𝜈 are empirical parameters determined by laboratory experiments, which applies to typical cohesionless surfaces such as coarse or loose sand (14). The calculated expected radius is from 6.9 m to 7.7 m, which is close to the apparent radius of the observed SCI crater of 7.3 ± 0.4 m. We conclude the SCI crater formed on a cohesionless surface and that the broader surface of Ryugu must be composed of sand-like cohesionless material. According to the conventional scaling law, we can estimate the optimal 𝑘 % value applicable to the surface of Ryugu from the SCI crater radius R = 7.3 m (14); we obtain 𝑘 % = 0.62 if we assume that other parameters have typical values for sand, and can be extrapolated to cohesionless targets at high impact velocities: 𝜇 = 0.41 and 𝜈 = 0.4 (14, 18). We derive the crater size scaling law for the surface of Ryugu; 𝜋 $ = 0.62 • 𝜋 & '1.%3 𝜋 - 1.11%-. Our results affect the interpretation of the crater chronology on Ryugu. Based on collision frequency models for the asteroid Main Belt (19,[START_REF] O'brien | The collisional and dynamical evolution of the main-belt and NEA size distributions[END_REF], two estimates of the surface age of Ryugu have been obtained: (1.58±0.47)×10 8 years for a surface composed of dry-soils with cohesion or (8.9±2.5)×10 6 years for a cohesionless surface [START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF]. Our conclusion that the SCI crater was formed in the gravity-dominated regime on a cohesionless surface supports the younger age estimate. It also suggests that the crater retention age of the top 1 m of Ryugu's surface is ~1×10 5 years or younger, assuming crater production functions for near-Earth orbits (14,19,[START_REF] O'brien | The collisional and dynamical evolution of the main-belt and NEA size distributions[END_REF].

Supplementary Materials:

Materials and Methods

Figures S1-S11

Table S1 References [START_REF] Saiki | The small carry-on impactor (SCI) and the Hayabusa2 impact experiment[END_REF][START_REF] Ogawa | System configuration and operation plan of Hayabusa2 DCAM3-D camera system for scientific observation during SCI impact experiment[END_REF][START_REF] Ishibashi | Performance of Hayabusa2 DCAM3-D camera for short-range imaging of SCI and ejecta curtain generated from the artificial impact crater formed on asteroid 162137 Ryugu (1999JU3)[END_REF][START_REF] Gault | Experimental studies of oblique impact[END_REF][START_REF] Szeliski | Computer vision: Algorithms and applications[END_REF][START_REF] Schwanghart | Short communication: TopoToolbox2-MATLAB-based software for topographic analysis and modeling in Earth surface sciences[END_REF][START_REF] Al-Hashemi | A review on the angle of repose of granular materials[END_REF][START_REF] Maxwell | Simple Z model of cratering, ejection, and the overturned flap[END_REF][START_REF] Hermalyn | Scouring the surface: Ejecta dynamics and the LCROSS impact event[END_REF][START_REF] Yasui | Experimental study on impact-induced seismic wave propagation through granular materials[END_REF][START_REF] Tsujido | Ejecta velocity distribution of impact craters formed on quartz sand: Effect of projectile density on crater scaling law[END_REF][START_REF] Kiuchi | Experimental study on gravitational and atmospheric effects on crater size formed by low-velocity impacts into granular media[END_REF][START_REF] Barnouin | Impacts into coarse-grained spheres at moderate impact velocities: Implications for cratering on asteroids and planets[END_REF][START_REF] Güttler | Cratering experiments on the self armoring of coarsegrained granular targets[END_REF][START_REF] Takizawa | M. acknowledges funding support from the French space agency CNES and from Academies of Excellence: Complex systems and Space, environment, risk, and resilience, part of the IDEX JEDI of the Université Côte d'Azur. Author contributions: M.Ar. coordinated coauthor contributions; led the SCI and DCAM3 data acquisition, analyses, and interpretations; and wrote the paper, with contributions from K.W., 396 s (hyb2_dcam3d_img-f_01700[END_REF] Movies S1-S2 The typical reflectance factor of this area is about 0.018, thus the darkest regions in this map have decreased in reflectance by about 20% after the impact. Four discrete extended darkened areas correspond with the ejecta rays and are numbered 1 to 4. The same pre-impact shape model is used for both pre-and post-impact reflectance calculations (14), to emphasize weak reflectance changes due to ejecta deposit outside of the crater.

  K.O., T.Ka., K.S., K.I., R.H., N.S., P.M., M.J., and Y.S. SCI operations: T.Sa., K.W., T.Ka., K.S., Y.I., C.O., H.Y., Y.Taka., M.H., H.Im., and N.H. (Aizu) DCAM3 data acquisitions and reductions: K.W., K.O., K.S., H.Sa., K.I., R.H., N.S., Y.I., S.Kim., Y.M., T.To., S.N., and H.H. ONC data aqcuisitions and reductions: H.Sa., R.H., N.S., M.H., S.Su., T.M., S.Ka., E.T., Y.C., K.Yoshio., Y.Yo., M.M., M.Ya., T.Ko., and C.Ho. Image analysis: K.W., K.O., T.Ka., K.S., K.I., R.H., N.S., and Y.S. Spacecraft operations: T.Sa., K.O., K.S., H.Sa., R.H., N.S., H.Y., M.H., Y.M., S.N., Y.Ts., S.W., M.Yo., S.Tan., F.T., S.Kik., T.Y., N.O., G.O., K.Yoshiz., T.Ta., Y.Take., A.F., H.T., Y.Ya., T.O., C.Hi., S.H., O.M., T.Sh., S.So., R.T., T.I., M.O., M.Ab., N.N., K.K., S.Tac., and H.Ik. Shape modeling: K.S., Y.S., N.H. (Aizu), N.H. (Kobe), R.N., and A.M. Competing interests: Y.Ya. is also affiliated with Tokyo Metropolitan University. Data and materials availability: All images, and input data used in this study are available at the JAXA Data Archives and Transmission System (DARTS) at www.darts.isas.jaxa.jp/pub/ hayabusa2/paper/Arakawa_2020/. Additional data from Hayabusa2 ONC and DCAM3 data will be made available through the JAXA Data Archives and Transmission System website (https://www.darts.isas.jaxa.jp/).

Fig. 1 .

 1 Fig. 1. ONC images of the SCI crater. (A) The area around the SCI aiming point before the impact, taken from 1.72 km altitude (image hyb2_onc_20190321_192706_tvf_l2b. The Mobile Block (MB) and Stable Block (SB) are indicated with arrows. (B) The same location after the impact, observed from the same altitude (hyb2_onc_20190425_031226_tvf_l2b). (C) Orthocorrected context image of the impact location (hyb2_onc_20190516_023615_tvf_l2c) with the SCI crater shown by a yellow dashed semicircle. The white box indicates the region in panel D. (D) Close up image of the SCI crater taken from 115 m altitude (hyb2_onc_20190613_020217_tvf_l2c). A small pit is shown by the dashed curve. For panels A

Fig. 2 .

 2 Fig. 2. Topography of the SCI crater. (A) Digital Elevation Map of the SCI impact area. The black dashed semicircle shows the crater rim with diameter 17.6 m, and the red long dasheddotted semicircle shows the crater wall with diameter 14.5 m. The green dotted circle shows the central pit. The local height is measured relative to point F, which is labeled northwest of the

Fig. 3 .Fig. 4 .

 34 Fig. 3. DCAM3 images of the SCI impact. The images were taken at elapsed time after the impact of (A) 5 s (hyb2_dcam3d_img-f_01309_l2a), (B) 36 s (hyb2_dcam3d_img-f_01340_l2a), (C) 100 s (hyb2_dcam3d_img-f_01404_l2a), (D) 192 s (hyb2_dcam3d_img-f_01496_l2a), (E)