Masatoshi Hirabayashi

Ryota Nakano

Eri Tatsumi

Kevin J Walsh

Olivier S Barnouin

Patrick Michel

Christine S Hartzell

Daniel T Britt

Seiji Sugita

Sei-Ichiro Watanabe

Christine M Hartzell

William F Bottke

Daniel J Scheeres

Ronald-Louis Ballouz

Yuichiro Cho

Tomokatsu Morota

Ellen S Howell

Spin-driven evolution of asteroids' top-shapes at fast and slow spins seen

Keywords: Asteroids, Asteroids, rotation, Asteroids, surfaces, Near-Earth objects, Regolith 1. Introduction Small asteroids are

Ryugu is a top-shaped asteroid with an equivalent radius of 500 m and is currently rotating at a spin period of 7.6 h [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu -A spinning top-shaped rubble pile[END_REF]]. Ryugu's bulk density is similar to Bennu's,i.e.,1.19 g/cm 3 [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu -A spinning top-shaped rubble pile[END_REF], which implies that this asteroid is also a rubble pile. Observed OH-bearing materials (but at low level) [START_REF] Kitazato | The surface composition of asteroid 162173 Ryugu from Hayabusa2 near-infrared spectroscopy[END_REF] and highly porous boulders with low thermal inertia [START_REF] Okada | Highly porous nature of a primitive asteroid revealed by thermal imaging[END_REF] imply that Ryugu may be a remnant after catastrophic disruption of an undifferentiated parent body. Also, the Small Carry-on Impact experiment showed that the artificial crater formation is in the gravity regime, implying that the surface is covered by mechanically weak layers [Arakawa et al., 2020]. Unlike Bennu, however, Ryugu's slope does not reflect the current rotational condition. The slope variation becomes the smallest when the rotation period is 3.5 h -3.75 h. This finding led to an interpretation that there may have been reshaping processes at this spin period. It was found that the surface morphology is less cratered in the western region than that in the eastern region (the east-west dichotomy) [START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF]. The boulder distribution was also observed to be heterogeneous [START_REF] Michikami | Boulder size and shape distributions on asteroid Ryugu[END_REF], while mass movements may be limited and correlate with the gravitational potential [START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF]. Based on these observations, assuming that Ryugu's structure is uniform, [START_REF] Hirabayashi | Rotationally induced failure of irregularly shaped asteroids[END_REF] hypothesized that large deformation occurred at a spin period of 3.5 -3.75 h to produce the eastwest dichotomy.

These studies give some discrepancy of Bennu's and Ryugu's reshaping processes while they have top shapes with a similar bulk density. Specifically, these two bodies are considered to have been affected by rotationally driven deformation but may have evolved differently. Bennu's top shape may have evolved due to surface processing at present. The equatorial ridge may be developed due to mass wasting on the surface while the interior is intact (no mass flows occur in the interior). On the other hand, Ryugu's shape may result from large deformation, perhaps caused by internal failure. The interior becomes the most sensitive to structural failure at high spin, leading to internal mass flows (inward flows along the spin axis and outward flows in the radial direction) that make the shape oblate and push materials on the equatorial plane outward.

We fill a gap between these hypothesized reshaping processes. We revisit the FEM results from earlier studies [START_REF] Hirabayashi | Rotationally induced failure of irregularly shaped asteroids[END_REF][START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF][START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu -A spinning top-shaped rubble pile[END_REF] and explore the variations in the structural failure mode of these objects at different spin periods. To support our discussion, we newly employ a semi-analytical approach proposed by [START_REF] Nakano | Mass shedding activities of Asteroid (3200) Phaethon enhanced by its rotation[END_REF], who extended an approach by [START_REF] Hirabayashi | Failure modes and conditions of a cohesive, spherical body due to YORP spin-up[END_REF]. This semi-analytical approach computes Bennu's and Ryugu's yield conditions by assuming their shapes to be uniformly rotating triaxial ellipsoids and characterizes their global failure modes. We finally acknowledge [START_REF] Cheng | Surface boulders as clues for reconstructing the formation history of topshaped asteroids[END_REF], who recently proposed that top-shaped asteroids like Ryugu and Bennu result from mass movements in shallow surface layers at a fast rotation period. They assumed a large, rigid core in a target sphere to reproduce a top shape. In the present study, the semi-analytical approach and the FEM approach discuss the variations in the structural failure modes under the uniform structure condition.

We outline the present manuscript. Section 2 briefly reviews the proposed evolution processes of top-shaped asteroids. Section 3 shows the analysis approaches employed in this manuscript: the plastic FEM approach and the semi-analytical approach. Section 4 summarizes the results from our analyses. The results from the FEM approach are already introduced in [START_REF] Hirabayashi | Rotationally induced failure of irregularly shaped asteroids[END_REF], [START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF], and [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu -A spinning top-shaped rubble pile[END_REF], and we compare these results in this manuscript. Section 5 discusses the evolution of top-shaped asteroids based on the results from Section 4. Our FEM approach assumes the structure of a body to be uniform. On the other hand, the semi-analytical approach provides the distribution of the minimum cohesive strength that can avoid structure failure (see the definition in Section 3) and give insights into the dependence of the failure mode on the internal structure.

Rotational evolution of top-shaped asteroids

While Section 1 summarized the rotational evolution of top-shaped bodies by focusing on Bennu and Ryugu, this section reviews the general background of such objects. Ground-based radar observations show that spheroidal shapes are dominant in observed objects [e.g., Taylor et al., 2012;Benner et al., 2015]. Among them are reported to be top shapes, some of which also have smaller satellites [e.g., Ostro et al., 2006;[START_REF] Brozović | Radar and optical observations and physical modeling of triple near-Earth Asteroid (136617[END_REF]Becker et al., 2015;[START_REF] Busch | Radar observations and the shape of near-Earth asteroid[END_REF][START_REF] Naidu | Radar imaging and characterization of the binary near-earth asteroid (185851) 2000 DP107[END_REF][START_REF] Naidu | Radar observations and a physical model of binary near-Earth asteroid 65803 Didymos, target of the DART mission[END_REF][START_REF] Walsh | Formation and Evolution of Binary Asteroids[END_REF]. Some top shapes are rotating close to their spin limits at which surface materials may be about to be shed [START_REF] Scheeres | Dynamical Configuration of Binary Near-Earth Asteroid[END_REF][START_REF] Harris | On the shapes and spins of "rubble pile[END_REF], while others like Ryugu are not [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu -A spinning top-shaped rubble pile[END_REF]. Findings of top shapes in different taxonomy types imply that the evolution of top shapes are not linked to the material compositions of asteroids [e.g., [START_REF] Walsh | Rubble Pile Asteroids[END_REF]. After the formation of top shapes and their precursor bodies during re-accumulation after catastrophic disruption [START_REF] Michel | Formation of top-shaped asteroids by collisional disruption: implications for the origin of Ryugu and Bennu[END_REF], these bodies may have evolved by experiencing fast rotation. There are two competing scenarios for rotationally driven reshaping processes that may have induced top shapes.

First, surface processing may be the primary contributor to developing top shapes [e.g., [START_REF] Harris | On the shapes and spins of "rubble pile[END_REF][START_REF] Minton | The topographic limits of gravitationally bound, rotating sand piles[END_REF][START_REF] Scheeres | Landslides and Mass shedding on spinning spheroidal asteroids[END_REF][START_REF] Walsh | Rotational breakup as the origin of small binary asteroids[END_REF][START_REF] Walsh | Spin-up of rubble-pile asteroids: Disruption, satellite formation, and equilibrium shapes[END_REF][START_REF] Cheng | Surface boulders as clues for reconstructing the formation history of topshaped asteroids[END_REF]. When the centrifugal force becomes dominant at fast rotation, surface mass movements may occur towards lower potential regions around the equator, while the interior is intact [e.g., [START_REF] Sánchez | Asteroid Evolution: Role of Geotechnical Properties, Asteroids: New Observations[END_REF]Scheeres, 2018, 2020]. The equatorial ridge evolves by the accumulation of surface materials from middle and high latitude regions. When materials slide down to the equatorial region, they may experience the Coriolis force and move toward the longitudinal direction [START_REF] Harris | On the shapes and spins of "rubble pile[END_REF][START_REF] Statler | Mechanism of self-reinforcing YORP acceleration for fast-rotating asteroids, Asteroids, Comets[END_REF]. Also, when surface processing is significant, the contacts of particles may tend to induce more mass ejection on irregular surfaces [START_REF] Yu | The Dynamical Complexity of Surface Mass Shedding from a Top-shaped Asteroid Near the Critical Spin Limit[END_REF].

Second, if the stress field in the interior is considered, another top-shape evolution process may be possible [START_REF] Hirabayashi | Stress and Failure Analysis of Rapidly Rotating Asteroid (29075) 1950 DA[END_REF]. With fast rotation, they should have higher stress condition around the center of mass [START_REF] Hirabayashi | Failure modes and conditions of a cohesive, spherical body due to YORP spin-up[END_REF]. This mechanism results from force vectors acting on each element that is affected by gravity and rotation. In a top shape, the direction of the gravity force is radially inward, while that of the centrifugal force is horizontally outward. When the body is rotating at fast spin, the centrifugal force at lower latitudes exceeds the gravity force, leading to a horizontal, outward loading. At higher latitudes, on the other hand, the inward component of the gravitational force is always dominant. These loading modes induce outward deformation in the horizontal direction around the equator and inward deformation in the vertical direction near the spin axis, possibly enhancing the evolution process of a top shape.

A stand-alone explanation of each scenario sounds incomplete to discuss the top-shape evolution mechanism. If surface processing is a main driver that develops top shapes at fast spin, the interior should be mechanically stronger than the surface layers. On the other hand, if the mechanical strength of the interior is similar to that of the surface layer, the stress field in the interior should reach its yield condition before surface regions experience mass movements [START_REF] Hirabayashi | Failure modes and conditions of a cohesive, spherical body due to YORP spin-up[END_REF][START_REF] Hirabayashi | Stress and Failure Analysis of Rapidly Rotating Asteroid (29075) 1950 DA[END_REF]. Gravity measurements implied that Ryugu's structure may be homogeneous [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu -A spinning top-shaped rubble pile[END_REF]. The present study attempts to explain the correlation of surface and internal failure by considering the structural conditions at different spin periods.

Analysis approach 3.1. Results from the earlier plastic FEM approach

A plastic FEM approach was developed and designed to analyze the locations of structurally failed regions in an irregularly shaped body that is uniformly rotating. The developed code was compatible with ANSYS Mechanical APDL 18.1. In this study, instead of conducting new FEM analyses, we revisit the structural analysis results for Bennu and Ryugu from our earlier works [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu -A spinning top-shaped rubble pile[END_REF][START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF][START_REF] Hirabayashi | Rotationally induced failure of irregularly shaped asteroids[END_REF]. We show how Bennu and Ryugu experience rotationally driven failure modes at multiple spin periods when the structure is assumed to be uniform. While the present work does not provide new simulations, we briefly introduce our FEM simulation settings below. Note that the detailed model development is also provided in the earlier work [e.g., [START_REF] Hirabayashi | Fission and reconfiguration of bilobate comets as revealed by 67P/Churyumov-Gerasimenko[END_REF][START_REF] Hirabayashi | Rotationally induced failure of irregularly shaped asteroids[END_REF].

The FEM meshes of Bennu and Ryugu were generated from the mission-derived polyhedron shape models [Barnouin et al., 2019;[START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu -A spinning top-shaped rubble pile[END_REF]. We first used TetGen [START_REF] Si | TetGen, a Delaunay-based Quality Tetrahedral mesh Generator[END_REF] to convert polyhedron shape models to 4-node triangular FEM meshes. Then, these FEM meshes were modified to obtain 10-node meshes. Because the meshes developed by TetGen were not optimally compatible with the simulation environment in ANSYS, we refined them by performing manual operations and by using the ANSYS mesh refining functions. The derived mesh was further processed to provide the constraints on the spatial movements of nodes. To mimic the condition of an irregularly shaped body in space, we gave constraints on three degree of freedom at the center of gravity for the translational motion and another three degrees of freedom at two surface nodes for the rotational motion [START_REF] Hirabayashi | Rotationally induced failure of irregularly shaped asteroids[END_REF]. We computed a force loading acting on each node by accounting for the self-gravity force and the centrifugal force.

Our model assumes that structural deformation is small and follows linear elasticity below the yield condition and perfect plasticity on that condition. Perfect plasticity specifies that there is no material hardening and softening. To define the yield condition, we used the Drucker-Prager yield criterion, which is written as [Chen and Han, 1988] 𝑓 = 𝛼𝐼 ! + &𝐽 " -𝑠 ≤ 0.

(1) 𝐼 ! and 𝐽 " are the stress invariants, which are given as

𝐼 ! = 𝜎 ! + 𝜎 " + 𝜎 # , (2)
𝐽 " = 1 6 {(𝜎 ! -𝜎 " ) " + (𝜎 " -𝜎 # ) " + (𝜎 # -𝜎 ! ) " }, ( 3
) where 𝜎 $ (𝑖 = 1, 2, 3) is the principal stress component. 𝛼 and 𝑠 are material constants that are defined as

𝛼 = 2 sin 𝜙 √3(3 -sin 𝜙) , 𝑠 = 6𝑌 cos 𝜙 √3(3 -sin 𝜙) , ( 4
) where 𝜙 and 𝑌 are the angle of friction and the cohesive strength, respectively. We note that 𝑌 is defined in the sense of continuum media, which is the cumulative effect of cohesion caused by grain-grain interactions in a given small element.

In our FEM simulations, the structural condition was assumed to be uniform over the entire body. Young's modulus (𝐸 % ), Poisson's ratio (𝜇), and the angle of friction (𝜙) were fixed at 1 × 10 & Pa, 0.25, and 35 ∘ , respectively [e.g., [START_REF] Lambe | Soil Mechanics[END_REF]. For simplicity, no dilatancy was allowed in these simulations. The necessary parameters for FEM simulations are available in Table 1. To determine the structural failure mode of an irregularly shaped body at a given spin period, we changed the cohesive strength under the constant parameters. If the actual cohesive strength is enough to resist failure, the body can keep the original shape. However, as the cohesive strength becomes smaller, some regions in the body start failing structurally. We computed the distribution of failure and the minimum cohesive strength that the body can avoid failure anywhere (later denoted as 𝑌 * ) [START_REF] Hirabayashi | Rotationally induced failure of irregularly shaped asteroids[END_REF].

[m 3 ] Ryugu's dimension [2𝑎, 2𝑏, 2𝑐] 1,004 × 1,004 × 876 [m 3 ]

Semi-analytical stress model

While our plastic FEM approach can describe the failure modes of irregularly shaped bodies, the results of structural failure depend on local topographic features [START_REF] Hirabayashi | Rotationally induced failure of irregularly shaped asteroids[END_REF]. Such local features may be too detailed to study global conditions that induce either surface processing or internal failure. Here, to characterize the mechanisms of such global failure modes, we apply a semi-analytical approach that computes the variations in the stress distributions in a uniformly rotating triaxial ellipsoid at different spin periods. Note that the model assumption of a rotating triaxial ellipsoid gives a reasonable approximation for the global failure modes of a topshaped asteroid. The present approach is an extension of the previous approach, which focused on a spherical body [START_REF] Hirabayashi | Failure modes and conditions of a cohesive, spherical body due to YORP spin-up[END_REF]. This approach is first discussed by [START_REF] Nakano | Mass shedding activities of Asteroid (3200) Phaethon enhanced by its rotation[END_REF].

The semi-analytical approach consists of two aspects. The first aspect is to compute the stress distribution in a uniformly triaxial ellipsoid by assuming linear elasticity [START_REF] Dobrovolskis | Internal stresses in phobos and other triaxial bodies[END_REF][START_REF] Holsapple | Equilibrium Configurations of Solid Cohesionless Bodies[END_REF][START_REF] Love | A Treatise on the Mathematical Theory of Elasticity[END_REF]. The second aspect is to use this stress distribution and the Drucker-Prager yield criterion to derive the minimum cohesive strength, 𝑌 ) * , at given elements. We distinguish 𝑌 ) * with 𝑌 * , which is given as the minimum cohesive strength that a considered body starts to experience failure somewhere (from Section 3.1). Importantly, we define the following relationship: 𝑌 * = max (𝑌 ) * ). This approach does not provide the deformation processes but shows the sensitivity of a body element to structural failure. We confirmed that the earlier semi-analytical approach and the FEM approach gave consistent results for the failure modes of a uniformly rotating spherical body [START_REF] Hirabayashi | Failure modes and conditions of a cohesive, spherical body due to YORP spin-up[END_REF]. Below, we describe how to compute 𝑌 ) * . The details are also found in [START_REF] Nakano | Mass shedding activities of Asteroid (3200) Phaethon enhanced by its rotation[END_REF].

For the first aspect, the stress distribution in a uniformly rotating triaxial ellipsoid is computed based on linear elasticity [START_REF] Dobrovolskis | Internal stresses in phobos and other triaxial bodies[END_REF][START_REF] Holsapple | Equilibrium Configurations of Solid Cohesionless Bodies[END_REF][START_REF] Love | A Treatise on the Mathematical Theory of Elasticity[END_REF]. Considering the equilibrium state of stress, we can write the stress distribution as

𝜕𝜎 !! 𝜕𝑥 ! + 𝜕𝜎 !" 𝜕𝑥 " + 𝜕𝜎 !# 𝜕𝑥 # + 𝜌𝑏 ! = 0, (5) 𝜕𝜎 "! 𝜕𝑥 ! + 𝜕𝜎 "" 𝜕𝑥 " + 𝜕𝜎 "# 𝜕𝑥 # + 𝜌𝑏 " = 0, (6) 𝜕𝜎 #! 𝜕𝑥 ! + 𝜕𝜎 #" 𝜕𝑥 " + 𝜕𝜎 ## 𝜕𝑥 # + 𝜌𝑏 # = 0, ( 7
)
where 𝜎 $* (𝑖, 𝑗 = 1,2,3) is the stress component, 𝜌 is the bulk density, which is assumed to be uniform over the entire volume, and 𝑏 $ (𝑖 = 1,2,3) is the body force, which includes a centrifugal force and a self-gravity force [START_REF] Holsapple | Equilibrium Configurations of Solid Cohesionless Bodies[END_REF]. In these equations, 𝜎 !" = 𝜎 "! , 𝜎 #" = 𝜎 "# , and 𝜎 !# = 𝜎 #! . Because a triaxial ellipsoid is axisymmetric, we define the 𝑥 ! , 𝑥 " , and 𝑥 # axes along the minimum, intermediate, and maximum moment of inertia axes, respectively.

From Equations ( 5) through ( 7), the displacement 𝑢 $ in the 𝑥 $ direction is expressed with 12 parameters [START_REF] Dobrovolskis | Internal stresses in phobos and other triaxial bodies[END_REF],

𝑢 ! = 𝑥 ! W𝐴 + 𝐵 𝑥 ! " 𝑎 " + 𝐶 𝑥 " " 𝑏 " + 𝐷 𝑥 # " 𝑐 " \ , (8)
𝑢 " = 𝑥 " W𝐸 + 𝐹 𝑥 ! " 𝑎 " + 𝐺 𝑥 " " 𝑏 " + 𝐻 𝑥 # " 𝑐 " \ , (9)
𝑢 # = 𝑥 # W𝐼 + 𝐽 𝑥 ! " 𝑎 " + 𝐾 𝑥 " " 𝑏 " + 𝐿 𝑥 # "
𝑐 " \ , (10) where 𝑎, 𝑏, and 𝑐 are the semi-major axis, the semi-intermediate axis, and the semi-minor axis, respectively (𝑎 ≥ 𝑏 ≥ 𝑐). Also, the upper alphabets indicate coefficients that are determined by the constitutive law and the boundary condition. Using Equations ( 8) through (10), we can describe the strain components as

𝜖 $* = 1 2 d 𝜕𝑢 $ 𝜕𝑥 * + 𝜕𝑢 * 𝜕𝑥 $ e .
(11) We next consider Hooke's law (linear elasticity), which is given as

𝜎 $* = 𝜆𝜖 ++ 𝛿 $* + 2𝜇𝜖 $* , ( 12
)
where 𝜖 ++ = 𝜖 !! + 𝜖 "" + 𝜖 ## and 𝛿 $* is the Kronecker delta. 𝜆 and 𝜇 are Lame's constants, which are given as

𝜆 = 𝐸 % 𝜈 (1 + 𝜈)(1 -2𝜈) , ( 13
) 𝜇 = 𝐸 % 2(1 + 𝜈) . ( 14
)
The traction boundary condition is given as 𝜎 $* 𝑛 * = 0, (15) where 𝑛 * is written as

𝑛 ! = 𝑥 ! 𝑎 " 𝜒 , 𝑛 " = 𝑥 " 𝑏 " 𝜒 , 𝑛 # = 𝑥 # 𝑐 " 𝜒 , ( 16
) and

𝜒 = d 𝑥 ! " 𝑎 , + 𝑥 " " 𝑏 , + 𝑥 # " 𝑐 , e ! "
. ( 17) Combining all these conditions, we obtain 12 linear equations to uniquely determine the 12 coefficients appearing in Equations ( 8) through (10).

For the second aspect, we convert these stress components to the principal stress components, 𝜎 ! , 𝜎 " , and 𝜎 # . Then, we substitute these principal components into Equation ( 1). Using the constant parameters defined in Table 1, we compute the distribution of 𝑌 ) * in a uniformly triaxial ellipsoid at different spin periods. Specifically, 𝑌 ) * is computed by considering 𝑓 = 0 in Equation ( 1).

Comparison of the structural conditions between Bennu and Ryugu 4.1. Results from finite element analysis approach

We first summarize the failure modes of Bennu and Ryugu at different spin periods by introducing the earlier results from our plastic FEM approach [START_REF] Hirabayashi | Rotationally induced failure of irregularly shaped asteroids[END_REF][START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF][START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu -A spinning top-shaped rubble pile[END_REF]. The present study shows the dependence of the failure mode on the spin period, which was also predicted for irregularly shaped bodies by [START_REF] Hirabayashi | Rotationally induced failure of irregularly shaped asteroids[END_REF]. If Bennu and Ryugu are rotating at longer spin periods, their failure mode is characterized as surface failure. On the other hand, if they are rotating at shorter spin periods, internal failure becomes significant. We describe the contrasts of these processes by referring to the plastic FEM results [START_REF] Hirabayashi | Rotationally induced failure of irregularly shaped asteroids[END_REF][START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF][START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu -A spinning top-shaped rubble pile[END_REF].

Figure 2 shows the failure modes of Bennu at spin periods of 4.3 h (the current spin period) and 3.5 h [START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF]. At the 4.3 h spin period, the interior does not fail structurally, while the surface regions at middle and high latitudes are sensitive to structural failure (Panels a and b). This failure mode appears when the cohesive strength is ~0.2 Pa. If the cohesive strength of surface materials is below this level, mass movements should occur. The flow processes observed by OSIRIS-REx [Barnouin et al., 2019;[START_REF] Scheeres | The geophysical environment of Bennu[END_REF][START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF][START_REF] Walsh | Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface[END_REF][START_REF] Jawin | The geology of (101955) Bennu from the first year of OSIRIS-REx observations, 51th Lunar and Planetary Science Conference[END_REF] imply that tin top-layers may be covered with cohesionless materials. Thus, our results support that Bennu's top shape may have evolved by surface processing at the present spin period. At the 3.5 h spin period, on the other hand, Bennu should experience internal failure if the cohesion is less than 0.8 Pa. The central region is the most sensitive to structural failure; the predicted failed region is almost axisymmetric and reaches the top surface.

Ryugu's structural behavior is similar to Bennu's, although the variations in the failed regions appear due to local topography at some level (Figure 3). At a spin period of 7.6 h (the current spin period), the structure does not reach the yield condition in the almost entire volume when the cohesive strength is fixed at 1 Pa (Panels a and b). This trend results from self-gravity, which leads the interior to have compression with limited shear and thus not to experience structural failure significantly. At a spin period of 3.5 h, internal failure becomes significant [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu -A spinning top-shaped rubble pile[END_REF]. Unlike the symmetric failed region appeared in Bennu at the 3.5 h spin period, Ryugu's failed region is asymmetric (Panels c and d) [START_REF] Hirabayashi | Rotationally induced failure of irregularly shaped asteroids[END_REF]. Structural failure is observed to be more concentrated in the central area and in the bottom-right part of the plots. This asymmetric distribution of the failed region may imply that in the current shape, Ryugu's internal structure at the 3.5 h spin period is divided into a structurally relaxed region (the upper-left part) and an unrelaxed region (the bottom-right part).

From these structural conditions at different spin periods, the evolution process of Ryugu's top shape was interpreted to be different from that of Bennu's. Ryugu's top shape may have evolved by large deformation when this asteroid was rotating at a spin period of 3.5 ~ 3.75 h [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu -A spinning top-shaped rubble pile[END_REF]. Because of Ryugu's east-west dichotomy, the surface in the western region is smoother than that in the eastern region [START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF]. Also, the equatorial ridge in the western region is sharper than that in the eastern region [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu -A spinning top-shaped rubble pile[END_REF]. The locations of the failed regions predicted from the plastic FEM approach correspond to those of the eastern region [START_REF] Hirabayashi | Rotationally induced failure of irregularly shaped asteroids[END_REF]. Therefore, the relatively less cratered condition of the western region implies that this region may be structurally relaxed due to a large deformation process at a spin period of 3.5 h ~ 3.75 h [START_REF] Hirabayashi | Rotationally induced failure of irregularly shaped asteroids[END_REF]. We note that when Ryugu's spin has slowed down to the current spin period, 7.6 h, the body may have also had surface processing as Bennus has.

Results from semi-analytical approach

The FEM approach gave a significant interpretation that Bennu's top shape may result from surface processing at the current spin period, 4.3 h, while Ryugu's may be developed by large deformation, which is due to internal failure, at a spin period of 3.5 h ~ 3.75 h. This contrast between Bennu's and Ryugu's failure modes implies that top shapes may be able to evolve by both surface processing and internal deformation. This section uses the semi-analytical approach, which can describe global failure modes, to further explore how these failure modes appear at different spin periods. We compute the distribution of 𝑌 ) * over slices of Bennu and Ryugu along their spin axis. We incorporate the parameter conditions defined in Table 1 into the semi-analytical model.

First, we discuss the distribution of 𝑌 ) * for Bennu in Figure 4. We observe that the distribution of 𝑌 ) * varies as a function of the spin period. When the spin period is shorter than 3.5 h, 𝑌 ) * becomes the maximum at the center of the body. At a spin period of 3 h, for example, 𝑌 ) * at the central region reaches higher than 12 Pa. 𝑌 ) * on the surface in the equatorial region is lower than that in the central region, but still higher than 0 Pa. In the polar regions, 𝑌 ) * becomes 0 Pa, implying that there is no strength necessary to avoid structural failure. This result indicates that the interior is the most sensitive to structural failure if the actual cohesive strength is uniform over the body. At a spin period longer than 3.75 h, nonzero 𝑌 ) * is concentrated in the equatorial surface region but does not appear in the central region. As the spin period continues to become longer, the regions with nonzero 𝑌 ) * become smaller. At a spin period of 4.3 h (the current spin period), nonzero 𝑌 ) * only appears in shallow surface regions at the equator. In this region, 𝑌 ) * is less than 1 Pa. This result indicates that if the actual cohesive strength is less than that level, surface materials should slide down to regions at lower latitudes. This interpretation is consistent with what was observed from the plastic FEM approach (Figure 2).

