INTERNAL STRUCTURE OF ASTEROIDS HAVING SURFACE SHEDDING DUE TO ROTATIONAL INSTABILITY

The Yarkovsky–O'Keefe–Radzievskii–Paddack effect, so-called the YORP effect, may play a crucial role in providing asymmetric asteroids with rotational torque due to the solar radiation pressure (Rubincam 2000). This implies that asteroids with several 100 m in diameter may be gradually able to spin up to their spin limits (Rubincam 2000). Observational surveys have shown that there is a spin period threshold of ∼2.3 hr for asteroids larger than a few 100 m in size (Pravec et al. 2007). Recent observations detected spin-up/down of asteroids due to this effect (Lowry et al. 2007; Kaasalainen et al. 2007; Taylor et al. 2007; Ďurech et al. 2008). These reports imply that spin-up processes could cause asteroid structures and shapes to evolve over their life time.

Some recently observed disruption events of active asteroids are considered to have resulted from rotational instability; however, there are other proposed mechanisms for these events as well. Jewitt et al. (2010) and Snodgrass et al. (2010) independently reported the disruption event of asteroid P/2010 A2. Jewitt et al. (2010) concluded that it was likely that the formation of the debris tail resulted from rotational instability, while Snodgrass et al. (2010) noted that a collisional impact is a likely event driving this event. Jewitt et al. (2014) described the recent breakup event of P/2013 R3, suggesting that this event came from rotational instability. Hsieh et al. (2004, 2010) reported that the dust tail of asteroid 133P/Elst-Pizarro, spinning with a spin period of 3.471 hr, may come from seasonal activity and rotational instability. Sheppard & Trujillo (2015) also observed a thin dust tail of asteroid (62412) 2000 SY178 whose spin period is 3.33 hr from 2014 March 28 to May 2. Although the activities of these bodies could be driven by sublimation of water ice that might be buried in the bodies, the slow particles in their dust tails also imply that rotational instability may also play a role in mass ejection from the surface. Jewitt et al. (2013, 2015) conducted observations of P/2013 P5 (PANSTARRS) with the Hubble Space Telescope to investigate its episodic dust ejection, implying that this event might result from rotational instability. Hainaut et al. (2014) and Hainaut & Snodgrass (2014) proposed an alternative scenario of this event in which the dust tail might be generated by the soft contact in a binary system gently rubbing each other.

Over the last decade, there have been many important analytical (Holsapple 2001, 2004, 2007, 2010) and numerical (Walsh et al. 2008, 2012; Sánchez & Scheeres 2012) investigations of failure conditions for asteroids. Of specific relevance for this study, it has been found that for a homogeneous strength distribution in an asteroid, the internal structure is more sensitive to failure than the surface region (Hirabayashi & Scheeres 2015). Another active research area is an understanding of mass ejection from a fast rotating body due to centrifugal forces. A popular explanation for the origin of mass ejection is that surface materials flow down to the equator, and if the rotation is fast enough, they fly away. Guibout & Scheeres (2003) investigated the surface stability due to the shape and spin period of a uniformly rotating, self-gravitating ellipsoid and found that based on the body configuration, the surface materials may accumulate at different ends. Using a hard-sphere discrete element method (HSDEM), Walsh et al. (2008, 2012) showed the formation of a binary system due to mass ejection from a cohesionless spheroid. By considering a constant angle of repose, Minton (2008) and Harris et al. (2009) numerically obtained equilibrium shapes of a uniformly rotating symmetric body. Scheeres (2015) analytically investigated granular surface flows on rotating asteroids and their equilibrium shapes, demonstrating the formation of unique surface slopes of oblate shapes such as 1999 KW4 Alpha and 2008 EV5. Also, Hirabayashi & Scheeres (2014) compared surface mass ejection with structural failure due to rapid rotation to determine a possible failure mode of (216) Kleopatra.

However, the mechanism behind surface mass ejection is still poorly understood, and thus we will try to shed some light on it in this paper. Developing a two-bulk-density-layer model, Hirabayashi (2014) showed that a high-density core makes the whole body structurally strong. Assuming a homogenous structure, Hirabayashi & Scheeres (2015) conducted plastic finite element analysis to explain that the failure mode of asteroid (29075) 1950 DA is plastic deformation of the internal region. These studies imply that in order for the surface region of an asteroid to fail first at an elevated spin rate, its interior must be stronger than a surface layer. Based on these studies, we propose one way to create surface shedding from a rubble pile asteroid. Here, we introduce a simplified model to discuss it. Our rubble pile asteroid is modeled as a self-gravitating spheroid composed of two cohesive strength layers, a spherical shell and an internal core. The bulk density and the friction angles are assumed to be constant to simplify our discussion. We note that cohesion may depend on porosity (Lambe & Whitman 1969), but we will leave consideration of these quantities for a future study.

In the following discussion, we will use the term "surface shedding." So, before we do, we will try to give it a proper definition. Hirabayashi & Scheeres (2014) introduced this term to explain the mode that particles just resting on the surface are shed due to strong centrifugal forces. However, this definition is still ambiguous. First, if particles on the surface mechanically interact with other particles by cohesion, it is difficult to determine the lofting condition by simple balance between gravitational and centrifugal forces. Second, if cohesion is uniformly distributed over the whole volume, the internal structure deforms before the force balance point, i.e., the dynamical equilibrium point, reaches the surface (Holsapple 2010; Hirabayashi & Scheeres 2015). To define this mode clearly, we describe "surface shedding" as a mode in which while the internal structure is below the yield condition, the surface region fails structurally and particles in this region are shed due to strong centrifugal forces. We hypothesize that this mode may play a role in the disruption events of 133P/Elst-Pizarro (Hsieh et al. 2004, 2010), P/2013 P5 (Jewitt et al. 2013, 2015), and 62412 (Sheppard & Trujillo 2015).

