SPH Simulations for Shape Deformation of Rubble-Pile Asteroids Through Spinup: The Challenge for Making Top-Shaped Asteroids Ryugu and Bennu

Asteroid Ryugu and asteroid Bennu, which were recently visited by spacecraft Hayabusa2 and OSIRIS-REx, respectively, are spinning top-shaped rubble piles. Other axisymmetric top-shaped near-Earth asteroids have been observed with ground-based radar, most of which rotate near breakup rotation periods of ~ 3 hours. This suggests that rotation-induced deformation of asteroids through rotational spinup produces top shapes. Although some previous simulations using the Discrete Element Method showed that spinup of rubble piles may produce oblate top shapes, it is still unclear what kinds of conditions such as friction angles of constituent materials and spinup timescales are required for top-shape formation. Here we show, through Smoothed Particle Hydrodynamics simulations of granular bodies spinning-up at different rates, that the rotation-induced deformation of spherical rubble piles before breakup can be classified into three modes according to the friction angle \phi_{d}: quasi-static and internal deformation for \phi_{d}<40 degrees, dynamical and internal deformation for 50 degrees<\phi_{d}<60 degrees, and surface landslides for \phi_{d}>70 degrees. Note that these apparent large values of friction angle can be acceptable if we consider the effect of cohesion among blocks of a rubble pile under weak gravity. Bodies with \phi_{d}<60 degrees evolve into oblate spheroids through internal deformation, but never form pronounced equators defining a top shape. In contrast, bodies with \phi_{d}>70 degrees deform into axisymmetric top shapes through an axisymmetric surface landslides if spinup timescales are1 month, bodies with \phi_{d}>70 degrees deform into non-axisymmetric shapes via localized landslides. We suggest that rapid spinup mechanisms are preferable for the formation of axisymmetric top shapes.


Abstract
Asteroid 162173 Ryugu and asteroid 101955 Bennu, which were recently visited by spacecraft Hayabusa2 and OSIRIS-REx, respectively, are spinning topshaped rubble piles. Other axisymmetric top-shaped near-Earth asteroids have been observed with ground-based radar, most of which rotate near breakup rotation periods of ∼ 3 hours. This suggests that rotation-induced deformation of asteroids through rotational spinup produces top shapes. Although some previous simulations using the Discrete Element Method showed that spinup of rubble piles may produce oblate top shapes, it is still unclear what kinds of conditions such as friction angles of constituent materials and spinup timescales are required for top-shape formation. Here we show, through Smoothed Particle Hydrodynamics simulations of granular bodies spinning-up at different rates, that the rotation-induced deformation of spherical rubble piles before breakup can be classified into three modes according to the friction angle φ d : quasi-static and internal deformation for φ d ≤ 40 • , dynamical and internal deformation for 50 • ≤ φ d ≤ 60 • , and a set of surface landslides for φ d ≥ 70 • . Note that these apparent large values of friction angle can be acceptable if we consider the effect of cohesion among blocks of a rubble pile under weak gravity. Bodies with φ d ≤ 60 • evolve into oblate spheroids through internal deformation, but never form pronounced equators defining a top shape. In contrast, bodies with φ d ≥ 70 • deform into axisymmetric top shapes through an axisymmetric set of surface landslides if spinup timescales are a few days. In addition, through slow spinups with timescales 1 month, bodies with φ d ≥ 70 • deform into non-axisymmetric shapes via localized sets of landslides. We suggest that rapid spinup mechanisms are preferable for the formation of axisymmetric top shapes.