Next, we show the distribution of 𝑌 ) * for Ryugu. The bulk density of Ryugu is the same as that of Bennu. Their global shapes are also similar. Therefore, we anticipate that Ryugu's failure mode is similar to Bennu's. Figure 5 illustrates the failure modes of Ryugu at spin periods from 2.5 h to 7.6 h. The magnitude of 𝑌 ) * for Ryugu is higher than that for Bennu because of the size difference. At a spin period shorter than 3.5 h, the region around the center needs higher 𝑌 ) * to resist structural failure. However, when the spin period is around 3.75 h, the region with nonzero 𝑌 ) * is separated towards the equatorial surface regions. As the spin period becomes longer, this region shrinks and eventually disappears. At a spin period of 7.6 h, there is no obvious region with nonzero 𝑌 ) * . This result implies that there is no significant deformation at the current spin period. Again, the results are consistent with those from the plastic FEM approach (Figure 3).

Our semi-analytical model describes that if the structure is uniform, structural failure in Ryugu and Bennu can be characterized by two modes. The first mode is that the interior structurally fails first at fast spin. This mode can induce vertical compression along the spin axis and radial deformation in the horizonal direction. The second mode is that surface processing is dominant, and top-surface layers at low and middle latitudes fail structurally at low spin. For these failure modes, the size of internal failure is larger than that of surface failure. As the spin period becomes longer, the transition occurs when nonzero 𝑌 * disappears from the central region and is split towards the surface regions. While these modes are apparently different, they induce mass movements towards the equatorial region and thus can generate top shapes. The similar finding was presented by [START_REF] Hirabayashi | Failure modes and conditions of a cohesive, spherical body due to YORP spin-up[END_REF] for a limited case in which a body is perfectly spherical.