2.1. Internal Core Model

We suppose that the test body is a uniformly rotating, self-gravitating spheroid that has a concentric, spherical core under the surface shell, later known as the internal core (Figure 1). It is assumed that these layers have different cohesion, while their bulk density and friction angle are the same over the whole volume. Surface shedding is induced not by the failure of the internal core but by that of the surface shell. If the test body is spinning fast enough, surface particles are shed from the test body; otherwise, they just move over the surface.

Zoom In Zoom Out Reset image size

The bulk density, the total radius, and the gravitational constant are denoted as ρ, R, and G, respectively. We normalize lengths, body forces, spin rates, and stress tensors by R, $\pi \rho {GR}$, $\sqrt{\pi \rho G}$, and $\pi {\rho }^{2}{{GR}}^{2}$, respectively. With this normalization, the sphere radius is denoted as 1 and the radius of the internal core is defined as R_b, so the thickness of the surface shell is given as $1-{R}_{b}$. The right plot in Figure 1 shows the spherical coordinate system $(r,\theta ,\phi )$, but we also use the cartesian coordinate system $(x,y,z)$ as well. The spin axis is assumed to be constant along the z axis. We consider that at a given spin rate, if the surface shell already fails structurally, the cohesion between the particles there is broken and they can move freely. For this case, the particle is lofted from the surface due to the centrifugal force at a spin rate of $\sqrt{4/3}\sim 1.15$, while it energetically escapes to infinity at that of $\sqrt{8/3}\sim 1.63$.

2.2. Failure Conditions of the Two Layers

In the analytical model developed here, the behavior of materials is assumed to be elastic-perfectly plastic and to follow an associated flow rule. The yield condition of a material is assumed to be characterized by the Mohr–Coulomb yield criterion, a pressure-shear dependent yield criterion, which is given as

$$\begin{eqnarray}&&({\sigma }_{1}-{\sigma }_{3})\mathrm{sec}\phi +({\sigma }_{1}+{\sigma }_{3})\mathrm{tan}\phi \leqslant 2Y,\end{eqnarray} \tag{ 1 }$$

where ${\sigma }_{i}\;(i=1,2,3)$ is the principal stress component, ϕ is a friction angle, and Y is cohesion. Since we consider a sphere, the stress state is always ${\sigma }_{1}={\sigma }_{2}\geqslant {\sigma }_{3}$.

The yield conditions of the surface shell and the internal core are calculated by using an elastic stress solution. Solving the elastic equations, including the equilibrium equations, the constitutive equations (Hooke's law), the strain–displacement equations, and proper boundary conditions (zero-traction condition), we obtain the solution as (Dobrovolskis 1982)

$$\begin{eqnarray}&&{\sigma }_{{xx}}={k}_{1}(1-{r}^{2})-{k}_{13}{x}^{2}-{k}_{14}{z}^{2},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{\sigma }_{{yy}}={k}_{5}(1-{r}^{2})-{k}_{13}{y}^{2}-{k}_{15}{z}^{2},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\sigma }_{{zz}}={k}_{9}(1-{r}^{2})-{k}_{14}{x}^{2}-{k}_{15}{y}^{2},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{\sigma }_{{xy}}={k}_{13}{xy},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\sigma }_{{xz}}={k}_{14}{xz},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\sigma }_{{yz}}={k}_{15}{yz},\end{eqnarray} \tag{ 7 }$$

where ${r}^{2}={x}^{2}+{y}^{2}+{z}^{2}$ and the coefficients, ks, are defined as follows:

$$\begin{eqnarray}&&{k}_{1}=\displaystyle \frac{2{b}_{x}(-12-\nu +5{\nu }^{2})+{b}_{z}(3-6\nu -5{\nu }^{2})}{10(-1+\nu )(7+5\nu )},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{k}_{9}=\displaystyle \frac{2{b}_{x}(3-6\nu -5{\nu }^{2})+{b}_{z}(-27+4\nu +15{\nu }^{2})}{10(-1+\nu )(7+5\nu )},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{k}_{13}=-\displaystyle \frac{2{b}_{x}(-4+3\nu +5{\nu }^{2})+{b}_{z}(1+3\nu )}{5(-1+\nu )(7+5\nu )},\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{k}_{14}=-\displaystyle \frac{{b}_{x}(-3+6\nu +5{\nu }^{2})+{b}_{z}(-4+3\nu +5{\nu }^{2})}{5(-1+\nu )(7+5\nu )}.\end{eqnarray} \tag{ 11 }$$

Also, ${k}_{5}={k}_{1}$ and ${k}_{14}={k}_{15}$. ν is Poisson's ratio. The coefficients of the body force components along the x, y, and z axes, denoted as b_x, b_y, and b_z, respectively, are given as

$$\begin{eqnarray}&&{b}_{x}={b}_{y}={\omega }^{2}-\displaystyle \frac{4}{3},\;{b}_{z}=-\displaystyle \frac{4}{3}.\end{eqnarray} \tag{ 12 }$$

