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

The numerical code utilized in this study is based on a version of the SPH code described in Sugiura et al., (2018, 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, 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:

We mainly investigated the rotational deformation of a uniform and spherical body with a radius of $500\,{\rm 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 $\approx 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(\phi_{d,{\rm real}})p+c_{{\rm coh}}$, where $\phi_{d,{\rm real}}$ and $c_{{\rm coh}}$ are the “real” friction angle and the cohesive strength, respectively. For simplicity, we adopted the relation in a form $Y_{d}=\tan(\phi_{d})p$ with the “effective” friction angle $\phi_{d}$ by interpreting that the relation $\phi_{d}=\tan^{-1}(Y_{d}/p)$ would be valid for cohesive granular materials. The effective friction angle $\phi_{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_{{\rm coh}}=670\,{\rm 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\,{\rm m}$ and mean density of $\rho_{0}=1.19\,{\rm g/cm^{3}}$ is only $p_{{\rm c}}=(2/3)\pi G\rho_{0}^{2}R^{2}\approx 50\,{\rm Pa}$, so that the effective friction angle is up to $\phi_{d}\approx\tan^{-1}(c_{{\rm coh}}/p_{{\rm c}})\approx 86^{\circ}$. To mimic large effective friction angles due to cohesion, we varied $\phi_{d}$ from $20^{\circ}$ to $80^{\circ}$. This treatment allows us to understand the detailed dependence of deformation modes on one parameter $\phi_{d}$. The validity of using such large effective friction angles instead of introducing the cohesion $c_{{\rm 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_{{\rm min}}$ required for keeping the initial configuration gets shorter for larger $\phi_{d}$ (Holsapple, 2001). We set the initial rotation period of the body in each simulation with an effective friction angle $\phi_{d}$ to be longer than $P_{{\rm min}}$ for the corresponding $\phi_{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 $\Delta\bm{v}_{i}$ in the rotational direction to the $i$-th SPH particle comprising the body at each numerical step as

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\times 10^{5}\,{\rm 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\pi 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^{\circ}.$ 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.

In simulations with the nominal spinup rate $\beta=1$, we set the initial rotation periods for $\phi_{d}\leq 40^{\circ}$ to be $5.5\,{\rm h}$ and those for $\phi_{d}\geq 50^{\circ}$ to be $3.5\,{\rm 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\,{\rm h}$. Figure 1 shows side-view snapshots of the $\beta=1$ spinup of the granular spherical bodies with the effective friction angles $\phi_{d}=40^{\circ}$, $60^{\circ}$, and $80^{\circ}$. 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 $\phi_{d}=40^{\circ}$ is $\approx 2.0\times 10^{5}\,{\rm s}$ later than those for $\phi_{d}=60^{\circ}$ and $80^{\circ}$. For better comparison, elapsed time from $2.0\times 10^{5}\,{\rm s}$ after the start of the simulation is shown for $\phi_{d}=40^{\circ}$ in Figs. 1a-e and 2.

All the deformed bodies in Fig. 1 have kept almost axisymmetric shapes, indicating that the spinup of a spherical body with $\beta=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)¹¹1The stability of granular asteroids has been investigated by Sharma, (2013) in detail. Although the work is relevant to our simulations, direct comparison between the work and our simulations is difficult because setups are different. Especially, our present study simulates bodies in which global plastic deformation hardly occurs, so that analysis of Sharma, (2013) on plastic modulus is not expected to affect our present study. In future work, we are planning to perform simulations that use the same setup as the work and discuss application for the stability of granular asteroids.. The non-axisymmetric deformation into triaxial ellipsoids with small intermediate-to-major axis ratio $<0.5$ occurs only for $\phi_{d}\leq 20^{\circ}$. This suggests that relatively small values of $\phi_{d}$ ($\sim 30^{\circ}$) of granular bodies will prevent the deformation into triaxial ellipsoids. Note that non-axisymmetric sets of landslides occurring in slower spinup of bodies with $\phi_{d}\geq 70^{\circ}$ (see Section 3.2) belong to a different deformation mode from the internal deformation forming non-axisymmetric triaxial ellipsoids for $\phi_{d}\leq 20^{\circ}$.

For $\phi_{d}=40^{\circ}$, the axis ratio $c/a$ decreases from unity (Fig. 1a) to $\approx 0.5$ (Fig. 1d) during $1.0\times 10^{4}\,{\rm s}-1.9\times 10^{5}\,{\rm s}$ (Fig. 2). Timescale of the deformation is comparable to the spinup timescale $\approx 9.3\times 10^{4}\,{\rm 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\times 10^{5}\,{\rm 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. 1c-e), which is probably caused by initial asymmetry of SPH particle distribution along $z$-direction.

For $\phi_{d}=60^{\circ}$, the initial spherical shape remains unchanged until $t\approx 7\times 10^{4}\,{\rm s}$ when a rapid deformation occurs (Fig. 2). The deformation timescale of $\sim 10^{4}\,{\rm s}$ is much shorter than the spinup timescale $\approx 9.3\times 10^{4}\,{\rm 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\approx 7\times 10^{4}\,{\rm s}$, the decrease rate of the axis ratio slows until $t\approx 1.4\times 10^{5}\,{\rm s}$. Significant mass ejection occurs at $t\approx 1.4\times 10^{5}\,{\rm s}$ (Fig. 2). The mass ejection under the slow-deformation phase seems similar to the case with $\phi_{d}=40^{\circ}$.

For $\phi_{d}=80^{\circ}$, the axis ratio suddenly decreases at $t\approx 9.0\times 10^{4}\,{\rm s}$, which is a similar dynamical deformation to the case with $\phi_{d}=60^{\circ}$ (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\,{\rm 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 $\sim 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\,{\rm 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.87$ and Bennu $c/a=0.88$ (Watanabe et al., 2019; Lauretta et al., 2019).

Figure 3 shows the instantaneous spin rate $\omega$ 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 $\phi_{d}=20^{\circ}$, $30^{\circ}$, and $40^{\circ}$, $\omega$ initially increases while the bodies maintain their spherical shapes ($c/a=1$) until $\omega$ approaches the maximum equilibrium angular speed. Then $\omega$ 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 spinning-up 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 $\phi_{d}=20-30^{\circ}$. 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 $\phi_{d}=20-30^{\circ}$. 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\approx 0.8$, $0.3$, and $0.2$ for the HSDEM, the SSDEM, and our simulation with $\phi_{d}=30^{\circ}$, respectively. This suggests that SSDEM and SPH simulations provide more similar results than HSDEM simulations.

For $\phi_{d}=50^{\circ}$, $60^{\circ}$, $70^{\circ}$, and $80^{\circ}$, $\omega$ exceeds the maximum equilibrium angular speed without deformation followed by sudden dynamical deformation in which both $c/a$ and $\omega$ decrease (Fig. 3b). These are quite different from the quasi-equilibrium evolution path for lower $\phi_{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 $\phi_{d}=50^{\circ}$ and $60^{\circ}$, while during the dynamical deformation for $\phi_{d}=70^{\circ}$ and $80^{\circ}$.

Figure 4 shows the cross-sectional snapshots of the main bodies with different $\phi_{d}$ during the dynamical deformation. The simulations with $\phi_{d}=50^{\circ}$ and $60^{\circ}$ (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 $\phi_{d}=60^{\circ}$ (Fig. 4b), where streamlines are similar to rectangular hyperbolae that tend to flatten the body. We define this deformation mode as internal deformation.

For $\phi_{d}=70^{\circ}$ and $80^{\circ}$ (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 $\phi_{d}=80^{\circ}$ (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\,{\rm h}$ for $\phi_{d}=70^{\circ}$ and $P=3.03\,{\rm h}$ for $\phi_{d}=80^{\circ}$. These rotation periods are very close to the rotational breakup period $P_{{\rm lim}}=2\pi/\sqrt{(4/3)\pi G\rho_{0}}=3.027\,{\rm 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 $\phi_{d}$ are identified through our simulations with $\beta=1$: quasi-static and internal deformation for $\phi_{d}\leq 40^{\circ}$ (Fig. 1a-e and Fig 3a); dynamical and internal deformation for $50^{\circ}\leq\phi_{d}\leq 60^{\circ}$ (Fig. 4a and b); and an axisymmetric set of surface landslides for $\phi_{d}\geq 70^{\circ}$ (Fig. 4c and d).

Figure 5 shows the cross-sections of the main bodies at the end of the simulations. For $\phi_{d}=50^{\circ}$ and $60^{\circ}$, oblate spheroidal shapes without equatorial ridges are formed (Fig. 5a and b). The surface tilt angles relative to the $z$-direction are $\sim 30^{\circ}$ at low latitudes (latitudes of $0^{\circ}-20^{\circ}$) and $\sim 60^{\circ}$ at mid latitudes (latitudes of $40^{\circ}-60^{\circ}$). The large tilt angle differences of $\sim 30^{\circ}$ between the low- and mid-latitudes show that the main body has a near-spheroidal surface.

For $\phi_{d}=70^{\circ}$ and $80^{\circ}$, 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 $\sim 20^{\circ}$ for $\phi_{d}=70^{\circ}$ and $\sim 7^{\circ}$ for $\phi_{d}=80^{\circ}$. Ryugu has differences of the surface tilt angles between low- and mid-latitudes $\lesssim 20^{\circ}$ at almost all longitudes (Watanabe et al., 2019). Thus the final shapes resulting from our simulations with $\phi_{d}\geq 70^{\circ}$ 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 $\phi_{d}\geq 70^{\circ}$. Therefore, top shapes can be produced through an axisymmetric set of landslides caused by spinup with $\beta=1$ and $\phi_{d}\geq 70^{\circ}$ (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 $\delta_{V}$ as

Figure 6 shows $\delta_{V}$ of the main body resulting from the spinup of the spherical body with $\beta=1$ and $0.25$. For both $\beta$ values, $\delta_{V}$ is $\approx 0.1$ for $\phi_{d}=40^{\circ}$, $50^{\circ}$, and $60^{\circ}$ and increases with increasing $\phi_{d}$ for $\phi_{d}>60^{\circ}$. For $\phi_{d}=70^{\circ}$ and $80^{\circ}$, $\delta_{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 $\phi_{d}\leq 60^{\circ}$ but angular shapes for $\phi_{d}\geq 70^{\circ}$ (Fig. 5). Thus, we quantitatively confirmed that our simulations produce top shapes for $\phi_{d}\geq 70^{\circ}$.

It should be noted that the top shape produced in the simulation with $\phi_{d}=80^{\circ}$ 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.

The deformation processes of the body with $\phi_{d}\leq 60^{\circ}$ under spinup with $\beta=0.05$ are similar to those with $\beta=1$, i.e., the quasi-static, axisymmetrical, and internal deformation occurs for $\phi_{d}<40^{\circ}$, and the dynamical, axisymmetrical, and internal deformation occurs for $50^{\circ}\leq\phi_{d}\leq 60^{\circ}$. The resultant shapes are axisymmetric. Thus, the spinup rate $\beta$ does not seem to affect the deformation modes of the body with a low effective friction angle $\phi_{d}\leq 60^{\circ}$.

In contrast, simulations of the body with $\phi_{d}\geq 70^{\circ}$ under spinup with $\beta=0.05$ result in a different deformation mode from those under $\beta=1$. Figure 7 shows the snapshots of the simulation with $\phi_{d}=80^{\circ}$ and $\beta=0.05$. The initial rotation period of the main body is set to be $3.25\,{\rm h}$. At $t\approx 5.60\times 10^{5}\,{\rm s}$, the rotation period becomes $P\approx 3.08\,{\rm h}$ and the mass ejection mainly starts from an equatorial region (Fig. 7b and e). At $t\approx 5.65\times 10^{5}\,{\rm s}$, a significant amount of mass ($\sim 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\approx 4.5\,{\rm h}$ due to angular momentum loss through the mass ejection.

Figure 8 shows the cross-sections of the main body at different longitudes. The landslides occur at a longitude of $45^{\circ}$ 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^{\circ}$ E (Fig. 8b), and no mass ejection occurs at longitudes around $270^{\circ}$ E (Fig. 8c). The body’s surface profiles in the cross-sections at longitudes around $180^{\circ}$ E and $270^{\circ}$ 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 $\beta=1$, major landslides occur when $P=3.03\,{\rm 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 $\beta=0.05$, major landslides occur at an earlier time when $P=3.08\,{\rm 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 non-axisymmetric 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.

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_{{\rm ridge}}$ on a spherical body to represent an equatorial ridge (see Fig. 9a). We set the effective friction angle to $\phi_{d}=80^{\circ}$ and the initial rotation period to $4.0\,{\rm h}$.

We found that the simulation of spinup with $H_{{\rm ridge}}=25\,{\rm m}$ and $\beta=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^{\circ}$ 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^{\circ}$ 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_{{\rm ridge}}=25\,{\rm m}$ and faster spinup rate $\beta=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_{{\rm ridge}}=0\,{\rm m}$ and $\beta=1$ results in an axisymmetric set of landslides, but that with $H_{{\rm ridge}}=0\,{\rm m}$ and $\beta=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.

Our simulations demonstrated that the effective friction angle $\phi_{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 $\phi_{d}\leq 40^{\circ}$ induces quasi-static and internal deformation, whereas that with $\phi_{d}\geq 50^{\circ}$ results in dynamical deformation. In the dynamical regime, internal deformation occurs in the cases with $\phi_{d}\leq 60^{\circ}$, whereas surface landslides occur in $\phi_{d}\geq 70^{\circ}$. These transitions of the deformation modes seem to be irrespective of the spinup rate $\beta$ according to our simulations with $0.05\leq\beta\leq 1$. However, for $\phi_{d}\geq 70^{\circ}$, axisymmetricity of the resultant shape depends on $\beta$; landslides occur almost axisymmetrically through fast spinup ($\beta=1$), while local landslides with a non-axisymmetric distribution occur through slow spinup ($\beta=0.05$).

A body with $\phi_{d}\leq 60^{\circ}$ 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 $\phi_{d}\geq 70^{\circ}$ 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 $\phi_{d}\geq 70^{\circ}$. As we discussed in Section 2.2, we interpreted $\phi_{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

We conducted numerical simulations with the same setup shown in Section 3.1 but with the shear strength of Eq. (4). We set $\phi_{d,{\rm real}}=40^{\circ}$ and $c_{{\rm coh}}=30,\,50$, and $100\,{\rm Pa}$, which, respectively, correspond to the effective friction angles of $\phi_{d}=55^{\circ},\,61^{\circ}$, and $71^{\circ}$ for the central pressure $p=(2/3)\pi G\rho_{0}^{2}R^{2}\approx 50\,{\rm Pa}$ of an asteroid with radius $R=500\,{\rm m}$ and density $\rho_{0}=1.19\,{\rm 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_{{\rm coh}}=30,\,50$, and $100\,{\rm 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_{{\rm coh}}=30,50$, and $100\,{\rm Pa}$ are similar to those with effective friction angles of $\phi_{d}=60^{\circ}$ (Fig. 1j and Fig. 4b), $70^{\circ}$ (Fig. 4c), and $80^{\circ}$ (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(\phi_{d})$ and corresponding effective friction angles $\phi_{d}$.

While we find the agreement, obviously Eq. (4) is not completely the same as $Y_{d}=p\tan(\phi_{d})$ even with corresponding effective friction angles $\phi_{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(\phi_{d})$. The surface profile of the body shown in Fig. 10a is more bumpy than that with the effective friction angle $\phi_{d}=60^{\circ}$ (Fig. 5), of which difference is quantitatively represented by the difference of $\delta_{V}$; $\delta_{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 cohesive 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_{{\rm coh}}=140-670\,{\rm Pa}$ (Arakawa et al., 2020) results in the effective friction angle $\phi_{d}\approx 75-86^{\circ}$ with a typical confining pressure for a kilometer sized body $p\approx 50\,{\rm Pa}$ and a friction angle $\phi_{d,{\rm real}}=40^{\circ}$. Thus, large effective friction angles of $\phi_{d}\geq 70^{\circ}$ are plausible for $1\,{\rm km}$-sized asteroids. In contrast, a larger object has larger confining pressure and smaller effective friction angle. The central pressure of a $10\,{\rm km}$-sized asteroid is $p_{{\rm c}}\approx 5\,{\rm kPa}$, which leads to an effective friction angle $\phi_{d}\approx 44^{\circ}$ even with $c_{{\rm coh}}=670\,{\rm Pa}$. Thus, our simulations imply that it would be difficult to form a top-shaped body with a diameter larger than $10\,{\rm km}$ through spinup, which is consistent with the fact that all top-shaped asteroids found so far are smaller than $10\,{\rm km}$.

Our results show that top shapes are formed by the fast spinup ($\beta\geq 1$) of spherical bodies with effective friction angles $\geq 70^{\circ}$ (Fig. 1k-o). However, the slow spinup ($\beta=0.05$) of spherical bodies with $\geq 70^{\circ}$ induces a non-axisymmetric set of landslides and results in non-axisymmetric shapes (Fig. 7). Moreover, surface roughness or topography suppresses top-shape formation. The fast spinup ($\beta=1$) of a body with a $25\,{\rm 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 $\beta$, we use spinup timescale $t_{{\rm su}}$, which is defined as the elapsed time to spin up the body from rotation period $P=3.5\,{\rm h}$ to $P=3.0\,{\rm h}$: $t_{{\rm su}}\approx 9.3\times 10^{4}\beta^{-1}\,{\rm s}$.

Our simulations showed that an axisymmetric set of landslides that produces a top-shaped asteroid ejects fragments with mass $\sim 0.1M_{{\rm b}}$ (Section 3.1), where $M_{{\rm b}}$ is the mass of the central body. The ejected fragments produce satellites with $\sim 10\%$ of the total fragments (Hyodo et al., 2015). The mass of the satellites $\sim 1\%$ of host asteroids is consistent with those of satellites around observed top-shaped 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., Walsh & Richardson, 2008). The YORP effect is probably the dominant satellite-evolution mechanism for a binary near-Earth asteroids with binary separation $\sim$ 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 (McMahon & Scheeres, 2010a; Walsh & Jacobson, 2015)