Discussion

Surface processing and internal deformation can develop a top-shaped body

The results in Section 4 showed that both surface processes and internal deformation can contribute to the evolution of Bennu's and Ryugu's top shapes. At a longer rotation period, surface layers are the most sensitive to structural failure, while an interior is structurally intact. If surface layers do not have enough cohesive strength, they should be modified. As the spin period becomes shorter, the area of elements sensitive to structural failure increases. Eventually, the central region becomes the most sensitive. At this stage, the interior fails structurally, leading to a deformation mode that consists of a horizontal, outward flow and a vertical, inward flow. The development of a top shape due to these two processes was confirmed by studies that used Discrete Element Model [e.g., [START_REF] Hirabayashi | Stress and Failure Analysis of Rapidly Rotating Asteroid (29075) 1950 DA[END_REF][START_REF] Sánchez | Disruption patterns of rotating self-gravitating aggregates: A survey on angle of friction and tensile strength[END_REF], 2018, 2020;[START_REF] Walsh | Rotational breakup as the origin of small binary asteroids[END_REF][START_REF] Walsh | Spin-up of rubble-pile asteroids: Disruption, satellite formation, and equilibrium shapes[END_REF]Yu et al., 2020;[START_REF] Zhang | Creep stability of the proposed AIDA mission target 65803 Didymos: I. Discrete cohesionless granular physics model[END_REF][START_REF] Zhang | Rotational Failure of Rubble-pile Bodies: Influences of Shear and Cohesive Strength[END_REF].

The dependence of deformation on the spin period controls the evolution of a top shape. Because the centrifugal force is the main contributor to reshaping, the shape evolves so that its moment of inertia along the spin axis always increases (Section 4). Because the angular momentum is constant during such a deformation process, the spin period becomes longer. On the other hand, YORP may accelerate/decelerate the spin state, depending on the shape and how these effects have affected the body [e.g., Bottke et al., 2015;[START_REF] Statler | Extreme sensitivity of the YORP effect to small-scale topography[END_REF][START_REF] Statler | Obliquity of "top-shaped" asteroids may not imply reshaping by YORP spinup[END_REF]. If top shapes spin up, they again become sensitive to another deformation process.

The magnitude of reshaping may be controlled by the evolution of the spin period. In Figure 6, we show the schematics of possible paths of the evolution process of a top shape. When the body is spinning at a spin period at which internal failure occurs (the upper path), deformation is large. This deformation process leads to a higher moment of inertia under constant angular momentum, making the spin period longer. After this event, deformation stops. To have another deformation process, the body needs to be rotationally accelerated for a long time. On the other hand, if the body is spinning at relatively slow spin, surface processes can still occur (the lower path), and the magnitude of deformation may be small. Because the rotational change may be small during such an event, it is not necessary for the body to wait for a long spin-up duration to have another small deformation event. This implies that surface processing events occur more frequently than internal deformation. The paths can interact with each other. This evolution process indicates a complex coupling relationship between the shape and the spin period evolution [START_REF] Cotto-Figueroa | Coupled spin and shape evolution of small rubble-pile asteroids: self-limitation of the YORP effect[END_REF]. Also, we point out that deformation is not necessarily axisymmetric. Structural and shape heterogeneity even at small level can induce asymmetric deformation processes [START_REF] Hirabayashi | Stress and Failure Analysis of Rapidly Rotating Asteroid (29075) 1950 DA[END_REF][START_REF] Sánchez | Asteroid Evolution: Role of Geotechnical Properties, Asteroids: New Observations[END_REF][START_REF] Sánchez | Disruption patterns of rotating self-gravitating aggregates: A survey on angle of friction and tensile strength[END_REF]. When a structurally weak location experiences deformation, that region should be relaxed structurally. Then, the body that had deformation in one area now has structural sensitivities in other areas, which may deform later.

Bennu has complex trends of mass movements [START_REF] Walsh | Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface[END_REF][START_REF] Jawin | The geology of (101955) Bennu from the first year of OSIRIS-REx observations, 51th Lunar and Planetary Science Conference[END_REF] and unique morphologies such as grooves [Barnouin et al., 2019]. The energy analysis showed that the boundary of the Roche lobe at the current spin period would be consistent with the surface slope variation observed [START_REF] Scheeres | The geophysical environment of Bennu[END_REF][START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF]. These findings imply that Bennu's surface topography is currently evolving. The present study implies that at the current spin period, Bennu's interior should be structurally intact, while its surface is sensitive to structural failure (Figure 4). By "intact", we mean that 𝑌 ) * in the interior is zero, thus deformation driven by rotation may not occur. This condition leads to surface processing from the middle latitudes to the equatorial regions if the strength is low. Each mass movement event is small and asymmetric. Because the change at the spin period is small, deformation processes occur occasionally. Also, Bennu's not-perfectlycircular equatorial ridge imply that the spin period has not reached 3.5 h, and the ridge has evolved within the gravity force-dominant region, i.e. the inside of the Roche lobe. If Bennu had experienced a shorter spin period than 3.5 h, the ridge would have been modified to be circular due to a strong centrifugal force [START_REF] Scheeres | The geophysical environment of Bennu[END_REF][START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF]].

Ryugu's top shape, on the other hand, may have evolved due to large deformation. This idea was supported by its almost-perfectly-circular equatorial ridge. [START_REF] Watanabe | Hayabusa2 arrives at the carbonaceous asteroid 162173 Ryugu -A spinning top-shaped rubble pile[END_REF] pointed out that such a circular ridge might indicate the centrifugal force-gravity force boundary that has occurred at a critical spin period, 3.5 h. At this spin period, the interior of this asteroid becomes the most sensitive to structural failure (Section 4). Hayabusa2 showed that the western region looks apparently younger than the eastern region [START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF]. This east-west dichotomy implies different time scales for the reshaping process [START_REF] Hirabayashi | Rotationally induced failure of irregularly shaped asteroids[END_REF]. First, a large-scaled deformation may have occurred to create the global feature of Ryugu's top shape (such as its equatorial ridge). Then, another large-scaled deformation may have happened in the western region and developed its young, smooth terrain.

From these analyses, we conclude that both internal deformation and surface processing can lead to top shapes. After re-accumulation processes due to catastrophic disruptions [START_REF] Michel | Formation of top-shaped asteroids by collisional disruption: implications for the origin of Ryugu and Bennu[END_REF], precursor bodies experience these processes at different spin periods in complex ways to evolve into top shapes. Again, it was previously argued that top shapes would be attributed to one of these processes only.

Insights into reshaping and top shapes' interiors from surface morphologies

While the present study assumed that the structure of Ryugu and Bennu are homogeneous, small bodies may have local heterogeneity at some level, given their formation and evolution processes. Many bodies smaller than 100 km in diameter may result from the collisional evolution of larger bodies [Bottke et al., 2005[Bottke et al., , 2015;;[START_REF] Morbidelli | Asteroids were born big[END_REF][START_REF] Walsh | Rubble Pile Asteroids[END_REF]. The following reaccumulation processes may cause the heterogeneities of rubble pile asteroids [START_REF] Michel | Formation of top-shaped asteroids by collisional disruption: implications for the origin of Ryugu and Bennu[END_REF].

Even after re-accumulation, deformation processes can change their volumes, changing the density distribution [e.g., [START_REF] Hirabayashi | Stress and Failure Analysis of Rapidly Rotating Asteroid (29075) 1950 DA[END_REF][START_REF] Scheeres | The dynamic geophysical environment of (101955) Bennu based on OSIRIS-REx measurements[END_REF]Bagatin et al., 2020]. Granular convections may modify the structure of asteroids [e.g., [START_REF] Murdoch | Granular Convection in Microgravity[END_REF]Yamada and Katsuragi, 2014]. Seismic shaking may be a driver that induces granular convection. We note that seismic waves may be significantly attenuated in granular media in a microgravity environment [START_REF] Yasui | Experimental study on impact-induced seismic wave propagation through granular materials[END_REF]Matsue et al., 2020], as seen by the Hayabusa2 SCI impact experiment [Arakawa et al., 2020[START_REF] Nishiyama | An absence of seismic shaking on Ryugu induced by the impact experiment on the Hayabusa2 mission, 51th Lunar and[END_REF]. Hypervelocity impact processes also resurface and reorganize the interior, depending on their energy levels [e.g., [START_REF] Wünnemann | A strain-based porosity model for use in hydrocode simulations of impacts and implications for transient crater growth in porous targets[END_REF].

The failure mode of a top-shaped asteroid depends on such structural heterogeneities [START_REF] Hirabayashi | Failure modes and conditions of a cohesive, spherical body due to YORP spin-up[END_REF]. Based on the distributions of 𝑌 ) * in Ryugu and Bennu (Figures 4 and5), Figure 7 describes the dependence of the failure mode on the cohesive strength. This schematic was first introduced by [START_REF] Hirabayashi | Failure modes and conditions of a cohesive, spherical body due to YORP spin-up[END_REF]; in the present study, we also add the case of a weak interior proposed by [START_REF] Sánchez | Rotational evolution of self-gravitating aggregates with cores of variable strength[END_REF]. If the structure is uniform, both surface processing and internal failure can contribute to reshaping. When structurally weaker layers are placed on top of a stronger interior, only surface particles can move at high spin [START_REF] Hirabayashi | Failure modes and conditions of a cohesive, spherical body due to YORP spin-up[END_REF][START_REF] Sánchez | Disruption patterns of rotating self-gravitating aggregates: A survey on angle of friction and tensile strength[END_REF], 2018, 2020;Yu et al., 2020]. If the interior is weaker than the surface layer, complex reshaping may occur to generate irregularly shaped bodies rather than top shapes [START_REF] Sánchez | Rotational evolution of self-gravitating aggregates with cores of variable strength[END_REF].

The observed crater morphologies may imply local heterogeneity of Bennu and Ryugu. On Bennu, some small craters exhibit the discontinuities of the boulder size distribution inside and outside the rims [START_REF] Walsh | Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface[END_REF]. Similar features are also observed on Ryugu [START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF][START_REF] Watanabe | Paucity of boulders in shallow craters on asteroid[END_REF][START_REF] Morota | Surface evolution of carbonaceous asteroid 162173 Ryugu revealed from global mapping and touchdown operation of[END_REF]. Small particles are distributed inside such craters, while the outside exhibits larger boulders [START_REF] Walsh | Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface[END_REF]. One explanation may be grain segregation due to granular convection such as the Brazil nut effect [START_REF] Watanabe | Paucity of boulders in shallow craters on asteroid[END_REF]. Another explanation may be that small particles move on top surface layers more favorably than large boulders and are depleted [START_REF] Szalay | Dust Phenomena Relating to Airless Bodies[END_REF]. If this is the case, top surface layers can be considered to be mechanically weak and have evolved rapidly. Such transportation may be driven not only by combinations of the centrifugal force and the gravitational force but also by other mechanisms. An example may be electrostatic lofting, the amount of which is dependent on cohesion and large electrostatic forces driven by photoemission or secondary electron emission [START_REF] Hartzell | Dynamics of 2D electrostatic dust levitation at asteroids[END_REF][START_REF] Hartzell | An Evaluation of Electrostatic Lofting as an Active Mechanism on Bennu[END_REF][START_REF] Wang | Dust charging and transport on airless planetary bodies[END_REF][START_REF] Zimmerman | Grain-scale supercharging and breakdown on airless regoliths[END_REF]. Micrometeoroid impacts [e.g., [START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF][START_REF] Walsh | Craters, boulders and regolith of (101955) Bennu indicative of an old and dynamic surface[END_REF]] and thermal cracking [e.g., Moralo et al., 2015, 2017;[START_REF] Hazeli | The origins of Asteroidal rock disaggregation: Interplay of thermal fatigue and microstructure[END_REF] may also result in small particle transportation.

Necessary technical improvements

We finally discuss that the current study has two issues that need to be solved in our future work. First, our discussions assumed that the body has a uniform structure. However, as discussed in Section 5.2, there may exist some structural heterogeneity. If this is the case, we expect the variations in the internal structure. Such variations affect our simulation settings, e.g., the friction angle and the density distribution. Second, our model does not give detailed local failure modes, which may highly correlate with the timescale that materials have been exposed to space weathering. In carbonaceous asteroids, space weathering may make materials bluer and redder, depending on their conditions [START_REF] Gillis-Davis | Incremental laser space weathering of Allende reveals non-linear like space weathering effects[END_REF][START_REF] Hendrix | C-complex asteroids: UV-visible spectral characteristics and implications for space weathering effects[END_REF][START_REF] Lantz | Ion irradiation of carbonaceous chondrites: A new view of space weathering on primitive asteroids[END_REF][START_REF] Matsuoka | Pulse-laser irradiation experiments of Murchison CM2 chondrite for reproducing space weathering on C-type asteroids[END_REF][START_REF] Thompson | Spectral and chemical effects of simulated space weathering of Murchison CM2 carbonaceous chondrite[END_REF]. Recent studies pointed out that on both Bennu and Ryugu, space weathering may redden materials [START_REF] Sugita | The geomorphology, color, and thermal properties of Ryugu: Implications for parent-body processes[END_REF][START_REF] Morota | Surface evolution of carbonaceous asteroid 162173 Ryugu revealed from global mapping and touchdown operation of[END_REF][START_REF] Dellaguistina | Properties of rubble-pile asteroid (101955) Bennu from OSIRIS-REx imaging and thermal analysis[END_REF][START_REF] Dellaguistina | Reflectance and color heterogeneity on asteroid Bennu, 51st Lunar and Planetary Science Conference[END_REF] while space weathering on Bennu needs further interpretations. However, the distribution of bluer materials is highly localized, and our models are too coarse to provide consistent results. These issues need to be resolved to give stronger insights into the evolution mechanisms of top-shaped asteroids.

Conclusion

The present study explored the rotationally driven evolution of Bennu's and Ryugu's top shapes. We started by revisiting the previously reported results derived from a plastic FEM approach. The semi-analytical approach was then introduced to describe the minimum cohesive strength distribution that can avoid structural failure at a given element, under the assumption that a top shape is a uniformly rotating triaxial ellipsoid. We then quantified how the failure modes of Bennu and Ryugu would vary at different spin periods. The results showed that the failure modes depended on the spin period. For these bodies, surface regions would be sensitive to structural failure at a long spin period, while interiors could be the most vulnerable at a higher spin period. Therefore, we suggest that top shapes can evolve by various reshaping processes that depend on both rotation and internal structure. Bottke, W. F. et al. (2005), Linking the collisional history of the main asteroid belt to its dynamical excitation of depletion, Icarus, 179, 1, p. 63-94, doi:10.101663-94, doi:10. /j.icarus.200563-94, doi:10. .05.017. Bottke, W. F. et al. (2015)), In search of the source of asteroid (101955) Bennu: Applications of the stochastic YORP model, Icarus, 247, p.191 -217, doi:10.1016Icarus, 247, p.191 -217, doi:10. /j.icarus.2014.09.046.09.046.

Figures and Tables

  80NSSC19K0548) and Auburn University/Intramural Grant Program. Also, the FEM results shown in this manuscript are obtained by using ANSYS Mechanical APDL 18.1. Reference Arakawa, M. et al. (2020), An artificial impact on the asteroid 162173 Ryugu formed a crater in the gravity-dominated regime, Accepted in Science. Bagatin, A. C. et al. (2020), Gravitational re-accumulation as the origin of most contact binaries and other small body shapes, Icarus, 339, 113603, doi:10.1016/j.icarus.2019.113603. Barnouin, O. S. et al. (2019), Shape of (101955) Bennu indicative of a rubble pile with internal stiffness, Nature Geosciences, 12, p.247-252, doi:10.1038/s41561-019-0330-x. Becker, T. M. et al. (2015), Physical modeling of triple near-Earth Asteroid (153591) 2001 SN263 from radar and optical light curve observations, Icarus, 248, p.499-515, doi:10.1016/j.icarus.2014.10.048. Benner, L. A. M. et al. (2015), Radar observations of near-Earth and main-belt asteroids. In Asteroid IV, p.165 -182, University of Arizona, Tucson, doi:10.2458/azu-uapress_9780816532131-ch009.

1 . 25 Figure 5 .

 1255 Figure 1. Images of Bennu (Panel a) and Ryugu (Panel b). Panel a. Image credit: NASA/Goddard/University of Arizona. Panel b. Image credit: JAXA/UTokyo/Kochi U./Rikkyo U./Nagoya U./Chiba Ins. Tech/Meiji U./U. Aizu/AIST.

Table 1 . Physical parameter settings in simulations. The dimensions are from Barnouin et al. [2019] for Bennu and Watanabe et al. [2019] for Ryugu.

 1

	Physical Parameters	Symbols	Value	Units
	Young's modulus	𝐸 %	1 × 10 &	[Pa]
	Poisson's ratio	𝜈	0.25	[-]
	Angle of friction	𝜙	35	[deg]
	Bulk density	𝜌	1.19	[g/cm3]
	Bennu's dimension	[2𝑎, 2𝑏, 2𝑐] 506 × 492 × 457


Nina Miolane
email: nmiolane@stanford.edu

Nicolas Guigui
email: nicolas.guigui@inria.fr

Johan Mathe

Benjamin Hou
email: benjamin.hou11@imperial.ac.uk

Yann Thanwerdas
email: yann.thanwerdas@inria.fr

Stefan Heyder
email: stefan.heyder@tu-ilmenau.de

Olivier Peltre
email: opeltre@gmail.com

Niklas Koep
email: niklas.koep@gmail.com

Hadi Zaatiti
email: hadi.zaatiti@irt-systemx.fr

Hatem Hajri
email: hatem.hajri@irt-systemx.fr

Yann Cabanes
email: yann.cabanes@gmail.com

Thomas Gerald
email: thomas.gerald@lip6.fr

Paul Chauchat
email: pchauchat@gmail.com

Christian Shewmake
email: cshewmake2@gmail.com

Daniel Brooks

Bernhard Kainz
email: b.kainz@imperial.ac.uk

Claire Donnat
email: cdonnat@stanford.edu

Susan Holmes

Xavier Pennec
email: xavier.pennec@inria.fr

Soon Cheng

Ong

Alice Le Brigant
email: alice.le-brigant@univ-paris1.fr

Geomstats: A Python Package for Riemannian Geometry in Machine Learning

Keywords: differential geometry, Riemannian geometry, statistics, machine learning, manifold

We introduce Geomstats, an open-source Python package for computations and statistics on nonlinear manifolds such as hyperbolic spaces, spaces of symmetric positive definite matrices, Lie groups of transformations, and many more. We provide object-oriented and extensively unit-tested implementations. Manifolds come equipped with families of Riemannian metrics with associated exponential and logarithmic maps, geodesics, and parallel transport. Statistics and learning algorithms provide methods for estimation, clustering, and dimension reduction on manifolds. All associated operations are vectorized for batch computation and provide support for different execution backends-namely NumPy, Py-Torch, and TensorFlow. This paper presents the package, compares it with related libraries, and provides relevant code examples. We show that Geomstats provides reliable building blocks to both foster research in differential geometry and statistics and democratize the use of Riemannian geometry in machine learning applications.

Introduction

Data on manifolds naturally arise in different fields. Hyperspheres model directional data in molecular and protein biology [START_REF] Kent | Using the Fisher-Bingham distribution in stochastic models for protein structure[END_REF] and are used through Principal Nested Spheres to model some aspects of 3D shapes [START_REF] Jung | Analysis of principal nested spheres[END_REF][START_REF] Hong | Non-Euclidean Classification of Medically Imaged Objects via s-reps[END_REF]. Computations on hyperbolic spaces arise for impedance density estimation [START_REF] Huckemann | Möbius deconvolution on the hyperbolic plane with application to impedance density estimation[END_REF], geometric network comparison [START_REF] Asta | Geometric Network Comparison[END_REF], and the analysis of reflection coefficients extracted from a radar signal [START_REF] Chevallier | Probability density estimation on the hyperbolic space applied to radar processing[END_REF]. Symmetric positive definite (SPD) matrices characterize data from diffusion tensor imaging (DTI) [START_REF] Pennec | A Riemannian Framework for Tensor Computing[END_REF][START_REF] Yuan | Local polynomial regression for symmetric positive definite matrices[END_REF] and functional magnetic resonance imaging (fMRI) [START_REF] Sporns | The human connectome: A structural description of the human brain[END_REF]. Covariance matrices, which are also SPD matrices, appear in applications such as automatic speech recognition (ASR) systems [START_REF] Shinohara | Covariance clustering on Riemannian manifolds for acoustic model compression[END_REF], image and video descriptors [START_REF] Mehrtash | Sparse Coding on Symmetric Positive Definite Manifolds using Bregman Divergences[END_REF], or air traffic complexity representation [START_REF] Le | Quantization and clustering on riemannian manifolds with an application to air traffic analysis[END_REF]. The Lie groups of transformations SO(3) and SE(3) appear naturally when dealing with articulated objects like the human spine [START_REF] Arsigny | Processing Data in Lie Groups: An Algebraic Approach. Application to Non-Linear Registration and Diffusion Tensor MRI[END_REF][START_REF] Boisvert | Principal spine shape deformation modes using riemannian geometry and articulated models[END_REF] or the pose of a camera [START_REF] Kendall | Geometric Loss Functions for Camera Pose Regression with Deep Learning[END_REF][START_REF] Hou | Computing CNN loss and gradients for pose estimation with Riemannian geometry[END_REF]. Stiefel manifolds are used to process video action data or to analyze two vector-cardiograms [START_REF] Chakraborty | Statistics on the (compact) Stiefel manifold : Theory and Applications[END_REF]. Grassmannians appear in computer vision to perform video-based face recognition and shape recognition [START_REF] Turaga | Statistical analysis on Stiefel and Grassmann Manifolds with applications in Computer Vision[END_REF]. Statistics on landmark spaces are used in anthropology and many applied fields to describe biological shapes [START_REF] Richtsmeier | Advances in Anthropological Morphometrics[END_REF]. More generally, Kendall shape spaces gave rise to a very important literature in shape statistics [START_REF] Dryden | Statistical shape analysis, with Applications in R[END_REF]. A variety of applications use tools from infinite-dimensional Riemannian geometry to study the shapes of discretized curves such as those sampled from closed two or three-dimensional curves defining the contours of organs in computational anatomy [START_REF] Younes | Spaces and manifolds of shapes in computer vision: An overview[END_REF]. In addition, open curves that describe the temporal evolution of physical phenomena are embedded in various manifolds, for example in a Lie group [START_REF] Celledoni | Shape analysis on Lie groups with applications in computer animation[END_REF] or in the hyperbolic plane [START_REF] Le | Computing distances and geodesics between manifold-valued curves in the SRV framework[END_REF]).