To determine lower and upper bound conditions for the actual failure, we apply the upper and lower bound theorems in limit analysis (Chen & Han 1988; Chen & Liu 1990). At a lower bound condition, the target volume should not fail structurally. At an upper bound condition, on the other hand, the volume should fail. Thus, the actual failure condition should always be between these conditions.

The lower bound theorem describes that when (1) an equilibrium stress solution is found and (2) there is a small element whose stress reaches its yield condition, the body may not widely fail. This means that the yield condition of a small element does not guarantee further propagation of failure regions. Therefore, we say that this condition is below the actual failure condition. Practically, this condition is given as the yield condition of the stress at the most sensitive element. Because the test body is spherical, such an element always appears on the equatorial plane.

For the lower bound condition of the surface shell, the boundary between the surface shell and the internal core becomes the most sensitive to failure. Considering y = z = 0, we obtain the principal stress components at this location as

$$\begin{eqnarray}&&{\sigma }_{e1}={k}_{1}(1-{r}^{2}),\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&{\sigma }_{e2}={k}_{1}(1-{r}^{2})-{k}_{13}{r}^{2},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\sigma }_{e3}={k}_{9}(1-{r}^{2})-{k}_{13}{r}^{2}.\end{eqnarray} \tag{ 15 }$$

It is not always satisfied that ${\sigma }_{e1}\gt {\sigma }_{e2}\gt {\sigma }_{e3}$. Thus, to have the correct order, we calculate ${\sigma }_{e{\rm{s}}1}=\mathrm{max}({\sigma }_{e1},{\sigma }_{e2},{\sigma }_{e3})$ and ${\sigma }_{e{\rm{s}}3}=\mathrm{min}({\sigma }_{e1},{\sigma }_{e2},{\sigma }_{e3})$. On the other hand, for the lower bound condition of the internal core, since the center is always the most sensitive location, we obtain the principal stress components as

$$\begin{eqnarray}&&{\sigma }_{e{\rm{c}}1}={k}_{1},\;{\sigma }_{e{\rm{c}}2}={k}_{1},\;{\sigma }_{e{\rm{c}}3}={k}_{9}.\end{eqnarray} \tag{ 16 }$$

The lower bound condition provides the highest approximated value of cohesion for the actual failure condition, later known as the highest cohesion. Based on Equation (1), this is expressed as

$$\begin{eqnarray}&&{Y}_{{ek}}=\displaystyle \frac{1}{2}({\sigma }_{{ek}1}-{\stackrel{}{\sigma }}_{{ek}3})\mathrm{sec}\phi +\displaystyle \frac{1}{2}({\stackrel{}{\sigma }}_{{ek}1}+{\stackrel{}{\sigma }}_{{ek}3})\mathrm{tan}\phi ,\end{eqnarray} \tag{ 17 }$$

where k is either s or c. k = s indicates cohesion of the surface shell, while k = c provides that of the internal core.

For the upper bound theorem, on the other hand, using the theory by Chen & Han (1988) and by Chen & Liu (1990), Holsapple (2008) developed a practical technique. His theorem guarantees that on the assumption that the velocity field is linear, for a given volume, the yield condition of the averaged stress over it is equal to an upper bound condition for the actual failure. The averaged principal stress over a given volume V_k is given as

$$\begin{eqnarray}&&{\bar{\sigma }}_{{jk}}=\displaystyle \frac{1}{{V}_{k}}{\displaystyle \int }_{{V}_{k}}{\sigma }_{j}{{dV}}_{k}.\end{eqnarray} \tag{ 18 }$$

Here, we average the stresses over the surface shell and the internal core. Using Equation (18), we calculated the averaged stresses of the surface shell as

$$\begin{eqnarray}&&{\bar{\sigma }}_{{\rm{s}}1}={\bar{\sigma }}_{{\rm{s}}2}=-{k}_{1}\displaystyle \frac{{R}_{b}^{3}-{R}_{b}^{5}}{1-{R}_{b}^{3}}+\displaystyle \frac{{b}_{x}}{5}\displaystyle \frac{1-{R}_{b}^{5}}{1-{R}_{b}^{3}},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{\bar{\sigma }}_{{\rm{s}}3}=-{k}_{9}\displaystyle \frac{{R}_{b}^{3}-{R}_{b}^{5}}{1-{R}_{b}^{3}}+\displaystyle \frac{{b}_{z}}{5}\displaystyle \frac{1-{R}_{b}^{5}}{1-{R}_{b}^{3}},\end{eqnarray} \tag{ 20 }$$

and those of the internal core as

$$\begin{eqnarray}&&{\bar{\sigma }}_{{\rm{c}}1}={\bar{\sigma }}_{{\rm{c}}2}={k}_{1}\left(1-{R}_{b}^{2}\right)+\displaystyle \frac{{b}_{x}{R}_{b}^{2}}{5},\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{\bar{\sigma }}_{{\rm{c}}3}={k}_{9}\left(1-{R}_{b}^{2}\right)+\displaystyle \frac{{b}_{z}{R}_{b}^{2}}{5}.\end{eqnarray} \tag{ 22 }$$