Introduction
The Japan Aerospace Exploration Agency (JAXA) spacecraft Hayabusa2 arrived at asteroid 162173 Ryugu in June 2018 and investigated its detailed physical and geometrical properties via remote-sensing instruments. The proximity observation revealed that Ryugu is a rubble-pile asteroid with a so-called "spinning-top" shape having a prominent equatorial ridge (Watanabe et al. 2019). Radar observations already found some top-shaped near-Earth asteroids such as 1999 KW4, 2001 SN263, and65803 Didymos (Ostro et al. 2006;Becker et al. 2015;Benner et al. 2010), and the recent proximity observation with the NASA spacecraft OSIRIS-REx confirmed that 101955 Bennu also has a top shape (Lauretta et al. 2019; Barnouin et al. 2019). Top-shaped asteroids have several common geometrical characteristics: polar-to-equatorial axis ratios of 0.8, pronounced equatorial ridges, nearly conical surfaces extending from the equators to the mid-latitudes, and almost axisymmetric shapes.
The fast rotation of a rubble-pile asteroid results in its deformation or breakup because the centrifugal force exceeds self-gravity at the equatorial radius, and the critical rotation period at which the centrifugal force equals selfgravity at the equator is 2π/ (4/3)πGρ ≈ 3.3 h using a typical bulk density of C-type asteroids ρ ≈ 1 g/cm 3 (Carry 2012), where G is the gravitational constant. The spin periods of the top-shaped asteroids 1999 KW4, 2001 SN263, Didymos, and Bennu are 2.8 h, 3.4 h, 2.3 h, and 4.3 h, respectively (Ostro et al. 2006;Becker et al. 2015;Richardson et al. 2016;Lauretta et al. 2019). Although the current rotation period of Ryugu is 7.6 h (Watanabe et al. 2019), surface slope distribution and failure analysis of Ryugu suggest that Ryugu once had a rotation period of ∼ 3.5 h (Watanabe et al. 2019;Hirabayashi et al. 2019). This strongly suggests that Ryugu was reshaped by centrifugally induced deformation during the period of rapid rotation (Watanabe et al. 2019). Thus, other top-shaped asteroids may have also formed through deformation due to rapid rotation.
All the top-shaped asteroids found so far are near-Earth asteroids and have diameters several kilometers (Ostro et al. 2006;Benner et al. 2010;Brozović et al. 2011;Busch et al. 2011;Rozitis et al. 2014;Becker et al. 2015;Lauretta et al. 2019;Watanabe et al. 2019). Their top shapes were identified through radar or proximity observations, which are mainly possible only for near-Earth asteroids. About half of the top-shaped asteroids have satellites (2001NE263:Becker et al. 2015, 1999KW4:Ostro et al. 2006, 1994CC:Brozović et al. 2011, Didymos:Benner et al. 2010. The orbits and sizes of the satellites and the fraction of the top-shaped asteroids with satellites may further constrain the formation process of top shapes. Torque induced by the absorption of sunlight and its thermal re-emission changes the rotation of asteroids. The spinup or spindown mechanism is called the Yarkovsky-O'Keefe-Radzievskii-Paddack (YORP) effect (e.g., Rubincam 2000;Walsh 2018), and the YORP spinup of several asteroids has already been confirmed (Taylor et al. 2007;Kaasalainen et al. 2007;Hergenrother et al. 2019). On the other hand, impacts with other small asteroids also change rotation states. Merging collisions and accretion of fragments after catastrophic collisions tend to increase spin rates due to angular momentum supply (e.g., Sugiura et al. 2018Sugiura et al. , 2020Michel et al. 2020). The spinup of asteroids caused by such mechanisms eventually leads to global deformation, and sometimes probably results in top shapes . The spinup also induces breakup of bodies to form satellites, which is consistent with the omnipresence of satellites around top-shaped asteroids.
How rapidly rotating rubble piles deform has been investigated in various studies through analytical approaches (Holsapple 2001(Holsapple , 2004Sharma et al. 2009;Holsapple 2010;Sharma 2013Sharma , 2018Hirabayashi 2015) and numerical simulations (Walsh et al. , 2012Sánchez & Scheeres 2016;Zhang et al. 2017Zhang et al. , 2018Sánchez & Scheeres 2018;Leisner et al. 2020). Most of the previous simulations about spinup of bodies have been conducted using the Discrete Element Methods (DEMs), in which usually spherical particles subject to gravity and contact forces represent the rubble-pile constituents. The simulations of spinup show that a body with a closely packed particle configuration experiences surface landslides and deforms into a top shape, while a body with a loosely packed particle configuration experiences internal deformation Zhang et al. 2017). This suggests that the interlocking of particles may affect deformation modes. DEM simulations correctly represent bulk friction of rubble piles if the sizes of DEM particles are comparable to those of blocks or regolith particles of asteroids. Practical sizes of DEM particles for asteroid simulations, however, are much larger than those of regolith particles due to computational limitations. The difference of the particle sizes possibly affects details of deformation modes of spinning-up bodies. Using a plastic finite element model (FEM) technique, some researches analyzed detailed relationship between cohesive strength and failure modes of rapidly rotating asteroids (e.g., Hirabayashi 2015;, while they could not investigate resultant shapes produced through the failure due to the limitation of FEM. Although the axisymmetric top-shaped figures are shown to be in equilibrium as a granular body (Harris et al. 2009), the formation process for such bodies is still unclear. Thus, the quantitative conditions of friction angle and acceleration rate required to form top shapes are still unclear.
The Smoothed Particle Hydrodynamics (SPH) method was recently applied to numerical simulations of shape deformation of asteroids (Jutzi 2015). In the SPH simulations, a rubble-pile body is considered as a continuum of granular material and we explicitly set a specific value of friction angle of the material. Thus, we expect that the SPH method has the advantage of precise control of bulk friction of rubble piles. So far the SPH method has mainly been used to investigate shapes of impact outcomes (e.g., Jutzi & Asphaug 2015;Jutzi & Benz 2017;Leleu et al. 2018;Sugiura et al. 2018Sugiura et al. , 2019Sugiura et al. , 2020, and this is the first work applying the SPH method to the rotational deformation of small bodies.
In this paper, we present the results of SPH simulations about spinups of spherical rubble-pile bodies with various values of friction angles of constituent material. We then examine the friction angles required for the formation of top shapes. To clarify the spinup rate dependence, we also present the results of simulations with a slower spinup rate. Although the adopted rate is even much faster than the typical YORP spinup rate owing to the limitation of computational time, the comparison of these simulations with different spinup rates will help us to consider what kinds of acceleration mechanisms are appropriate for the top-shape formation.

Method
The numerical code utilized in this study is based on a version of the SPH code described in Sugiura et al. (2018Sugiura et al. ( , 2019. The code solves hydrodynamic equations (e.g., the equations of continuity and motion) with self-gravity and yield processes with the aid of a friction model that determines shear strength of granular material (Jutzi 2015). For the surface of a granular body, the time evolution of the density distribution of the body is calculated by solving the equation of continuity. The code is parallelized for distributed memory computers by using the Framework for Developing Particle Simulator (Iwasawa et al. 2015;Iwasawa et al. 2016). Sugiura (2020) verified the validity of our SPH approach for granular materials through the comparison with laboratory experiments of cliff collapse (Lajeunesse et al. 2005).
We investigated the rotational deformation of a kilometer-sized granular body, where compressional stresses derived from self-gravity are too small for granular material to experience plastic compression. This enables us to use the following elastic equation of state: where p is the pressure, ρ is the bulk density, C s is the bulk sound speed, and ρ 0 is the bulk density at uncompressed states. We set ρ 0 = 1.19 g/cm 3 , which is the common value of the measured bulk densities of Ryugu and Bennu (Watanabe et al. 2019;Lauretta et al. 2019). We used a typical sound speed of granular material ∼ 100 m/s for C s (Teramoto & Yano 2005). Simulation results only weakly depend on the adopted sound speed or the equations of state (see Appendix A). Note that we ignored the direct effect of cohesion, i.e., tensile forces, and thus we set p = 0 if p becomes negative.

Simulation setup and analysis method
We mainly investigated the rotational deformation of a uniform and spherical body with a radius of 500 m, which is similar to Ryugu's equatorial radius (Watanabe et al. 2019). Note that the deformation process and resultant shape are almost independent of the size of the body in our simulations (see Appendix A). We briefly mention the dependence of initial shapes on the deformation process by simulating bodies slightly different from the sphere in Section 3.3. The total number of SPH particles comprising the body is about 25,000 (or ≈ 20 SPH particles in the radial direction), which has a sufficient resolution to capture shape-related features of a resultant body. The SPH particles were initially not placed on a regular lattice, but distributed randomly around the lattice points with a uniform density based on a procedure described by Sugiura et al. (2018).
In general, the shear strength Y d of cohesive granular material as a function of the confining pressure p can be represented as Y d = tan(φ d,real )p+c coh , where φ d,real and c coh are the "real" friction angle and the cohesive strength, respectively. For simplicity, we adopted the relation in a form Y d = tan(φ d )p with the "effective" friction angle φ d by interpreting that the relation φ d = tan −1 (Y d /p) would be valid for cohesive granular materials. The effective friction angle φ d used in our simulations can be expected to mimic an aspect of cohesion. The cohesion of a subsurface layer of asteroid Ryugu is estimated to be up to c coh = 670 Pa from an artificial impact experiment done in the Hayabusa2 mission (Arakawa et al. 2020), whereas the central pressure of a spherical body with a radius of R = 500 m and mean density of ρ 0 = 1.19 g/cm 3 is only p c = (2/3)πGρ 2 0 R 2 ≈ 50 Pa, so that the effective friction angle is up to φ d ≈ tan −1 (c coh /p c ) ≈ 86 • . To mimic large effective friction angles due to cohesion, we varied φ d from 20 • to 80 • . This treatment allows us to understand the detailed dependence of deformation modes on one parameter φ d . The validity of using such large effective friction angles instead of introducing the cohesion c coh is discussed in Section 4.1.
We adopted an inertial coordinate system with origin at the center of mass of the body and the z-axis directed along the initial rotation axis of the body. The initial spherical shape of the body is stable for slow rotations, and the minimum rotation period P min required for keeping the initial configuration gets shorter for larger φ d (Holsapple 2001). We set the initial rotation period of the body in each simulation with an effective friction angle φ d to be longer than P min for the corresponding φ d . To represent the spinup process of a body, we accelerated the rotation of the body around the z-axis as follows. We added an increment of velocity ∆v i in the rotational direction to the i-th SPH particle comprising the body at each numerical step as whereω 0 is the predetermined angular acceleration of rotation, ∆t ≈ 0.1 s is the time step, e z is the unit vector parallel to the z-axis, and r i is the position vector of the i-th SPH particle. The incremental velocity would increase the rotation rate byω 0 ∆t if its shape did not change. However, deformation actually changes the moment of inertia of the body, which results in a further change of the rotation rate. We setω 0 = βω n , whereω n = 8.954 × 10 −10 rad/s 2 is the angular acceleration in the nominal case defining the dimensionless parameter β. Here, we define the spinup timescale of a body as the elapsed time of spinup from the rotation period of 3.5 h to 3.0 h. In the nominal case (β = 1), the spinup timescale is 9.3 × 10 4 s, which corresponds to 8.0 rotations. The nominal rotational acceleration of β = 1 is much faster than that due to the actual YORP effect for a kilometer-sized asteroid, which corresponds to β ∼ 10 −8 (Kaasalainen et al. 2007;Hergenrother et al. 2019). We present the results of rotational deformation through spinups with β = 1 in Section 3.1 and with much smaller β in Section 3.2. The spinup of the body eventually results in structural deformation and the ejection of surface masses. We accelerated only the SPH particles that comprise the main body, i.e., the largest clump of SPH particles, each of which has a distance from the nearest member of less than 1.5 times the smoothing length of a SPH particle. We identified clumps using a friends-of-friends algorithm (e.g., Huchra & Geller 1982).
Major deformation of the initial spherical body is followed by or accompanied by major mass ejection. In this study, we focused on the analysis of the shape of the main body formed through the major deformation event. To check the stability of the final shape of a deformed body, we stopped acceleration after 1% of the total mass was ejected and continued the simulations for more than 1.0 × 10 5 s, which is longer than a typical deformation timescale of the body.
The ejected mass at some stage of a simulation is defined as the total mass of SPH particles that do not belong to the main body at that time. We measured the ratio c/a of the minor to major axis of the main body using the top-down method (e.g., Capaccioni et al. 1984), where the two axis lengths a or c were measured as the distances between two parallel plates that bound the main body. The rotation period P of the main body is calculated as P = 2πI z /L z , where I z and L z are, respectively, the moment of inertia and the angular momentum of the main body around the z-axis. Note that the angle between the angular momentum vector and the z-axis is < 0.3 • . Although we do not do any corrections for the rotation axis of the body, the rotation axis is almost aligned with the z-axis. Thus, P calculated from I z and L z is a reliable measure for the rotation period.

Spinup with β = 1
In simulations with the nominal spinup rate β = 1, we set the initial rotation periods for φ d ≤ 40 • to be 5.5 h and those for φ d ≥ 50 • to be 3.5 h in order to save the computational time because in the latter cases deformation of the body never occurs during the periods when P > 3.5 h. Figure 1 shows side-view snapshots of the β = 1 spinup of the granular spherical bodies with the effective friction angles φ d = 40 • , 60 • , and 80 • . Figure 2 shows a time evolution of the axis ratio c/a of the main body obtained from the same simulations. Because of the difference of the initial rotation periods, time when deformation occurs for φ d = 40 • is ≈ 2.0 × 10 5 s later than those for φ d = 60 • and 80 • . For better comparison, elapsed time from 2.0 × 10 5 s after the start of the simulation is shown for φ d = 40 • in Figs. 1a-e and 2. Figure 1: Side-view snapshots of spinning-up bodies with effective friction angles φ d = 40 • (a-e), 60 • (f-j), and 80 • (k-o). Each panel includes elapsed time t and the instantaneous rotation period of the main body P . Note that t for φ d = 60 • and 80 • is the elapsed time from the start of the simulations, while that for φ d = 40 • is from 2.0 × 10 5 s after the start of the simulation.
All the deformed bodies in Fig. 1 have kept almost axisymmetric shapes, indicating that the spinup of a spherical body with β = 1 results in axisymmetric deformation. This behavior is quite different from a self-gravitating fluid body, for which rapid rotation induces a non-axisymmetric instability to form a series of Jacobi triaxial ellipsoids (Chandrasekhar 1969) 1 . The non-axisymmetric deformation into triaxial ellipsoids with small intermediate-to-major axis ratio < 0.5 occurs only for φ d ≤ 20 • . This suggests that relatively small values of φ d (∼ 30 • ) of granular bodies will prevent the deformation into triaxial ellipsoids. Note that non-axisymmetric sets of landslides occurring in slower spinup of bodies with φ d ≥ 70 • (see Section 3.2) belong to a different deformation mode from the internal deformation forming non-axisymmetric triaxial ellipsoids for φ d ≤ 20 • .
For φ d = 40 • , the axis ratio c/a decreases from unity ( Fig. 1a) to ≈ 0.5 ( Fig. 1d) during 1.0 × 10 4 s − 1.9 × 10 5 s ( is comparable to the spinup timescale ≈ 9.3 × 10 4 s, and the deformation of the body almost halts once we stop the spinup. We define this spinup-driven deformation mode as a quasi-static deformation. Further acceleration of the rotation during t > 1.9 × 10 5 s leads to mass ejection from the body's surface in the equatorial plane (Figs. 1e and 2). Note that the body spins down because of the increase in the moment of inertia due to the deformation of the body, although the angular momentum increases owing to the imposed spinup condition. In addition, the vertically asymmetric shape is formed (see Figs. 1ce), which is probably caused by initial asymmetry of SPH particle distribution along z-direction.
For φ d = 60 • , the initial spherical shape remains unchanged until t ≈ 7 × 10 4 s when a rapid deformation occurs (Fig. 2). The deformation timescale of ∼ 10 4 s is much shorter than the spinup timescale ≈ 9.3 × 10 4 s. Once the significant deformation occurs, the deformation is promoted even without the spin acceleration. We define this deformation mode as a dynamical deformation. After t ≈ 7 × 10 4 s, the decrease rate of the axis ratio slows until t ≈ 1.4 × 10 5 s. Significant mass ejection occurs at t ≈ 1.4 × 10 5 s (Fig. 2). The mass ejection under the slow-deformation phase seems similar to the case with φ d = 40 • .
For φ d = 80 • , the axis ratio suddenly decreases at t ≈ 9.0 × 10 4 s, which is a similar dynamical deformation to the case with φ d = 60 • (Fig. 2). However, the mass ejection from the surface in the equatorial plane is accompanied by the deformation of the main body (Figs. 1n and 2). Significant deformation begins at the time when P = 3.03 h (Fig. 1k). At the equator, the centrifugal forces are almost equal to the gravitational attraction of the body. Thus, the deformation pushes the material near the equator out of the surface of the body, which rapidly leads to mass ejection (Fig. 1n). The amount of the ejected mass is ∼ 10% of the main body. The angular momentum loss due to the mass ejection increases the rotation period, which stalls further deformation. The main body finally has axis ratio c/a = 0.76 (Fig. 2) and rotation period P = 4.8 h (Fig. 1o). The final shape seems similar to a top shape (Fig. 1o, see also the supplementary movie), but the axis ratio c/a = 0.76 is smaller than that of Ryugu c/a = 0. . The horizontal axis shows the ratio c/a of the minor to major axis. The right and left vertical axes show, respectively, the rotation period for ρ = 1.19 g/cm 3 and the angular speed ω scaled by √ πGρ, where G is the gravitational constant. (a) Red, blue, and green curves show the cases for the effective friction angle φ d = 20 • , 30 • , and 40 • , respectively. (b) Magenta, cyan, light green, and brown curves show the cases for φ d = 50 • , 60 • , 70 • , and 80 • , respectively. For comparison with the low φ d case shown on the panel (a), we plot evolutionary tracks of axis ratio and normalize spin rate obtained from two DEM spinup simulations: the black dotted curve is an HSDEM model from Walsh et al. (2012); the gray dotted curve is an SSDEM model from Sánchez & Scheeres (2012). The crosses and filled triangles on the solid curves show the initial states and the states at the times when the spinup is halted (when 1% of the total mass is ejected), respectively. For φ d = 20 • , the mass ejection occurs when P = 7.6 h and c/a = 0.16, so that the red triangle is out of the panel (a). Figure 3 shows the instantaneous spin rate ω of the main body as a function of axis ratio c/a obtained from our numerical simulations with various effective friction angles in comparison with the maximum equilibrium spin of oblate spheroids given by an analytic formula found by Holsapple (2001) and Sharma et al. (2009). For the cases with low effective friction angles φ d = 20 • , 30 • , and 40 • , ω initially increases while the bodies maintain their spherical shapes (c/a = 1) until ω approaches the maximum equilibrium angular speed. Then ω and c/a vary mostly according to the maximum equilibrium spin state. The mass ejection begins after significant deformation with c/a < 0.5. Figure 3a also shows the normalized spin rates and axis ratios of spinningup bodies simulated using the DEMs for cohesionless granular materials: the Hard-Sphere DEM (HSDEM) simulation in Walsh et al. (2012) (red solid curve in Fig. 9 of the paper) and the Soft-Sphere DEM (SSDEM) simulation in Sánchez & Scheeres (2012) (red solid curve in Fig. 26 of the paper). We only plot the data of the DEM simulations until major mass ejection events. The precise evaluation of the bulk friction is difficult for the DEM simulations, but the DEM results roughly correspond to the maximum equilibrium spin curves with φ d = 20 − 30 • . Although not only the methods but also the particle numbers, spinup rates, and spinup procedures are different from those of our simulations, the evolution paths of the DEM simulations are consistent with those of our simulations with φ d = 20 − 30 • . On the other hand, the axis ratio c/a of the HSDEM simulation at the epoch of mass ejection is largely different from those of the SSDEM and our simulations: c/a ≈ 0.8, 0.3, and 0.2 for the HSDEM, the SSDEM, and our simulation with φ d = 30 • , respectively. This suggests that SSDEM and SPH simulations provide more similar results than HSDEM simulations.
For φ d = 50 • , 60 • , 70 • , and 80 • , ω exceeds the maximum equilibrium angular speed without deformation followed by sudden dynamical deformation in which both c/a and ω decrease (Fig. 3b). These are quite different from the quasiequilibrium evolution path for lower φ d (Fig. 3a). The dynamical deformation timescale is much shorter than the spinup timescale (see Fig. 2). The sudden deformation allows the spin rate to be smaller than the maximum spin rate in equilibrium. Mass ejection occurs after the dynamical deformation for φ d = 50 • and 60 • , while during the dynamical deformation for φ d = 70 • and 80 • . Figure 4 shows the cross-sectional snapshots of the main bodies with different φ d during the dynamical deformation. The simulations with φ d = 50 • and 60 • (Fig. 4a and b) show that the deformation speeds at the polar regions and the equatorial regions are almost the same, and not only the surface but also the deep interior of the bodies deform. A typical flow pattern can be seen in the case with φ d = 60 • (Fig. 4b), where streamlines are similar to rectangular hyperbolae that tend to flatten the body. We define this deformation mode as internal deformation.
For φ d = 70 • and 80 • (Fig. 4c and d), deformation speeds in the surface layers from equator to mid-latitudes are faster than those in polar and internal regions, and deformation vectors in the surface layers point toward the equators. The typical flow pattern can be seen for φ d = 80 • (Fig. 4d), where deformation is negligible in the internal region and the deformation speeds are significant only at the surface layers. We define this deformation mode as surface landslides. The deformation occurs at the time when P = 3.07 h for φ d = 70 • and P = 3.03 h for φ d = 80 • . These rotation periods are very close to the rotational breakup period P lim = 2π/ (4/3)πGρ 0 = 3.027 h for a spherical body without tensile strength. Thus, the deformation directly causes significant mass ejection (see also Fig. 1k-m), and the sudden mass ejection from the equatorial surface triggers a set of surface landslides. The mass ejection simultaneously occurs at various longitudes, which induces an axisymmetric occurrence of multiple landslides at almost the same time. We refer to the axisymmetric and simultaneous occurrence of multiple landslides as an axisymmetric set of landslides.
In summary, three deformation modes according to the effective friction angle φ d are identified through our simulations with β = 1: quasi-static and internal deformation for φ d ≤ 40 • (Fig. 1a-e and Fig 3a); dynamical and internal deformation for 50 • ≤ φ d ≤ 60 • (Fig. 4a and b); and an axisymmetric set of surface landslides for φ d ≥ 70 • (Fig. 4c and d). Figure 5 shows the cross-sections of the main bodies at the end of the simulations. For φ d = 50 • and 60 • , oblate spheroidal shapes without equatorial ridges are formed ( Fig. 5a and b). The surface tilt angles relative to the z-direction are ∼ 30 • at low latitudes (latitudes of 0 • − 20 • ) and ∼ 60 • at mid latitudes (latitudes of 40 • − 60 • ). The large tilt angle differences of ∼ 30 • between the low-and mid-latitudes show that the main body has a near-spheroidal surface.
For φ d = 70 • and 80 • , the final shapes of the main bodies seem to have conical surfaces rather than spheroidal surfaces ( Fig. 5c and d). The difference of the surface tilt angles between the low-and mid-latitudes is ∼ 20 • for φ d = 70 • and ∼ 7 • for φ d = 80 • . Ryugu has differences of the surface tilt angles between low-and mid-latitudes 20 • at almost all longitudes (Watanabe et al. 2019). Thus the final shapes resulting from our simulations with φ d ≥ 70 • have conical surfaces similar to that of Ryugu. Such conical surfaces are most likely produced by landslides. Moreover, a top shape has a characteristic that the topographic distance from the center of the body has a minimum at a mid-latitude, and that characteristic is reproduced for the bodies formed through the simulations with φ d ≥ 70 • . Therefore, top shapes can be produced through an axisymmetric set of landslides caused by spinup with β = 1 and φ d ≥ 70 • (Figs.1o and 4d).
One cannot distinguish top shapes from spheroids by their axis ratios. Here, we adopt a quantitative method to distinguish top shapes from spherical shapes. We define the relative volume difference δ V as where V and V e are the volume of a body and an ellipsoid with the same triaxial dimensions as the body, respectively. We measured the axis length of the body by the top-down method (e.g., Capaccioni et al. 1984). Bodies with ellipsoidal or rounded shapes have δ V 1, while bodies with angular shapes such as top shapes have large δ V . Ryugu has δ V ≈ 0.3 (Watanabe et al. 2019).  and 60 • and increases with increasing φ d for φ d > 60 • . For φ d = 70 • and 80 • , δ V is close to or larger than that of Ryugu. This tendency is consistent with the results of our simulations: the body has spheroidal shapes for φ d ≤ 60 • but angular shapes for φ d ≥ 70 • (Fig. 5). Thus, we quantitatively confirmed that our simulations produce top shapes for φ d ≥ 70 • .
It should be noted that the top shape produced in the simulation with φ d = 80 • has c/a = 0.76, while almost all of the top shapes found so far have c/a > 0.85 (e.g., Ostro et al. 2006;Busch et al. 2011;Watanabe et al. 2019;Lauretta et al. 2019). Landslides mainly induce material transfer from high latitudes to equatorial regions (see Fig. 4d), and thus it is basically difficult for bodies to keep c/a of about unity through landslides. The discrepancy regarding c/a between the observed top shapes and that produced through our simulations should be considered in future work.

Spinup with β = 0.05
The deformation processes of the body with φ d ≤ 60 • under spinup with β = 0.05 are similar to those with β = 1, i.e., the quasi-static, axisymmetrical, and internal deformation occurs for φ d < 40 • , and the dynamical, axisymmetrical, and internal deformation occurs for 50 • ≤ φ d ≤ 60 • . The resultant shapes are axisymmetric. Thus, the spinup rate β does not seem to affect the deformation modes of the body with a low effective friction angle φ d ≤ 60 • . In contrast, simulations of the body with φ d ≥ 70 • under spinup with β = 0.05 result in a different deformation mode from those under β = 1. Figure  7 shows the snapshots of the simulation with φ d = 80 • and β = 0.05. The initial rotation period of the main body is set to be 3.25 h. At t ≈ 5.60 × 10 5 s, the rotation period becomes P ≈ 3.08 h and the mass ejection mainly starts from an equatorial region (Fig. 7b and e). At t ≈ 5.65 × 10 5 s, a significant amount of mass (∼ 10% of the total mass) is ejected, but it does not come from an equatorial region around the +y-direction (Fig. 7c and f). Therefore, non-axisymmetric mass ejection occurs. The rotation period becomes P ≈ 4.5 h due to angular momentum loss through the mass ejection.  The landslides occur at a longitude of 45 • E (Fig. 8a), where the straight surface in the meridional cross-section appears in the northern hemisphere. However, only a small portion of mass is ejected from longitudes around 180 • E (Fig. 8b), and no mass ejection occurs at longitudes around 270 • E (Fig. 8c). The body's surface profiles in the cross-sections at longitudes around 180 • E and 270 • E maintain almost circular profiles.
Although the main body initially has an almost uniform density, small fluctuations in the SPH particle distribution cause slight locational variations in surface roughness and local slopes, which cause different stability levels against local landslides at different surface positions. In the case of the spinup with β = 1, major landslides occur when P = 3.03 h (Fig. 1k). Under such a fast rotation, landslides are induced everywhere by strong centrifugal forces regardless of local surface slopes. However, in the case of slower spinup with β = 0.05, major landslides occur at an earlier time when P = 3.08 h (Fig. 7). Owing to the weaker centrifugal forces at that time, landslides seem to selectively occur on surfaces with locally steeper slopes, allowing surfaces with locally shallower slopes to remain mostly intact. This may result in localized landslides or a nonaxisymmetric set of landslides, forming a non-axisymmetric shape. The angular momentum of the main body is lost due to the landslides in regions with locally steeper slopes, which leads to a decrease of the rotation rate. Thus, further landslides in the regions with locally shallower slopes do not occur.

Spinup of a body with an equatorial ridge
As shown in Section 3.2, even small fluctuations of SPH particle distribution cause non-axisymmetric distribution of landslides. Thus, we expect that a visible surface topographic high affects occurrence of landslides. We examined the spinup of a body with an equatorial ridge as an example of spinup of a body with surface topographic high. For simplicity, we put a half torus with a height of H ridge on a spherical body to represent an equatorial ridge (see Fig. 9a). We set the effective friction angle to φ d = 80 • and the initial rotation period to 4.0 h. We found that the simulation of spinup with H ridge = 25 m and β = 1 results in a non-axisymmetric set of landslides. Figure 9b and c show cross-sections of the body at two different longitudes. At a longitude of 90 • E (Fig. 9b), landslides occur and a straight surface in this meridional cross section that resembles those seen in cross-sections of top shapes is produced. However, at the longitude of 270 • E (Fig. 9c), no landslides occur and the resultant cross-section remains spherical. Thus, the spinup of a body with the surface topographic highs results in a non-axisymmetric set of landslides.
In contrast, the simulation with H ridge = 25 m and faster spinup rate β = 2 results in an axisymmetric set of landslides, producing a top shape. In Sections 3.1 and 3.2, we see that the simulation with H ridge = 0 m and β = 1 results in an axisymmetric set of landslides, but that with H ridge = 0 m and β = 0.05 results in a non-axisymmetric set of landslides. These results suggest that critical spinup rates that are required for an axisymmetric set of landslides become slower for bodies with less surface roughness.

Discussion
Here, we discuss implications for the formation of real top-shaped asteroids in light of the results of our simulations using the SPH method.

Cohesion and large effective friction angle
Our simulations demonstrated that the effective friction angle φ d of the constituent material mainly determines the deformation modes of the rotating bodies. We found that spinup of a spherical body with effective friction angles φ d ≤ 40 • induces quasi-static and internal deformation, whereas that with φ d ≥ 50 • results in dynamical deformation. In the dynamical regime, internal deformation occurs in the cases with φ d ≤ 60 • , whereas surface landslides occur in φ d ≥ 70 • . These transitions of the deformation modes seem to be irrespective of the spinup rate β according to our simulations with 0.05 ≤ β ≤ 1. However, for φ d ≥ 70 • , axisymmetricity of the resultant shape depends on β; landslides occur almost axisymmetrically through fast spinup (β = 1), while local landslides with a non-axisymmetric distribution occur through slow spinup (β = 0.05).
A body with φ d ≤ 60 • deforms into an oblate spheroid through internal deformation, but never into a top shape that has conical surfaces extending from the equator to the northern and southern mid-latitudes. In contrast, a body with φ d ≥ 70 • can deform into an axisymmetric top shape through an axisymmetric set of landslides. Thus, we suggest that the formation of top shapes requires large effective friction angles of φ d ≥ 70 • . As we discussed in Section 2.2, we interpreted φ d used in our simulations as the effective friction angle that mimics an aspect of cohesion. Here, we show the validity of using large effective friction angles instead of introducing the cohesion.
Recall, the shear strength of granular material Y d is expressed as where φ d,real is the friction angle, p is the confining pressure, and c coh is the cohesion. Thus, the effective friction angle φ d , i.e., the arctangent value of the ratio of the shear strength to the confining pressure, is expressed as We conducted numerical simulations with the same setup shown in Section 3.1 but with the shear strength of Eq. (4). We set φ d,real = 40 • and c coh = 30, 50, and 100 Pa, which, respectively, correspond to the effective friction angles of φ d = 55 • , 61 • , and 71 • for the central pressure p = (2/3)πGρ 2 0 R 2 ≈ 50 Pa of an asteroid with radius R = 500 m and density ρ 0 = 1.19 g/cm 3 . Note that the central pressure is the maximum pressure in an asteroid and gives the minimum estimated value of effective friction angles. We found that the resultant shapes obtained from the simulations with c coh = 30, 50, and 100 Pa are an oblate spheroid (Fig. 10a), a slightly flattened top shape (Fig. 10b), and a top shape (Fig. 10c), respectively. The resultant shapes and axis ratios c/a with c coh = 30, 50, and 100 Pa are similar to those with effective friction angles of φ d = 60 • (Fig. 1j and Fig. 4b), 70 • (Fig. 4c), and 80 • (Fig. 1o and Fig. 4d), respectively. For example, c/a of the body in Fig. 1o and Fig. 10c is 0.76 and 0.75, respectively. Thus, the simulations with the shear strength of Eq. (4) agree well with those with the shear strength of Y d = p tan(φ d ) and corresponding effective friction angles φ d .
While we find the agreement, obviously Eq. (4) is not completely the same as Y d = p tan(φ d ) even with corresponding effective friction angles φ d . For example, the confining pressures vary in post-yield flow of grain material, which may be minor effect for the spinup deformation of kilometer sized asteroids, and the shear strength with Eq. (4) somewhat deviates from that with Y d = p tan(φ d ). The surface profile of the body shown in Fig. 10a is more bumpy than that with the effective friction angle φ d = 60 • (Fig. 5), of which difference is quantitatively represented by the difference of δ V ; δ V for the former is 0.16 while that for the latter is 0.071. This may be caused by the difference of the constitutive laws. Moreover, we ignored the direct effect of cohesion, that is, tensile forces of co-hesive materials, which may affect the results of our simulations. Investigation of the detailed effect of cohesion is planned for future work.
The cohesion of a subsurface layer of Ryugu c coh = 140 − 670 Pa (Arakawa et al. 2020) results in the effective friction angle φ d ≈ 75 − 86 • with a typical confining pressure for a kilometer sized body p ≈ 50 Pa and a friction angle φ d,real = 40 • . Thus, large effective friction angles of φ d ≥ 70 • are plausible for 1 km-sized asteroids. In contrast, a larger object has larger confining pressure and smaller effective friction angle. The central pressure of a 10 km-sized asteroid is p c ≈ 5 kPa, which leads to an effective friction angle φ d ≈ 44 • even with c coh = 670 Pa. Thus, our simulations imply that it would be difficult to form a top-shaped body with a diameter larger than 10 km through spinup, which is consistent with the fact that all top-shaped asteroids found so far are smaller than 10 km.

Formation of a top shape through fast spinup
Our results show that top shapes are formed by the fast spinup (β ≥ 1) of spherical bodies with effective friction angles ≥ 70 • (Fig. 1k-o). However, the slow spinup (β = 0.05) of spherical bodies with ≥ 70 • induces a nonaxisymmetric set of landslides and results in non-axisymmetric shapes (Fig. 7). Moreover, surface roughness or topography suppresses top-shape formation. The fast spinup (β = 1) of a body with a 25 m high equatorial ridge results in a non-axisymmetric shape (Fig. 9).
We discuss the formation of axisymmetric top shapes through YORP spinup (Section 4.2.1). We also discuss fast spinup through the reaccumulation of fragments after catastrophic collisions (Section 4.2.2), although applying our spinup simulations to spinup through the accretion requires the assumption that mass increase and shape changes are negligible through the accretion, which is discussed in Section 4.2.2. In addition, we discuss possibilities of other spinup mechanisms (Section 4.2.3). Here, instead of spinup rate β, we use spinup timescale t su , which is defined as the elapsed time to spin up the body from rotation period P = 3.5 h to P = 3.0 h: t su ≈ 9.3 × 10 4 β −1 s.

YORP effect
The YORP effect for a kilometer-sized asteroid causes spinup with a timescale t su 100 kyr ≈ 3 × 10 12 s (e.g., Kaasalainen et al. 2007;Hergenrother et al. 2019). We showed that a non-axisymmetric set of landslides tends to occur under slower spinups even in the case with β = 0.05 or t su ≈ 1.9 × 10 6 s. This suggests that a non-axisymmetric set of landslides and the formation of a nonaxisymmetric shape are likely results from the YORP spinup of an asteroid.
We showed that a body with less surface roughness experiences an axisymmetric set of landslides even for slower spinup. The kilometer-sized spherical bodies with 25,000 SPH particles initially set in our simulations are estimated to have a surface roughness of ∼ 10 m due to discretization of the bodies using SPH particles. The numerical roughness is larger than typical surface roughness of asteroids (e.g., Sugita et al. 2019), and thus the spinup timescale of asteroids required for an axisymmetric set of landslides may be longer than what we obtained in the simulations. We extrapolated the results obtained in Sections 3.2 and 3.3 through a power-law fitting and found that an axisymmetric set of landslides for bodies with 1 cm surface roughness requires t su ∼ 10 7 s; the required spinup timescale is still much shorter than the timescale of the YORP effect.
However, we do not rule out the possibility that the YORP spinup will produce a top-shaped asteroid. The YORP effect may spin up an asteroid even after a major landslide event if deformation and topographical changes caused by the major landslides insignificantly alter the direction of the net YORP torque exerted on the body, which eventually causes further landslides. The place where a local landslide has occurred has a locally conical surface and has been stabilized against further landslides. Thus, the subsequent landslides due to further YORP spinup may occur at the longitudes where significant local landslides have not occurred yet, and eventually multiple occurrences of local landslides at different longitudes would produce an almost axisymmetric top shape. Since we only investigated a shape formed until a single major landslide event, we do not rule out the possibility of the formation of top shapes through multiple cycles of the YORP spinup and local landslides.

Reaccumulation of fragments
The reaccumulation phase after catastrophic disruption of the parent body of a rubble-pile asteroid by a collision is yet another practical situation for fast spinup of the asteroid through accretion Michel et al. 2020). Note that Michel et al. (2020) focus on the possibility of top-shape formation through deformation due to the accretion process, but we focus on top-shape formation through spinup due to accretion. The catastrophic collision produces a sheet-like structure composed of fragments, and dense parts of the sheet gravitationally collapse into small cores of remnants. Then, the remaining fragments accumulate on the remnants to supply angular momentum. Note that the supplied angular momentum originates from the small part of the sheet with not random but specific velocity field around it, so that the accretion of the fragments typically causes continuous spinup of the remnant. Moreover, the timescale of the reaccumulation is ∼ (Gρ 0 ) −1/2 ≈ 4 × 10 3 s, which is shorter than the spinup timescale t su required for the formation of an axisymmetric top shape t su ∼ 10 5 s. The fast spinup caused through the reaccumulation may form a top shape.
Accretion of fragments supplies not only angular momentum but also mass onto the core. Moreover, accretion may induce surface mass flow, which may affect how subsequent landslides occur. These effects were not represented in our simulations; we just accelerated the spin of the initial body. However, the mass of fragments required for spinup may be negligible compared to the core mass. We estimated angular momentum increase due to the accretion as M f Rv esc , where M f is the total mass of fragments accreting on the core, R = 500 m is the radius of the core, and v esc ≈ 40 cm/s is the escape speed from the core. Note that M f Rv esc roughly explains angular momentum increase through accretion in the simulations of the catastrophic collisions . The amount of fragments required to accelerate the core's rotation period from 3.5 h to 3.0 h is estimated to be only 4% of the core mass. Thus, our simulations of the spinup possibly mimic the spinup due to the reaccumulation of fragments. The possibility remains to be studied in future works.

Other mechanisms with short spinup timescales
Acceleration rates until the rotation periods of bodies reach the critical rotation period ∼ 3 h are not related to how landslides occur. Thus, the following formation procedure of axisymmetric top shapes is possible: the YORP effect accelerates the rotation periods of bodies up to ∼ 3 h and then other mechanisms with short spinup timescales further accelerate the bodies and cause an axisymmetric set of landslides. In this case, small amounts of angular momentum supply due to the fast spinup mechanisms are sufficient to cause an axisymmetric set of landslides.
A non-destructive impact of a small asteroid supplies angular momentum and may cause fast spinup. A prograde impact (i.e., angular momentum vector provided by an impactor is parallel to that of a rotating body) tends to cause tentative spinup while the impactor contacts with the rotating body. Since the impact may erode the body and take away angular momentum to some extent, the process is not straightforward and thus this phenomena should be investigated via simulations in future works.
An asteroid passing around the Roche radius of a planet acquires angular momentum via tidal torque (Hyodo et al. 2016), while an encounter much closer than the Roche radius results in elongation and disruption of the asteroid (e.g., Walsh & Richardson 2006). A moderate encounter may result in the spinup of an asteroid with keeping its shape. All the top-shaped asteroids found so far are near-Earth asteroids, so that a close encounter with the Earth is one possible spinup mechanism for the formation of top-shaped asteroids.

Subsequent evolution of ejected materials via landslides
Our simulations showed that an axisymmetric set of landslides that produces a top-shaped asteroid ejects fragments with mass ∼ 0.1M b (Section 3.1), where M b is the mass of the central body. The ejected fragments produce satellites with ∼ 10% of the total fragments (Hyodo et al. 2015). The mass of the satellites ∼ 1% of host asteroids is consistent with those of satellites around observed topshaped asteroids (Ostro et al. 2006;Brozović et al. 2011;Becker et al. 2015). However, about half of observed top-shaped asteroids do not have satellites (e.g., Ryugu or Bennu). This suggests that there is a mechanism to remove satellites from the top-shaped asteroids within their lifetimes.
An orbit of a satellite is mainly altered by tidal interaction with a host asteroid (e.g., Goldreich & Sari 2009), the binary YORP (BYORP) effect (e.g., Ćuk & Burns 2005;McMahon & Scheeres 2010a,b), and a close encounter with a planet (e.g., . The YORP effect is probably the dominant satellite-evolution mechanism for a binary near-Earth asteroids with binary separation ∼ the Roche radius of a host asteroid, although the orbital evolution through each mechanism highly depends on various physical parameters. The BYORP effect is an analog of the YORP effect for a binary system with a satellite spinning synchronously with its orbital motion, and torque induced by the absorption of sunlight by and the thermal re-emission from the irregularly shaped satellite changes the orbit of the satellite. The timescale of the secular evolution of the semimajor axis a due to the BYORP effect is estimated as where ρ 0 is the bulk density of the binary asteroids, ω d = (4/3)πGρ 0 is the critical spin frequency, R h is the radius of the host asteroid, M s /M h is the mass ratio of the satellite and the host asteroid, the heliocentric orbit factor H is 4 × 10 −5 g cm −1 s −2 at 1 au (Scheeres 2007), and the dimensionless parameter B s that depends on the satellite's shape is estimated to be ∼ 0.01 (McMahon & Scheeres 2010a;Walsh & Jacobson 2015). If the BYORP effect causes the outward migration of the satellite, the satellite is ejected from the binary system within ∼ 10 5 yr, which is much shorter than the dynamical lifetime of near-Earth asteroids ∼ 10 Myr (Gladman et al. 2000). Thus the BYORP effect can eject the satellite during its stay in the near-Earth space, which may explain the lack of satellites for several top-shaped asteroids.

Summary
Rubble-pile asteroids Ryugu and Bennu visited by the spacecraft Hayabusa2 and OSIRIS-REx, respectively, have a similar spinning-top shape with an axisymmetric equatorial ridge. Many top-shaped asteroids rotate near their breakup periods of ∼ 3 hours. The axisymmetricity and present-day high spin rates suggest that these top-shaped asteroids were formed through some spinup processes. However, how a rubble-pile body with a specific friction angle deforms under various spinup rates is still under debate.
Given the radius and bulk density, we numerically simulated the spinup of a spherical rubble-pile body with different effective friction angles φ d using the SPH method for granular material. The high φ d mimics the effect of shear cohesion, and we confirmed that the simulation results with φ d are similar to those with the real cohesive strength law (Eq. (4)) with corresponding values of cohesion. The non-dimensional acceleration rate of the rotation was defined as β =ω 0 /ω n , where a normalizing acceleration rateω n = 8.954 × 10 −10 rad/s 2 corresponds to a decrease in rotation period from 3.5 h to 3.0 h within 9.3×10 4 s (≈ 8.0 rotations). We also defined the spinup timescale as t su ≈ 9.3 × 10 4 β −1 s. We mainly investigated spinups of a spherical body with an acceleration rate β = 1 (t su ≈ 9.3 × 10 4 s) and effective friction angles φ d = 20 − 80 • . Contrary to the fluid cases (φ d = 0 • ), our simulations show that a moderate effective friction angle (φ d ≥ 30 • ) will suppress the deformation into a triaxial ellipsoid. In the cases with φ d ≤ 40 • , quasi-static and internal deformation occurred, forming oblate spheroids (Fig. 1a-e). The simulations with 50 • ≤ φ d ≤ 60 • resulted in dynamical internal deformation and formed oblate spheroids ( Fig. 1f-j, Fig. 4a, b, and Fig. 5a, b), indicating that internal deformation will not lead to topshape formation. In contrast, spinup of spherical bodies with φ d ≥ 70 • induced an axisymmetric set of surface landslides, producing axisymmetric top shapes ( Fig. 1k-o, Fig. 4c, d, and Fig. 5c, d).
We also investigated spinup with a slower acceleration rate β = 0.05 (t su ≈ 1.9 × 10 6 s) for comparison. We found that the change of β probably does not affect the principal deformation modes found in the nominal case (β = 1): quasi-static and internal deformation for φ d ≤ 40 • ; dynamical and internal deformation for 50 • ≤ φ d ≤ 60 • ; and a set of surface landslides for φ d ≥ 70 • . However, the slower β affects axisymmetricity of the distribution of landslides for φ d ≥ 70 • . The simulation with φ d = 80 • and β = 0.05 (t su ≈ 1.9×10 6 s) showed that a non-axisymmetric shape is formed through a set of local landslides.
We showed that the formation of axisymmetric top shapes through spinup requires an effective friction angle φ d ≥ 70 • , which is much larger than those of cohesionless granular materials. However, estimated small values of the cohesion for kilometer-sized asteroids (c coh ∼ 100 Pa) make effective friction angles larger than 70 • . Thus, the cohesion of asteroids plays an important role for the formation of axisymmetric top shapes.
In addition, our simulations showed that the formation of axisymmetric top shapes prefers fast spinup with timescales shorter than ∼ 10 6 s. The timescale of the YORP spinup ( 100 kyr) is much longer than the spinup timescales required for the formation of axisymmetric top shapes, suggesting a difficulty to produce axisymmetric top shapes via deformation through YORP spinup itself. On the other hand, catastrophic collisions of parent asteroids produce many small remnants and accretion of fragments onto these remnants may spins them up within the reaccumulation timescale ∼ 4×10 3 s. Thus, catastrophic collisions are one possible origin of axisymmetric top shapes. Besides reaccumulation of fragments after catastrophic collisions, coalescence of fragments that are ejected through the YORP spinup, a non-destructive and prograde impact of an asteroid, and tidal torque caused by a close encounter with a planet may spin up asteroids with short timescales and may lead to the formation of top shapes. Our simulations showed that the formation of top-shaped asteroids is accompanied by surface mass ejection, which leads to the formation of satellites. About half of observed top-shaped asteroids do not have satellite, and this can be explained by the ejection of satellites such as outward migration of satellites' orbits due to the BYORP effect.