Yet, the adoption of methods from differential geometry has been inhibited by the lack of a reference implementation. Code sequences are often custom-tailored for specific problems and are not easily reused. Some Python packages do exist, but they often focus on optimization: Pymanopt [START_REF] Townsend | Pymanopt: A python toolbox for optimization on manifolds using automatic differentiation[END_REF], Geoopt [START_REF] Bécigneul | Riemannian Adaptive Optimization Methods[END_REF][START_REF] Kochurov | Geoopt: Riemannian Adaptive Optimization Methods with pytorch optim[END_REF]McTorch (Meghwanshi et al., 2018). Others are dedicated to a single manifold: PyRiemann on SPD matrices [START_REF] Barachant | PyRiemann: Python package for covariance matrices manipulation and Biosignal classification with application in Brain Computer interface[END_REF], PyQuaternion on 3D rotations [START_REF] Wynn | PyQuaternions: A fully featured, pythonic library for representing and using quaternions[END_REF], and PyGeometry on spheres, toruses, 3D rotations and translations [START_REF] Censi | PyGeometry: Library for handling various differentiable manifolds[END_REF]. Lastly, others lack unit-tests and continuous integration: TheanoGeometry (Kühnel and Sommer, 2017). There is a need for an open-source implementation of differential geometry and associated learning algorithms for manifold-valued data.