The details of these equations are in the Appendix. To verify these stress components, we show that they recover the averaged stresses over the whole volume obtained by Holsapple (2007). If ${R}_{b}=0$, the averaged stress over the surface shell times $\pi {\rho }^{2}{{GR}}^{2}$ is identical to Equation (4.12) with $a=b=c=R$ in Holsapple (2007). If ${R}_{b}=1$, Equation (4.12) is recovered again because this case is reduced to aggregates with a homogeneous interior. This consistency allows us to compare the upper bound conditions of the surface shell and the internal core. The yield conditions of these stress components provide the lowest estimated value of the actual cohesion, later known as the lowest cohesion. Similar to Equation (17), this is described as

$$\begin{eqnarray}&&{Y}_{k}=\displaystyle \frac{1}{2}({\bar{\sigma }}_{k1}-{\bar{\sigma }}_{k3})\mathrm{sec}\phi +\displaystyle \frac{1}{2}({\bar{\sigma }}_{k1}+{\bar{\sigma }}_{k3})\mathrm{tan}\phi .\end{eqnarray} \tag{ 23 }$$

In Section 4.1, we discuss the upper and lower bound conditions for the surface shell and the upper bound condition for the internal core. For the internal core, since the surface shedding mode is described by a core that does not reach the yield, the lower bound condition for the internal core does not provide constraints on it. Therefore, we do not consider this condition in the present problem.

The simulation program that is used for this research applies a SSDEM to simulate a self-gravitating granular aggregate (Cundall 1971; Cundall & Hart 1992; Sánchez & Scheeres 2011). The particles, modeled as spheres that follow a predetermined size distribution, interact through a soft-repulsive potential when in contact. This method considers that two particles are in contact when they overlap. When this happens, normal and tangential contact forces are calculated (Herrmann & Luding 1998).

The calculation of the normal forces between colliding particles is modeled by a hertzian spring and a dashpot. The elastic force is modeled as

$$\begin{eqnarray}&&{{\boldsymbol{f}}}_{{\rm{e}}}={k}_{{\rm{n}}}{\xi }^{3/2}{\boldsymbol{}}\hat{n},\end{eqnarray} \tag{ 24 }$$

the damping force as

$$\begin{eqnarray}&&{{\boldsymbol{f}}}_{{\rm{d}}}=-{\gamma }_{{\rm{n}}}\dot{\xi }{\boldsymbol{}}\hat{n},\end{eqnarray} \tag{ 25 }$$

and the cohesive force between the particles is calculated as

$$\begin{eqnarray}&&{{\boldsymbol{f}}}_{{\rm{c}}}=-2\pi \displaystyle \frac{{r}_{1}^{2}{r}_{2}^{2}}{{r}_{1}^{2}+{r}_{2}^{2}}{\sigma }_{{yy}}{\hat{{\boldsymbol{r}}}}_{12}\end{eqnarray} \tag{ 26 }$$

where r₁ and r₂ are the radii of the two particles in contact, ${\sigma }_{{yy}}$ is the tensile strength of this contact, which is given by a cohesive matrix formed by the (non-simulated) interstitial regolith (Sánchez & Scheeres 2014), and ${\hat{{\boldsymbol{r}}}}_{12}$ is the branch vector between the centers of these two particles. Then, the total normal force is calculated as ${{\boldsymbol{f}}}_{{\rm{n}}}={{\boldsymbol{f}}}_{{\rm{e}}}+{{\boldsymbol{f}}}_{{\rm{c}}}+{{\boldsymbol{f}}}_{{\rm{d}}}$. In these equations, k_n is the elastic constant, ξ is the overlap of the particles, ${\gamma }_{{\rm{n}}}$ is the damping constant (related to the dashpot), $\dot{\xi }$ is the rate of deformation, and ${\boldsymbol{}}\hat{n}$ is the normal vector on the surfaces of the colliding particles. As the simulation only deals with spheres, ${\hat{{\boldsymbol{r}}}}_{12}$ has the direction as $\hat{{\boldsymbol{n}}}$. This dashpot models the energy dissipation that occurs during a real collision.

The tangential component of the contact force models surface friction statically, and dynamically. This is calculated by placing a linear spring attached to both particles at the contact point at the beginning of the collision (Herrmann & Luding 1998; Silbert et al. 2001) and by producing a restoring frictional force ${{\boldsymbol{f}}}_{{\rm{t}}}$. The magnitude of the elongation of this tangential spring is truncated in order to satisfy the local Coulomb yield criterion $| {{\boldsymbol{f}}}_{{\rm{t}}}| \leqslant \mu | {{\boldsymbol{f}}}_{{\rm{n}}}|$.

Rolling friction (Ai et al. 2011) has also been implemented in order to mimic the behavior of aggregates formed by non-spherical grains. That is, particles are subjected to a torque that opposes the relative rotation of any two particles in contact. This torque, similar to surface–surface friction, is implemented as linearly dependent on the relative angular displacement of any two particles in contact and has a limiting value of

$$\begin{eqnarray}&&{M}_{r}^{m}={\mu }_{r}{R}_{r}| {{\boldsymbol{f}}}_{{\rm{n}}}| ,\end{eqnarray} \tag{ 27 }$$

where ${\mu }_{r}$ is the coefficient of rolling resistance, ${R}_{r}={r}_{1}{r}_{2}/({r}_{1}+{r}_{2})$ is the rolling radius, and r₁ and r₂ are the radii of the two particles in contact. This allows our simulations to reach friction angles of up to $\approx 35^\circ$ as evaluated by the Druker–Prager yield criterion (Sánchez & Scheeres 2012). This value is typical of geological aggregates, though friction angles of ∼40° are not rare. This implementation of rolling friction is similar to that of surface friction, but is instead related to the relative angular displacement.

To obtain the initial configuration of the aggregates, we leave the particles (initially cohesionless and frictionless) to coalesce only under the influence of their mutual gravitational interactions. Then, the particles are encapsulated inside a perfectly solid sphere that shrinks to a predetermined size. This size has been pre-determined from previous simulations so that the particles are forced into a close-to-perfect spherical shape. When the set-up procedure is finished, friction and cohesive forces are applied and the simulation is ready to start. Once the simulation starts, the aggregate is spun up in small, discrete increments of a normalized spin rate of $5.5\times {10}^{-3}$ every 3000 s around the spin axis and the disruption process is observed. This time step proved to be long enough to allow the aggregates to safely reshape before the next spin-up event. The aggregates rotate in the −z direction to facilitate the observation of the reshaping and disruption process through the graphical interface of the simulation code.

For this research, to observe the effect of an internal core in a self-gravitating aggregate, we have chosen to have cores whose radii are 0.5, 0.6, 0.7, 0.8, and 0.9. The particles forming the core are subjected to cohesive forces that are five times as strong as those felt by the particles of the shell. All the aggregates used in the simulations are formed by 3000 perfectly spherical particles with normalized radii between 0.025 and 0.035. The normalized parameters of the model are as follow: ${k}_{{\rm{n}}}=6.1\times {10}^{7}$, ${k}_{{\rm{t}}}=2.1\times {10}^{6}$, ${\gamma }_{{\rm{n}}}=67$, $\mu =0.5$, and ${\mu }_{r}=0.8$.

4.1. Failure Condition of the Internal Structure

With the analytical model developed above, we now investigate the effect of the internal core on the surface failure and the surface shedding condition. The friction angle is fixed at 35°, a mean value of typical friction angles of geological materials (Lambe & Whitman 1969). For a choice of Poisson's ratio, we take account of the earlier studies by Holsapple (2008) and by Hirabayashi & Scheeres (2014). They confirmed that Poisson's ratio do not critically change the failure condition.⁵ Here, to take into account compressibility of a soil material, we define Poisson's ratio as 0.25. With these fixed values, we derive the upper and lower bound conditions for the surface shell and the upper bound condition for the internal core.

Figure 2 shows the upper and lower conditions for the normalized cohesion as a function of the normalized core radius. The spin rate is constant on each line. We have also plotted the 1.15 spin rate, at which a particle resting on the surface can be gently lofted, and the 1.63 spin rate, at which a particle escapes from the gravity of the test body. The former and latter spin rates are given by the dotted and dashed lines, respectively. Figure 2(a) gives the failure condition of the surface shell. The bold lines define the lowest cohesion at a given spin rate, while the narrow lines provide the highest cohesion at that spin rate. The actual failure spin rate should be between these two lines. At a constant spin rate, as the core radius increases, the highest and lowest cohesion decrease. On the other hand, though the surface shell should fail structurally, the internal core should be below the yield (see Figure 2(b)). The bold lines describe the lowest cohesion of the internal core at a given spin rate, while the narrower lines shows the highest cohesion of the surface shell at that spin rate. These results show that at a constant spin rate, the lowest cohesion that the internal core needs to remain stable is higher than the cohesion that the surface shell requires to be able to fail. This means that if the cohesion of the internal core is not high enough, it will tend to fail before the surface. In a homogeneous body, failure should start therefore at its center and propagate toward the surface. This scenario corresponds to the possible failure mode of asteroid 1950 DA (Hirabayashi & Scheeres 2015).

Zoom In Zoom Out Reset image size

Figures 2(a) and (b) imply that if surface shedding is observed, details of the process could provide information about internal structure. Given a constant spin period, we could give constraints on cohesive strength. Since the actual cohesion should be enclosed by the lowest and highest cohesion curves, it should be possible to obtain the range of cohesion based on a possible range of the core radius. In general, the core radius is not well determined; however, it is still possible to give a constraint on it. For example, let us consider a test body that experiences surface shedding at a spin rate of 1.3. For the surface shell (Figure 2(a)), taking the minimum value of the lowest cohesion (0.04 at a core radius of 1.0) and the maximum value of the highest cohesion (0.3 at that of 0), we roughly obtain the range of cohesion of the surface shell as 0.03–0.3. On the other hand, for the internal core (Figure 2(b)), we roughly give that of the internal core as $\gt 0.14$.

4.2. Dynamics of Failure Mode

This section investigates dynamics of failure modes by changing the size of the internal core. The initial conditions of the friction angle and the porosity are given based on experimental and observational evidences. The friction angle considered in the numerical simulation is 35° which is consistent with the analytical model. The porosity is calculated to be 34%, similar to that of the lunar regolith observed by the Apollo 15 mission (Mitchell et al. 1972). This is calculated by dividing the total mass of the particles into the volume obtained from semi-axes of a dynamically equivalent ellipsoid with the same moments of inertia as the aggregate.