We present Geomstats, an open-source Python package for computations and statistics on nonlinear manifolds. Geomstats has three main objectives: (i) foster research in differential geometry and geometric statistics by providing low-level code to gain intuition or test a theorem and a platform to share algorithms; (ii) democratize the use of geometric statistics by implementing user-friendly geometric learning algorithms using Scikit-Learn API; and (iii) provide educational support to learn "hands-on" differential geometry and geometric statistics, through its examples, notebooks and visualizations.

Implementation Overview

The package geomstats is organized into two main modules: geometry and learning. The module geometry implements concepts in Riemannian geometry with an objectoriented approach. Manifolds mentioned in the introduction are available as classes that inherit from the base class Manifold. The base class RiemannianMetric provides methods such as the geodesic distance between two points, the exponential and logarithm maps at a base point, etc. As examples, HyperbolicMetric, StiefelCanonicalMetric, the Lie groups' InvariantMetric s, or the curves' L2Metric and SRVMetric [START_REF] Srivastava | Shape analysis of elastic curves in euclidean spaces[END_REF] inherit from RiemannianMetric. Going beyond Riemannian geometry, the class Connection implements affine connections. The implementation uses automatic differentiation with autograd to allow computations on manifolds when closed-form formulae do not exist. The API of the module geometry uses terms from differential geometry to foster contributions of researchers from this field. The module learning implements statistics and learning algorithms for data on manifolds. The code is object-oriented and classes inherit from Scikit-Learn base classes and mixins: BaseEstimator, ClassifierMixin, RegressorMixin, etc. This module provides classes such as FrechetMean, KMeans, TangentPCA, for implementations of Fréchet mean estimators [START_REF] Fréchet | Les éléments aléatoires de nature quelconque dans un espace distancié[END_REF], K-means, and principal component analysis (PCA) designed for manifold data. The API of the module learning follows Scikit-Learn's API, being therefore user-friendly to machine learning researchers and engineers.