The normalized cohesion of the internal core and the surface shell is 0.5 and 0.1, respectively.⁶ To obtain these values, we conduct the following numerical settings; given the tensile strength between the particles of the shell, ${\sigma }_{\mathrm{ts}}$, the tensile strength between the particles in the core is described as ${\sigma }_{\mathrm{tc}}=5{\sigma }_{\mathrm{ts}}$ and that between the particles in the interface is obtained as ${\sigma }_{\mathrm{ti}}=1.5{\sigma }_{\mathrm{ts}}$. These values of cohesion were chosen as previous simulations of homogeneous aggregates showed their critical spin rates and disruption modes to be very different; a 0.1 cohesion homogeneous aggregate fails at 1.13 (Figure 5) and a 0.5 cohesion aggregate breaks at 1.44 (Figure 8). These figures are explained later. Also, as seen in Figure 2, from the analytical model, known are the spin rates at which homogeneous aggregates would deform and disrupt, and so the ratio of cohesive forces warrants that the core would not fail before the surface.

It is important to note that the yield condition and irreversible flow in the SSDEM are consistent with the Drucker–Prager yield criterion and an associated flow rule, respectively (Sánchez & Scheeres 2012), while the analytical mode considers the Mohr–Coulomb yield criterion and an associated flow rule. In the present study, we use these different yield criteria to show how they are consistent with each other in the present problem.

We first compare the failure conditions obtained from the SSDEM and the analytical models. Figure 3 shows the spin-up profile for the cases of normalized core radii of 0, 0.5, 0.6, 0.7, 0.8 and 0.9, respectively. The peak spin rates for these cases are 1.13, 1.16, 1.18, 1.21, 1.27 and 1.33, respectively. After experiencing their peak spin rates, the test body cannot reach this spin rate again. This comes from the fact that the shape has already changed irreversibly, causing the test body to be more susceptible to failure. An increase in spin rate induces further deformation with the consequent slow-down in spin rate to conserve angular momentum. Figure 4 shows the critical spin rate at a constant surface cohesion of 0.10 as a function of the core radius. The solid lines are obtained by the analytical model, while the dotted line is given by the SSDEM. The numerical simulation results are located between the lower and upper bound conditions calculated by the analytical model, showing the consistency of these two approaches

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next, we observe the deformation of the test body that occurs right after the rapid spin-rate drop to investigate the initial process of its failure mode. Figure 5 shows the deformation mode of the test body with a normalized core radius of 0. Figure 5(a) shows the three-dimensional depiction of the aggregate after failure. The color shows the latitude of the test body. Figures 5(b) through (d) show the deformation vectors over the thin slices that include the center of mass and which are perpendicular to the x through z axes, respectively. Figures 5(b) and (c) (views from +x and +y) show how the poles of the aggregate push inwards toward the center, helping in the formation of an equatorial bulge. Figure 5(d) (view from +z) shows how the aggregate deforms outwards from the center. This deformation mode corresponds to the deformation of 1950 DA (Hirabayashi & Scheeres 2015) and could lead to its catastrophic disruption of the internal core. A possible disruption mode for this case will be discussed in Section 5.

Zoom In Zoom Out Reset image size

As seen in Figures 6 and 7, on the other hand, the effect of the internal core on the deformation is obvious. The deformation of the internal core is minimal compared to the homogeneous (zero-core) case, while the surface shell fails structurally.⁷ For the case of a normalized radius of 0.5 (Figure 6), after the surface shell fails, the particles there flow over the surface to reach the equatorial region. It is found that the failure region is not global, but local. For the case of a normalized radius of 0.9 (Figure 7), the failure region becomes smaller than that for the case of a normalized radius of 0.5 and mainly occurs near the equatorial region.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Now we analyze Figures 5 through 7 more closely to understand the effect of the internal core. These cases show the following three characteristics. First, the flow of particles is asymmetric with respect to the rotation axis and is symmetric with respect to the the equatorial plane, causing a ridge-like shape. Second, the region with the greatest deformation is flat and is diametrically opposed to the region that is least deformed and that is still very round. Third, points on the equatorial region between these two extremes show intermediate degrees of deformation/stretching (oblateness). These characteristics become less and less pronounced as the core grows in size. For a radius of 0.9, there is no deformation toward an oblate shape and the particles of the equator are simply ejected from the main body. The deformation of the test body happens only due to the failure of the surface shell, so the thicker it is, the greater the deformation. This implies that for the case of a normalized core radius of 0, the deformation is global (Sánchez & Scheeres 2012).

Our simulations, in no case show a deformation toward a symmetric oblate body and instead favor a wedge-like shape. As homogeneous as we would like our aggregates to be, by construction they present internal weak points where the body (as any real aggregate) can fail more easily. This is the reason for local nature of the region of greater deformation. Once the aggregate deforms in one region, its internal stress is diminished as the angular velocity has to decrease to conserve angular momentum. In this paper we have only focused on the initial phase on deformation so that the dynamics of surface shedding can be understood. However, future efforts will be directed toward understanding of its stable shape.

The simulations also show the dynamics of surface shedding. In Figure 4, the lower dashed line represents the condition at which the particles are just lofted from the surface, i.e., 1.15. If the spin rate is faster than this line and the surface shell already fails, surface shedding should occur. The upper dashed line represents the condition at which particles escape to infinity, i.e., 1.63. According to this figure, if the core radius is 0, particles cannot be shed. For the cases of core radii of 0.5 and larger, particles are lofted. However, since all the cases are below the upper dashed line, particles cannot reach escape velocity. They would either stay in orbit around the primary body or come back to its surface. Potentially, an absolute increase in the cohesion of the shell and interior would enable shed particles to immediately escape from the primary body.

In this paper, we have proposed a possible structure for a rubble pile asteroid that experiences surface shedding. It is important to consider what we know about failure modes from ground observations. Holsapple (2010) found that if the internal structure is uniformly distributed, the original body completely deforms before surface shedding. This predicts that a homogeneous spheroid would never have surface shedding, but could experience a catastrophic breakup at high enough cohesion. Sánchez & Scheeres (2014a, 2014b) simulated the zero core case for 560,000 s, which is 10,000 s longer than the initial failure time (∼460,000 s). The cohesion was fixed at 0.5, which is the same cohesion used for the internal core in the previous sections. Figure 8 shows a snapshot of the failure mode at 560,000 s. The test body catastrophically breaks up into multi-components. This mode is clearly distinguished from surface shedding. The surface shedding mode originates from the failure of the surface shell, causing granular flow.

Zoom In Zoom Out Reset image size

Several asteroids recently observed are now considered to have undergone surface shedding. Asteroid P/2013 P5 (PANSTARRS) was reported to experience its episodic dust ejection (Jewitt et al. 2013, 2015). This asteroid is orbiting near the inner edge of the asteroid belt in the vicinity of the Flora family of S-type asteroids. It implies that sublimation of water ice is unlikely to occur on this asteroid (Jewitt et al. 2013). The size of dust particles was up to a few millimeters, negligibly smaller than the nucleus, and the ejection velocity was less than 1 m s⁻¹ (Jewitt et al. 2015). From these observational results, the episodic event of this asteroid is currently considered to come from surface shedding. Asteroids 133P/Elst-Pizarro and (62412) 2000 SY178 have relatively rapid spin periods of 3.47 and 3.33 hr, respectively, and had experienced mass ejection. Although the mass ejection might have been driven by the seasonal activity, i.e., sublimation of water ice, rotational instability is also considered to help their activity (Hsieh et al. 2004, 2010; Sheppard & Trujillo 2015). Some remarkable observations related to their activity include: (1) the size of ejecta is small, (2) the mass ejection is negligible compared to the nucleus, and (3) the ejection velocity is relatively small. These facts contradict the catastrophic breakup of the coreless spheroid observed in Figure 8, but are consistent with our model of an internal core. It implies that these active asteroids may have internal cores which prevent them from breaking up catastrophically.

Additionally, simulations allowed us to observe that the failure due to rotational deformation has a local nature. The local deformation of the surface shell is dependent on its structure, especially on the spatial distribution of particles which is translated in the spatial variation of porosity and the possible presence of fractures. Even though the settling procedure, devoided of cohesive and frictional forces, encourages structural homogeneity, this is not perfect and, at high enough spin rates, granular flows always occur locally. We consider this to be the initial deformation mode. Though long-term deformation modes are not the focus of this paper, there are two possible scenarios for them. In the first scenario, material flows occur globally, causing the formation of a uniform equatorial ridge. This scenario corresponds to the result by obtained Walsh et al. (2008, 2012). In the second scenario, a local deformation leads to further structural bias, causing the further emergence of local material flows. We will conduct further investigations for this problem in the future.

The formation process of a strong core is beyond our scope; thus, we do not quantitatively analyze it in this paper. However, there may be two likely scenarios. The first scenario is that at an accretion process after a catastrophic impact, due to the difference of energy dissipation, large boulders and rocks could accumulate first and then small particles might cover these large objects. Such a formation process could produce bodies with layers of different strength similar to the model we have analyzed (D. Durda & K. Walsh 2015, private communication). Note that since mechanical strength depends on many different parameters such as bulk density (Hirabayashi 2014), friction angle and cohesion, we cannot conclude that such a structure gives a strong core. We leave this problem for future works. The second scenario is based on the assumption that the distribution of loose materials is uniform initially. This condition leads to compression of loose particles in the interior. When such particles irreversibly move out of their original locations, the volumetric strain becomes negative as compression increases (Lambe & Whitman 1969; Nedderman 2005). Because of the self-gravity, the interior experiences larger compression than the surface. As a result, particles in the interior are highly packed. Since higher compaction gives higher cohesion (Lambe & Whitman 1969), this scenario could produce a strong core in an asteroid.

Finally, we compare the HSDEM by Walsh et al. (2008, 2012) with our model. Walsh et al. (2008, 2012) observed surface shedding from a rubble pile body with a homogeneous hexagonal closed packing (HCP) structure,⁸ while we concluded that a body should have an internal core to have surface shedding. This different result comes from the fact that the simulations carried out by Walsh et al. (2008, 2012) and our model may have different deformation features. If particles are arranged to form a crystalline lattice, each particle dislocation is highly controlled by its closest neighbors and the deformation is not isotropic (Ashby et al. 2014). However, even for such a configuration, particles sitting on the surface have open spaces where they can move freely; thus, their motion is less restricted than if they were in the interior. Consequently, the initial HCP configuration could cause anisotropic flows and provide their aggregates with a strong core and a weak surface shell implicitly. Note that Walsh et al. (2012) also investigated a near-fluid case that was modeled by the bi-modal size distribution (a randomly packed aggregate) and obtained a pancake shape, which might indicate an isotropic deformation mode.