The code follows international standards for readability and ease of collaboration, is vectorized for batch computations, undergoes unit-testing with continuous integration using Travis, relies on either Numpy, TensorFlow or PyTorch backends. The package comes with a visualization module to provide intuition on differential geometry (see Figure 1), and with a datasets module that provides toy data sets on manifolds. The repositories examples and notebooks provide convenient starting points to get familiar with geomstats.

The GitHub repository at github.com/geomstats/geomstats offers a convenient way to ask for help or request features by raising issues. The website geomstats.github.io provides documentation for users and important guidelines for those wishing to contribute to the project.

Comparison and Interaction with Existing Packages

This section compares the package geomstats with related Python implementations on differential geometry and learning. Table 1 compares the geometric operations and Table 2 compares the engineering infrastructures.

The library TheanoGeometry [START_REF] Kühnel | Computational Anatomy in Theano[END_REF] is the most closely related to Geomstats and provides nonlinear statistics and stochastic equations on Riemannian manifolds. The differential geometric tensors are computed with automatic differentiation. However, this library does not provide statistical learning algorithms and lacks engineering maintenance. Several other packages focus on optimization on Riemannian manifolds. Pymanopt [START_REF] Townsend | Pymanopt: A python toolbox for optimization on manifolds using automatic differentiation[END_REF] computes gradients and Hessian-vector products on Riemannian manifolds with automatic differentiation and provides the following solvers: steepest descent, conjugate gradient, the Nelder-Mead algorithm, particle swarm optimization, and the Riemannian trust regions. Geoopt [START_REF] Kochurov | Geoopt: Riemannian Adaptive Optimization Methods with pytorch optim[END_REF] focuses on stochastic adap- tive optimization on Riemannian manifolds, for machine learning problems. The library provides stochastic solvers, stochastic gradient descent and Adam, as well as the following samplers: Stochastic Gradient Langevin Dynamics, Hamiltonian Monte-Carlo, Stochastic Gradient Hamiltonian Monte-Carlo. Lastly, McTorch [START_REF] Meghwanshi | McTorch, a manifold optimization library for deep learning[END_REF] provides optimization on Riemannian manifold for deep learning by adding a "Manifold" parameter to PyTorch's network layers and optimizers. The library provides the following solvers: stochastic gradient descent, AdaGrad and conjugate gradients. As these libraries focus on optimization, they substitute potentially computationally expensive operations by practical proxies, for example, by replacing exponential maps by so-called retractions. However, they are less modular than Geomstats in terms of the Riemannian geometry. For example, each manifold comes with a single Riemannian metric, in contrast to Geomstats where families of Riemannian metrics are implemented. Furthermore, they do not provide statistical learning algorithms.

The optimization libraries are complementary to Geomstats and interact easily with it. Geomstats provides low-level implementations of Riemannian geometry that can be used to define optimization costs. In turn, an optimization library can provide an efficient solver to use within the implementation of Geomstats' learning algorithms. An example of such interactions, between Pymanopt and Geomstats can be found in Geomstats' examples folder.

Usage: Examples of Learning on Riemannian manifolds

Three steps are needed to run learning algorithms on manifolds with Geomstats: (i) instantiate the manifold of interest, (ii) instantiate the learning algorithm of interest, and (iii) run the algorithm. The following code snippet illustrates the use of K-means on the hypersphere. s p h e r e = Hypersphere ( dim=5) ) c l u s t e r i n g = OnlineKMeans ( m e t r i c=s p h e r e . m e t r i c , n _ c l u s t e r s =4) c l u s t e r i n g = c l u s t e r i n g . f i t ( data )

The following code snippet shows the use of tangent PCA on the 3D rotations. All geometric computations are performed behind the scenes. The user only needs a high-level understanding of Riemannian geometry. Each algorithm can be used with any of the manifolds and metric implemented in the package. The folders examples and notebooks provide many more code snippets that help users get started with Geomstats.

Conclusion

We presented the Python Package Geomstats, with the aim to provide the wider Machine Learning community with off-the-shelf geometric learning algorithms on a wide variety of manifolds, and flexibility in the choice of metrics, while being faithful to the mathematician's formulation of Riemannian geometry. This sometimes comes at cost of efficiency, and future contributions will be devoted to addressing this caveat.