This paper quantitatively analyzed the effect of the internal core of a rubble pile spheroid on the mechanism of surface shedding. On the assumption that the friction angle and the bulk density are uniform over the whole volume, we introduced a simple surface shedding model. It is composed of two different cohesion layers: the internal core and the surface shell. The analytical method was constructed by the lower and upper bound theorems in limit analysis, while the numerical method was based on the SSDEM by Sánchez and Scheeres. The failure conditions derived by the analytical model were consistent with those by the SSDEM. The results showed that an aggregate with homogeneous cohesion would always fail internally before it fails superficially. To have surface shedding, therefore, the progenitor body needed a strong core so that only the surface region failed structurally. The SSDEM also exposed the dynamics of the initial failure phase. The primary result was that for all the cases, a rubble pile spheroid would fail locally through a severe material flow on its surface. This came from the fact that since the initial body configuration was not perfectly uniform, which we believe is more realistic, some surface regions were susceptible to fail earlier than others. Further research would be able to provide a better understanding of the formation of unique shapes of asteroids and the mechanisms that originate the activity observed in active asteroids.

The authors wish to thank Dr. Kevin Walsh and Dr. Dan Durda at SwRI Boulder, Dr. Dave Jewitt at UCLA, and Dr. Tony Dobrovolskis at SETI Institute for useful comments on the present work.

The averaged stresses over the internal core and the surface shell are given by the integrations of Equations (2) through (7) over these volumes. Since the shear stress components given by Equations (5) through (7) are trigonometric functions, they become zero by the integration. On the other hand, the normal stress components, Equations (2) through (4), have non-zero values. Using the indices notations defined above, we obtain the averaged normal stresses for the surface shell as

$$\begin{eqnarray}{\bar{\sigma }}_{{\rm{s}}1} & = & \displaystyle \frac{3}{4\pi ({R}^{3}-{R}_{b}^{3})}{\displaystyle \int }_{{R}_{b}}^{1}{\displaystyle \int }_{-\displaystyle \frac{\pi }{2}}^{\displaystyle \frac{\pi }{2}}{\displaystyle \int }_{0}^{2\pi }{\sigma }_{1}\\ & & \times (r,\theta ,\psi ){r}^{2}\mathrm{cos}\theta d\psi d\theta {dr},\\ & = & \displaystyle \frac{3{k}_{1}}{1-{R}_{b}^{3}}\left(\displaystyle \frac{2}{15}-\displaystyle \frac{{R}_{b}^{3}}{3}+\displaystyle \frac{{R}_{b}^{5}}{5}\right)\\ & & +\displaystyle \frac{{b}_{x}-2{k}_{1}}{5}\displaystyle \frac{1-{R}_{b}^{5}}{1-{R}_{b}^{3}},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}{\bar{\sigma }}_{{\rm{s}}3} & = & \displaystyle \frac{3}{4\pi ({R}^{3}-{R}_{b}^{3})}{\displaystyle \int }_{{R}_{b}}^{1}{\displaystyle \int }_{-\displaystyle \frac{\pi }{2}}^{\displaystyle \frac{\pi }{2}}{\displaystyle \int }_{0}^{2\pi }{\sigma }_{3}\\ & & \times (r,\theta ,\psi ){r}^{2}\mathrm{cos}\theta d\psi d\theta {dr},\\ & = & \displaystyle \frac{3{k}_{9}}{1-{R}_{b}^{3}}\left(\displaystyle \frac{2}{15}-\displaystyle \frac{{R}_{b}^{3}}{3}+\displaystyle \frac{{R}_{b}^{5}}{5}\right)\\ & & +\displaystyle \frac{{b}_{z}-2{k}_{9}}{5}\displaystyle \frac{1-{R}_{b}^{5}}{1-{R}_{b}^{3}},\end{eqnarray} \tag{ 29 }$$

and those for the internal core as

$$\begin{eqnarray}{\bar{\sigma }}_{{\rm{c}}1} & = & \displaystyle \frac{3}{4\pi {R}_{b}^{3}}{\displaystyle \int }_{0}^{{R}_{b}}{\displaystyle \int }_{-\displaystyle \frac{\pi }{2}}^{\displaystyle \frac{\pi }{2}}{\displaystyle \int }_{0}^{2\pi }{\sigma }_{1}(r,\theta ,\psi ){r}^{2}\mathrm{cos}\theta d\psi d\theta {dr},\\ & = & {k}_{1}\left(1-{R}_{b}^{2}\right)+\displaystyle \frac{{b}_{x}{R}_{b}^{2}}{5},\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}{\bar{\sigma }}_{{\rm{c}}3} & = & \displaystyle \frac{3}{4\pi {R}_{b}^{3}}{\displaystyle \int }_{0}^{{R}_{b}}{\displaystyle \int }_{-\displaystyle \frac{\pi }{2}}^{\displaystyle \frac{\pi }{2}}{\displaystyle \int }_{0}^{2\pi }{\sigma }_{3}(r,\theta ,\psi ){r}^{2}\mathrm{cos}\theta d\psi d\theta {dr},\\ & = & {k}_{9}\left(1-{R}_{b}^{2}\right)+\displaystyle \frac{{b}_{z}{R}_{b}^{2}}{5}.\end{eqnarray} \tag{ 31 }$$

Note that ${\bar{\sigma }}_{{\rm{s}}1}={\bar{\sigma }}_{{\rm{s}}2}$ and ${\bar{\sigma }}_{{\rm{c}}1}={\bar{\sigma }}_{{\rm{c}}2}$. Since the shear stress components are zero, these components are considered to be the principal stress components.