Figure 1 :

 1 Figure 1: Left: Minimization of a scalar field on the sphere S 2 using Riemannian gradient descent. Middle: Regular geodesic grid on the hyperbolic space H 2 in Poincaré disk representation. Right: Geodesic on the Lie group SE(3) for the canonical left-invariant metric. These examples and more are available at geomstats.ai.

  s o 3 = S p e c i a l O r t h o g o n a l ( n=3, point_type= ' v e c t o r ' ) m e t r i c = s o 3 . b i _ i n v a r i a n t _ m e t r i c data = s o 3 . random_uniform ( n_samples =10) t p c a = TangentPCA ( m e t r i c=m e t r i c , n_components=2) t p c a = t p c a . f i t ( data ) t a n g e n t _ p r o j e c t e d _ d a t a = t p c a . t r a n s f o r m ( data )

Table 1 :

 1 Comparison of libraries in terms of geometric operations

		Manifolds	Geometry
	Pymanopt	Euclidean manifold, symmetric matri-ces, sphere, complex circle, SO(n), Stiefel, Grassmannian, oblique mani-fold, SPD(n), elliptope, fixed-rank PSD matrices	Exponential maps, retraction, vector transport, and logarithmic egrad2rgrad, ehess2rhess, inner product, distance, norm
	Geoopt	Euclidean manifold, sphere, Stiefel, Poincaré ball	Same as Pymanopt
	McTorch	Stiefel, SPD(n)	Same as Pymanopt
			Inner product, exponential and log-
	TheanoGeometry	Sphere, ellipsoid, SPD(n), Landmarks, GL(n), SO(n), SE(n)	arithmic maps, parallel transport, Christoffel symbols, Riemann, Ricci and scalar curvature, geodesics, Fréchet
			mean
	Geomstats	Euclidean, space, Poincaré polydisk, hypersphere, Minkowski, hyperbolic SO(n), SE(n), GL(n), Stiefel, Grass-mannian, SPD(n), symmetric matrices, skew-symmetric matrices, discretized curves on manifolds, landmarks on manifolds	Levi-Civita connection, Christoffel symbols, parallel transport, expo-nential and logarithmic maps, inner product, distance, norm, geodesics, group invariant metrics, Fréchet means and learning algorithms on manifolds
		Backends	Continuous integration (CI) and coverage
	Pymanopt	Autograd, PyTorch, TensorFlow, Theano CI, coverage 85%
	Geoopt	PyTorch	75%
	McTorch	PyTorch	CI, coverage 84%
	TheanoGeometry Theano	No CI, no unit tests
	Geomstats	NumPy, PyTorch, TensorFlow	CI, coverage 92% (NumPy), 76% (TensorFlow), 79% (PyTorch)

Table 2 :

 2 Comparison of code infrastructure data = s p h e r e . random_uniform ( n_samples =10

Acknowledgments

This work is partially supported by the National Science Foundation, grant NSF DMS RTG 1501767, the Inria-Stanford associated team GeomStats, the European Research Council (ERC) under the EU Horizon 2020 research and innovation program (grant agreement G-Statistics No. 786854) and by the French Government through the 3IA Côte d'Azur Investments in the Future project (National Research Agency ANR-19-P3IA-0002).