The Dynamical Complexity of Surface Mass Shedding from a Top-shaped Asteroid Near the Critical Spin Limit

This study focuses on top-shaped asteroids that have frequently been observed in recent decades, such as 1999 KW4, 101955 Bennu, 2008 EV5, and 65803 Didymos (Benner et al. 2010; Statler 2015). These asteroids have diameters ranging from hundreds of meters to kilometers, have relatively high spin rates, and similar overall shapes, i.e., moderate oblateness, elevated bulge at the equator, and small-scale topographic asymmetry about the spin axis. We treat the considered body as a top-shaped homogeneous asteroid uniformly rotating about the maximum principal axis of inertia, which is a long-term stable state assuming that there is some internal frictional dissipation. Mass shedding is expected to occur with sufficient rotation speed, but it can be strongly limited when cohesion is present, e.g., asteroids Didymos (primary) and 1950 DA both have very high spin rates (2.26 hr in period for Didymos, and 2.12 hr for 1950 DA), which are supposed to be fast enough to cause overall reshaping if no cohesion is present (Hirabayashi & Scheeres 2014; Zhang et al. 2018).

In this section, we develop a dynamical model of a regolith grain moving on and off the surface of an asteroid, which is considered to be a spinning gravitational body with arbitrary shape. To keep it a general discussion, we adopt dimensionless quantities in the formulation. The rotational period T is taken as the timescale, and the equivalent radius (the radius of a sphere with equal volume) is taken as the length scale. The accelerations of a unit-mass particle moving on (subscript s) and off (subscript o) the surface are represented in Equations (1) and (2), respectively,

$$\begin{eqnarray}&&{{\boldsymbol{a}}}_{s}=-2{\boldsymbol{\omega }}\times {\boldsymbol{v}}-{\rm{\nabla }}V+{\boldsymbol{f}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{{\boldsymbol{a}}}_{o}=-2{\boldsymbol{\omega }}\times {\boldsymbol{v}}-{\rm{\nabla }}V,\end{eqnarray} \tag{ 2 }$$

where ${\boldsymbol{f}}$ in Equation (1) indicates the acceleration of the particle affected by the contact with the surface, which consists of the normal contact force and the tangential stick-slip force (Schwartz et al. 2012). V is the potential of the centrifugal force and the gravity of the body, defined by Equation (3),

$$\begin{eqnarray}&&V=-\displaystyle \frac{1}{2}({\boldsymbol{\omega }}\times {\boldsymbol{r}})\cdot ({\boldsymbol{\omega }}\times {\boldsymbol{r}})-{{GT}}^{2}\sigma \iiint \displaystyle \frac{1}{| {\boldsymbol{r}}-{\boldsymbol{\rho }}| }d{\boldsymbol{\rho }},\end{eqnarray} \tag{ 3 }$$

where ${\boldsymbol{r}}$, ${\boldsymbol{v}}$ indicate the position and velocity of the particle with respect to the body-fixed frame. ${\boldsymbol{\omega }}$ is the vector of angular velocity, which is a constant magnitude of 2π in the scaled representation. $d{\boldsymbol{\rho }}$ indicates an infinitesimal volume with position vector ${\boldsymbol{\rho }}$, which is integrated over the three-dimensional region occupied by the body. G indicates the gravitational constant, and σ is the bulk density of the body. Noting that the coefficient of the second term is dimensionless, we introduce a factor κ below to describe the relative strength of the gravity compared to the centrifugal force:

$$\begin{eqnarray}&&\kappa =\displaystyle \frac{{{GT}}^{2}\sigma }{3\pi }.\end{eqnarray} \tag{ 4 }$$

The factor κ measures the dynamical environment near a homogenous rigid body rotating around a fixed axis at a constant spin rate. For a homogenous rigid sphere, as the body is spun up, the net gravitational and centripetal acceleration first becomes zero at the equator when κ = 1, which is hereafter defined as the reference value of the normalized critical spin limit. A non-cohesive grain will be gravitationally bound to the surface if κ > 1. Yu (2016) calculated the values of κ for 24 known asteroids based on measured or estimated density and period, showing that their κ values are distributed over a wide range from 1.68 to 193.16. We adopt κ as an intrinsic number of the dynamical system. In the following discussion, we use it as a key system parameter throughout all formulas; and in the discussion of some real asteroid systems, we give both the κ values and properties that have direct physical meaning.

The asteroid shape is described as a parametric surface S(θ, ϕ), in which θ, ϕ are the local angular coordinates of the surface. It offers a complete representation for arbitrary "star-shaped" body, which refers to a limited class of objects for which there exists a point in the interior of the object, standing on which the whole surface is visible. The system dynamical model is then completely determined by S(θ, ϕ) and κ. Although the YORP effect imparts a secular change in the rotational state (Rubincam 2000; Lowry et al. 2014), numerical experiments show that non-catastrophic reshaping occurs intermittently and the structural disturbance of each reshaping only lasts for a short period in comparison (Walsh et al. 2008). In reality, the lifetimes of shed-off regolith particles may cover a wide range, depending on their size, initial positions, and velocities, and other physical properties. We concentrate on hours to months after the particles detach from the body, which is considered to be a violent stage of the body's evolution. Within this period, the orbits of the particles will be restricted to the vicinity of the body, and the motion is dominated by the gravity and centrifugal force (we consider the case of a small total shedding mass compared to the mass of the body). The shape and the spin rate of the body are both assumed to be constant within this short period, so the orbital evolution is essentially governed by Equation (2). Given that a YORP spin-up/spin-down timescale for kilometer-sized near-Earth asteroids (NEAs) is 10⁴–10⁶ years (Rubincam 2000), we can see that the factor κ changes quite slowly versus the dynamical timescale of the particle in the vicinity of the surface. Thus in this study, we adopt a specific value of κ for each short-term simulation (<1 year), but consider different κ values in parallel to reflect representative moments of the whole geologic history of a spun-up asteroid.

We investigate the shape of the primary of the Didymos system as an example. Combined radar and photometric observations show that the primary has typical features of a top-like shape, and is rotating at a critical speed (∼2.26 hr in period), with the spin vector lying near (310°, −84°) in J2000 ecliptic coordinates (Pravec et al. 2006; Scheirich & Pravec 2009; Benner et al. 2010). Figure 1 shows the shape model of the primary in spherical harmonics, truncated to degree 20. The coefficients of the harmonics were derived by applying a fitting algorithm to the Cartesian coordinates of the vertices of the polyhedron model (Benner et al. 2010). The parametric surface of the basic shape can be formulated as

$$\begin{eqnarray}&&S(\theta ,\phi )=\displaystyle \sum _{l=0}^{\infty }\displaystyle \sum _{m=-l}^{l}{c}_{l}^{m}{Y}_{l}^{m}(\theta ,\phi ),\end{eqnarray} \tag{ 5 }$$

where c_l^m are complex coefficients of the harmonics and the bases Y_l^m take the form

$$\begin{eqnarray}&&{Y}_{l}^{m}(\theta ,\phi )=\sqrt{\displaystyle \frac{2l+1}{4\pi }\displaystyle \frac{(l-m)!}{(l+m)!}}{P}_{l}^{m}(\cos \theta ){e}^{{\rm{i}}m\phi },\end{eqnarray} \tag{ 6 }$$

in which P_l^m is the associated Legendre function. Equation (5) presents a general formulation for the star-shaped surfaces (Floater & Hormann 2004), i.e., it establishes an injective map from the spherical angles (θ, ϕ) to the radial distance of the projected point on the surface. So when we introduce a unit orientation vector $\hat{{\boldsymbol{r}}}(\theta ,\phi )$, the position vector of the projected point is ${\boldsymbol{S}}=S\cdot \hat{{\boldsymbol{r}}}$. Figure 1 also defines the geographic coordinate system: the x-axis points to longitude 0°, the y-axis points to the eastern hemisphere, and the z-axis points to the North Pole. We apply this general model to explore the role of asymmetric shape S(θ, ϕ), and in particular, that of the induced irregular gravitational field near the equator. Aiming for a generalized understanding of the dynamics near an asymmetric body close to its spin limit, we discuss the orbital behaviors of lofted particles as a function of the dimensionless factor κ. The reshaping effects of the body owing to a changing κ are ignored, because 1. the κ range explored in this study is narrow, i.e., "near the spin limit"; and 2. we seek to examine the influence of a single variable in this paper. An investigation of the influence of the shape modification is left to future work.

Zoom In Zoom Out Reset image size

We consider loose regolith particles resting on the surface of the considered body with no cohesive forces. As the spin rate is near the critical limit, instantaneous separation from the body can be triggered easily by various low-energy events, e.g., landslides that accelerate the particles to escape speed, or propagation of seismic quakes that give initial lifts. Otherwise, particles that gain no launching energy are trapped on the surface until the spin rate is high enough for the gradient direction of the geopotential to become outward (making the local surface unstable so that no loose particles can attach to it). This section takes into account the detailed shape model of the primary of Didymos and discusses the mass shedding conditions of regolith material in terms of local topography. The spherical harmonic shape S(θ, ϕ) and a polyhedral gravitational field are applied to analysis of the mass shedding scenario, which occurs in two cases: particles being levitated when the local surface becomes unstable, as stated above, and particles being involved in a slope failure and boosted to a speed that exceeds the local escape speed.

2.2.1. Shedding-off Due to the Unstable Surface

The normalized potential V is called geopotential in the vicinity of the considered body. The surface geopotential is defined as

$$\begin{eqnarray}&&{V}_{s}=V\,\circ \,{\boldsymbol{S}}(\theta ,\phi ),\end{eqnarray} \tag{ 7 }$$

which is a fundamental function that describes the tendency of migration of the dense flow moving on the asteroid surface, i.e., regolith particles tend to migrate from high geopotential to low geopotential on the surface. The operator ◦ indicates the composition of two functions. Note that in our model the distribution of V_s is determined by both S(θ, ϕ) and κ. The minimum V_s on an oblate shape is not always located near the equator, which occurs only if κ reduces below a limit. For very high κ values, as the latitude increases, the decrement of the gravitational potential is larger than the increment of centripetal potential, in which case the minimum V_s appears near the polar region. Scheeres (2015) explained the role of geopotential in trapping the shedding regolith in the vicinity of the body. It was shown that as the spin rate increases, the Roche lobe of a spherical body first touches the high-latitude region, and particles in these regions will have sufficient energy to start moving, and probably escape from the asteroid.

We first use a direct method to derive the levitating condition of a cohesionless particle on a surface with the shape of Didymos. When we assume that no sliding occurs before shedding, the levitation occurs when the local geopotential force −∇V switches from inward to outward, i.e., the normal components of the gravitational and centrifugal forces balance at the moment of the switching:

$$\begin{eqnarray}&&{\rm{\nabla }}{V}_{s}\cdot \hat{{\boldsymbol{n}}}=0,\end{eqnarray} \tag{ 8 }$$

where $\hat{{\boldsymbol{n}}}$ indicates the unit normal vector of the local surface (pointing outward). ∇V_s indicates the gradient of the geopotential at the surface. Equation (8) defines the levitating condition. Particles in the region subject to ${\rm{\nabla }}{V}_{s}\cdot \hat{{\boldsymbol{n}}}\leqslant 0$ will be lofted if there are no cohesive interactions with the surface, and in contrast, particles in the region ${\rm{\nabla }}{V}_{s}\cdot \hat{{\boldsymbol{n}}}\gt 0$ are still allowed to be bound to the surface without tensile strength. For the loose particles that were sitting on the surface with weak or no cohesion, the levitation can only lead to small hops of low relative speeds, which are termed "saltation" by Harris et al. (2009). The subsequent evolution relies on the coupling effects with the complex morphology and mechanics of the surface. Taking the primary of Didymos, for example, according to the bulk density 2.146 g/cc and rotation period 2.26 hr derived from observations (Michel et al. 2018), the dimensionless factor κ ≈ 1.0048, which we denote κ*, and which is slightly above the reference value 1. We confirm that for this fast rotator the surface geopotential V_s shows a largely monotonous decrease as the latitude decreases, and the minimum value is around the equator, which is consistent with the oblate nature of the shape.

To explore the levitation behavior, we change κ within a range close to the critical spin limit, i.e., from κ = 0.21 to κ = 1.06. The boundaries of stable regions are determined by Equation (8). Figure 2 illustrates the solutions for different κ values and gives the results as functions of the geographic coordinates. Areas enclosed by the solution curves correspond to ${\rm{\nabla }}{V}_{s}\cdot \hat{{\boldsymbol{n}}}\lt 0$. Note that for a value as low as κ = 0.21, the rotation period is 1.04 hr, which in observations is only found for very small asteroids (<200 m in size) that must have non-zero tensile strength. For comparison, the color map in Figure 2 shows the surface distribution of the scaled elevation, which measures the small-scale topographic features of Didymos shape model.

Zoom In Zoom Out Reset image size

We find that the value κ* derived from observations leads to the consequence that the equatorial region is mostly unstable. The narrow band around the equator, as outlined by the dashed curves in Figure 2, has larger centrifugal force than gravity along the local normal directions, except for a small area in the neighborhood of (135°W, 0°). This suggests that any cohesionless material should be rapidly cleared in this region, exposing unweathered subsurface material; and a loose regolith layer is only possible at the equator near longitude 135°W, which, as illustrated in Figure 1(b), is an obvious low-elevation area. To make this analysis more general, we intentionally excluded the effect of the secondary. Even so, we cannot rule out the possibility that our finding of an unstable equator could be caused by modeling errors, because we used the polyhedral methodology to estimate the gravity near the surface (Werner & Scheeres 1997). Moreover, the shape model based on remote observation does not perfectly reconstruct the small-scale topology, and the bulk density has estimated uncertainties as large as 20% (Zhang et al. 2017; Michel et al. 2018).

Figure 2 reveals some interesting dependencies. First, the unstable region expands from low latitude to high latitude as κ decreases, i.e., while the rotation is accelerating, loose regolith particles originally resting in low-latitude regions tend to be displaced or shed first. Since the equatorial ridge is interpreted as a structure formed in a YORP spin-up, it may also disrupt as the spin is further accelerated (in the absence of tensile strength). At κ = 0.85 (2.08 hr period), the unstable region occupies the region between 20°S and 20°N, which is essentially the entirety of the equatorial ridge. Second, the bulges of the shape correspond to favored regions for mass loss, compared with the flat and depressed regions. It is evident that on the equator the bulged terrains reach the levitation condition prior to other terrains; e.g., when κ ≥ 1.06 (period ≥2.32 hr), the unstable regions split into a chain of isolated islands, and each of these islands corresponds to a bulged place at the equator. It implies a basic mechanism in the formation of the top shape: regolith particles are easier shed from the bulged terrains, and are easier deposited on the depressed terrains, which helps the equator to evolve into a basically round shape.

2.2.2. Shedding Due to the Slope Failure

The analysis of Section 2.2.1 is based on the hypothesis that no tangential slips occur before the surface geopotential reverses, which is not always true because landslides may occur before the spin rate satisfies Equation (8). Nevertheless, we realize that landslides do not have to take place even when the local slope exceeds the angle of repose β_r, not only because of a plausible existence of cohesive forces, but also because of the diversity of local small-scale topographies, e.g., the regolith particles can be geometrically stuck in tangential directions.

On the other hand, considering that there is a tangential speed in regolith flows (landslides), the shedding-off condition will become more variable. The local slope angle is defined by

$$\begin{eqnarray}&&\beta ={\rm{acos}}\left(\displaystyle \frac{-{\rm{\nabla }}{V}_{s}\cdot \hat{{\boldsymbol{n}}}}{| -{\rm{\nabla }}{V}_{s}| }\right),\end{eqnarray} \tag{ 9 }$$

in which 0 ≤ β ≤ π, and if β ≥ π/2, the levitating condition is met. The tangential component of the geopotential force at the surface is

$$\begin{eqnarray}&&{\boldsymbol{l}}=-{\rm{\nabla }}{V}_{s}+({\rm{\nabla }}{V}_{s}\cdot \hat{{\boldsymbol{n}}})\hat{{\boldsymbol{n}}}.\end{eqnarray} \tag{ 10 }$$

The vector $\hat{{\boldsymbol{l}}}$ is the unit vector that determines the steepest descent direction of the geopotential on the surface (also known as the "slope direction"). We assume that the regolith flow is down the slope direction when a landslide happens. It is a simplified assumption that applies to a short-lived landslide, but for grains in a long flow, the direction could be altered strongly by Coriolis forces (Statler et al. 2014). The velocity of a particle moving on the surface is given in the surface coordinates (θ, ϕ) by

$$\begin{eqnarray}&&v\hat{{\boldsymbol{l}}}={{\boldsymbol{S}}}_{\theta }\dot{\theta }+{{\boldsymbol{S}}}_{\phi }\dot{\phi },\end{eqnarray} \tag{ 11 }$$

where the subscripts denote the partial derivative with respect to θ and ϕ. The operator $\dot{\ }$ indicates derivatives with respect to time. v is the magnitude of the velocity. A surface particle undergoing a landslide is governed by Equation (1), and the condition for a lifting off from the surface is

$$\begin{eqnarray}&&{\boldsymbol{f}}\cdot \hat{{\boldsymbol{n}}}=0,\end{eqnarray} \tag{ 12 }$$

which is, at the moment the normal component of the coupling force reaches zero, when the particle is viewed as being separated from the surface. In order to solve Equation (12), we have to introduce the local geometry of the shape model. Combining Equations (1) and (12), and using the parametric expression ${\boldsymbol{S}}(\theta ,\phi )$, we can derive the particle-lifting condition in terms of (θ, ϕ). The first and second fundamental forms of the parametric surface ${\boldsymbol{S}}(\theta ,\phi )$ are defined by

$$\begin{eqnarray}\{{h}_{{ij}}\}=\left[\begin{array}{cc}{{\boldsymbol{S}}}_{\theta }\cdot {{\boldsymbol{S}}}_{\theta } & {{\boldsymbol{S}}}_{\theta }\cdot {{\boldsymbol{S}}}_{\phi }\\ {{\boldsymbol{S}}}_{\theta }\cdot {{\boldsymbol{S}}}_{\phi } & {{\boldsymbol{S}}}_{\phi }\cdot {{\boldsymbol{S}}}_{\phi }\end{array}\right],\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}\{{g}_{{ij}}\}=\left[\begin{array}{cc}{{\boldsymbol{S}}}_{\theta \theta }\cdot \hat{{\boldsymbol{n}}} & {{\boldsymbol{S}}}_{\theta \phi }\cdot \hat{{\boldsymbol{n}}}\\ {{\boldsymbol{S}}}_{\theta \phi }\cdot \hat{{\boldsymbol{n}}} & {{\boldsymbol{S}}}_{\phi \phi }\cdot \hat{{\boldsymbol{n}}}\end{array}\right].\end{eqnarray} \tag{ 14 }$$

$\{{h}_{{ij}}\}$ and $\{{g}_{{ij}}\}$ are metric tensors that together determine the extrinsic invariants of ${\boldsymbol{S}}(\theta ,\phi )$, i.e., the curvature at arbitrary position on the surface. The detailed expansions of the first and second derivatives of the vector-valued function ${\boldsymbol{S}}(\theta ,\phi )$ are given in the Appendix. We can then define a 2 × 1 array $\{{b}_{i}\}$ consisting of the projections of the two tangential vectors onto the slope direction $\hat{{\boldsymbol{l}}}$:

$$\begin{eqnarray}&&\{{b}_{i}\}=\left[\begin{array}{c}{{\boldsymbol{S}}}_{\theta }\cdot \hat{{\boldsymbol{l}}}\\ {{\boldsymbol{S}}}_{\phi }\cdot \hat{{\boldsymbol{l}}}\end{array}\right].\end{eqnarray} \tag{ 15 }$$

Finally, a 3 × 3 orientation matrix {c_ij} is defined to describe the relation between the angular momentum of the asteroid, the local normal vector at the surface, and the local slope direction:

$$\begin{eqnarray}\{{c}_{{ij}}\}=\left[\begin{array}{ccc}{\boldsymbol{\omega }} & \hat{{\boldsymbol{l}}} & \hat{{\boldsymbol{n}}}\end{array}\right].\end{eqnarray} \tag{ 16 }$$

Substituting the expression of ${\boldsymbol{f}}$ (derived from Equations (1) to (12)), and combining Equations (13)–(16), the particle-lifting condition can be written as

$$\begin{eqnarray}&&\,{\rm{\nabla }}{V}_{s}\cdot \hat{{\boldsymbol{n}}}+2{\rm{\det }}(\{{c}_{{ij}}\})v\\ &&\,+\,{\{{b}_{i}\}}^{T}\,{\rm{inv}}{(\{{h}_{{ij}}\})}^{T}\{{g}_{{ij}}\}{\rm{inv}}(\{{h}_{{ij}}\})\{{b}_{i}\}{v}^{2}=0,\end{eqnarray} \tag{ 17 }$$

which is a quadratic polynomial of v. Equation (17) shows that the distribution of solutions depends on the local curvatures, the local slope, and the local orientation with respect to the rotating direction. Physically permitted landslide speeds that lead to a separation from the surface correspond to the positive real roots of Equation (17), which are defined hereafter as the separation speeds for regolith material sliding down in the slope direction.

Figure 3 illustrates the distribution of the separation speed for Didymos, with a contour plot of the slope angle, for four cases of different spin rates (κ values). The dotted, dashed, and solid curves indicate the boundaries of regions where the slope angle is below the specified value β_r (the angle of repose)—i.e., in the regions satisfying β < β_r loose regolith is allowed to be deposited—while beyond these regions (β > β_r), the regolith material tends to migrate because of local landslides.

Zoom In Zoom Out Reset image size

A nonlinear scale is adopted for the color map in Figure 3, in which separation speeds above 30 cm s⁻¹ are all marked in monocolor. The threshold value 30 cm s⁻¹ is sufficiently high to exceed, according to our calculation, the maximum expected speed of a true landslide. "Dead zones" are found around the polar sites in all four cases of Figure 3 (marked in black; the north polar zone is larger than the south polar zone); particles in these regions are required to have infinitely high speeds to separate from the surface ($v\to +\infty$), i.e., the local landslides would never lead to mass shedding in these areas (note that we consider material that carries only tangential velocity on the local surface). The "dead zones" have emerged around the polar sites of Didymos (primary), and their areas shrink as the body is spinning up. These depressed areas at the polar sites of Didymos (see Figure 1) show the capability of keeping the sliding material bound to the local surface, which is independent of the sliding speed down the slope.

Figure 3 reveals several explicit dependencies of the separation speed on κ. First, the surface separation speed generally increases with increasing latitude. For κ = 0.85 and κ = κ*, the minimum separation speed is distributed along the edge of the unstable band (see Section 2.2.1); for κ = 1.27 and κ = 1.70 (when the unstable strip disappears), the minimum separation speed is distributed along the equator. This means that during slope failure, the loose regolith at low latitudes requires lower speeds to be shed from the surface than that at high latitudes, which is consistent with the levitating condition. Second, the separation speed decreases rapidly as the body is spun up. According to the parametric model of Didymos, as the rotation period decreases from 2.94 to 2.08 hr, the minimum surface separation speed drops from 1.73 cm s⁻¹ to 0.12 mm s⁻¹, which was calculated from the dimensionless results in Figure 3. This suggests that the surface becomes increasingly unstable as the body approaches its critical spin limit. Third, a dependence on the local topography is noted, especially at the equatorial ridges. The bulges correspond to local minimum values of the separation speed, and the depressed terrains correspond to local maximum values. This result reflects the correlation of separating condition with local curvature, according to Equation (17), and it is also consistent with the levitating condition. In effect, material is likely to escape from the bulges and be deposited at the depressed terrains, which provides a secular mechanism to force the ridge into good average roundness.

In addition to the separating condition, Figure 3 also shows the contour lines of the slope angle. If we assume that the angle of repose of cohesionless regolith β_r = 45°, which is a rather conservative estimate of the angle of repose based on the values of terrestrial dry sands (a similar analysis applies to β_r = 15° and β_r = 30°), at κ = 1.70, the surface of the primary is globally below β_r, and at κ = 1.27, part of the mid-latitude region exceeds the angle of repose β_r where regional failures seem to occur first. For the rapid rotation cases κ = 8 and κ = κ*, the wide region between 45°S and 45°N has exceeded the angle of repose β_r. This means that a vast region of the surface beyond the polar areas is unstable according to the parametric model, so the loose regolith over this region should have been rearranged and unweathered subsurface material is probably exposed. Combining the steep slopes and low levitating/separating condition, we speculate that cohesionless material has been removed from the mid-low latitudes on the surface of Didymos (primary).

Once the surface shedding conditions are met, particles start to enter a complicated phase of evolution. Based on our model, the lofted particles seem very unlikely to enter a high-eccentric trajectory after the first detachment from the surface. Without cohesive forces, the particles can be levitated or slide down the local slope, but neither is capable of causing a sharp increment in kinetic energy. The ballistic motion of a particle departing from the surface is governed by the local gravity, and subsequent collisions with the surface might occur frequently. After this stage, which is highly coupled with the surface topography, it is possible that the particle becomes trapped by and settles down on the surface again (reaccretion), or gains enough energy that it maintains a temporary low-altitude cycling orbit, i.e., with the periapsis close to the surface of the asteroid. For our analysis, we define a particle to be in cycling motion if no collisions with the surface occur between two consecutive periapsis passages.

Here we consider the energy transported to a unit mass shedding during the scenario described above. Note that the analysis in Sections 2.2.1 and 2.2.2 shows that the shedding process first occurs at the equator, so we make the simple assumption that the lofted particles depart from the equatorial plane. To conduct a simplified analysis, the scaled total energy of the cycling orbit of a particle at the equator of a spinning spherical body is

$$\begin{eqnarray}&&E=2{\pi }^{2}(1-2\kappa ),\end{eqnarray} \tag{ 18 }$$

note E < −2π² if κ > 1. The Keplerian energy of the cycling orbit of a unit mass of the lofted material is

$$\begin{eqnarray}&&E^{\prime} =-\displaystyle \frac{2{\pi }^{2}\kappa }{a},\end{eqnarray} \tag{ 19 }$$

where a is the semimajor axis of the orbit, satisfying $E^{\prime} \geqslant E$. The considered particle does not permanently orbit in the vicinity of the surface, but returns to the surface sometime. We define the energy increment to a single change of the velocity of the particle Δv, and the increment in energy is given as the result of a sudden change in velocity:

$$\begin{eqnarray}&&E^{\prime} -E=\displaystyle \frac{1}{2}{\rm{\Delta }}{v}^{2}+2\pi {\rm{\Delta }}v,\end{eqnarray} \tag{ 20 }$$

and we scale the velocity change via

$$\begin{eqnarray}&&\eta =\displaystyle \frac{{\rm{\Delta }}v}{2\pi \sqrt{\kappa }}.\end{eqnarray} \tag{ 21 }$$

Substituting Equations (18), (19), and (21) into Equation (20), we have

$$\begin{eqnarray}&&a=\displaystyle \frac{\kappa }{2\kappa -{(\eta \sqrt{\kappa }+1)}^{2}}.\end{eqnarray} \tag{ 22 }$$

Because the eccentricity of the orbit is negligible, i.e., e ≈ 0, in this analysis, the dynamical fate of a test particle is only determined by two parameters, a and κ. The Roche lobe defines the spatial region that is energetically accessible to the particle moving around the body (Scheeres 2015). Considering an orbiting particle in the equatorial plane, the radius of the synchronous orbit for the spherical model is given by

$$\begin{eqnarray}&&q=\sqrt[3]{\kappa }.\end{eqnarray} \tag{ 23 }$$

The synchronous orbit is defined as a natural circular orbit about the target body with its period equal to the spinning period. The Roche lobe is represented in the rotating frame. The orbital motion bounded in the vicinity of the asteroid should satisfy the condition

$$\begin{eqnarray}&&\left\{\begin{array}{l}1\lt a\lt q\\ {J}_{a}\lt {J}_{q},\end{array}\right.\end{eqnarray} \tag{ 24 }$$

where J_a and J_q indicate the Jacobi integrals of the orbits of semimajor axis a and the synchronous orbit, respectively. The dimensionless form of the Jacobi integral is given as

$$\begin{eqnarray}&&J=-\displaystyle \frac{4{\pi }^{2}\kappa }{r}-2{\pi }^{2}{r}^{2}+\displaystyle \frac{1}{2}{v}^{2},\end{eqnarray} \tag{ 25 }$$

then substituting the orbital radii and velocities (in the rotating frame) into the second inequality of Equation (24) yields

$$\begin{eqnarray}&&{\left(\sqrt{\displaystyle \frac{\kappa }{a}}-a\right)}^{2}-2\kappa \left(\displaystyle \frac{1}{a}-\sqrt[3]{\displaystyle \frac{1}{\kappa }}\right)-({a}^{2}-\sqrt[3]{{\kappa }^{2}})\lt 0.\end{eqnarray} \tag{ 26 }$$

Figure 4 presents a κ versus a diagram that maps out the equivalent orbital altitude of particles lofted from the equator of the central body versus the rotation period. The semimajor axis a is represented as a function of κ as defined by Equation (22). The parametric function is varied by the energy transported via surface coupling, which is described using η in Equation (22). Figure 4 illustrates the function curves for three η values, 0.0, 0.1, and 0.2. As an estimate of the upper limit, η = 0.2 is capable of lofting particles originating from the equator of the target body with κ < 1.56 (T < 2.82 hr for Didymos), and η = 0.1 is the same for κ < 1.24 (T < 2.50 hr). η = 0.0 is the case of no extra energy transported so that only particles from the body with κ < 1 will be shed. Within the considered parameter range, we find that transported orbits are limited to a region closely above the surface (a < 3.0), and the modified energy is quite insufficient to allow escape from the vicinity of the asteroid.

Zoom In Zoom Out Reset image size

Figure 4 also illustrates the region confined by the Roche lobe in the equatorial plane (shadowed area), which is defined by Equations (24) and (26). Scheeres (2015) stated that the particles that cycle around the body in a region below the synchronous orbit are permanently trapped in the vicinity of the asteroid, and conversely, particles moving in regions beyond are allowed to escape to arbitrary distances away from the body. The shadowed area in Figure 4 shows that as an asteroid is spinning up (or κ is decreasing), the region bounded by the Roche lobe shrinks rapidly, and vanishes at the reference value κ = 1. Thus given the same ratio of energy transport η, a larger κ corresponds to a wider Roche lobe radius than a smaller κ. Statistically, this means that an asteroid in slow rotation has the capability of trapping the lofted particles in its vicinity. Note that this analysis is based on the assumption of equal energy transports via surface coupling. Based on Figure 4, we can infer that the mass leaking will be reinforced as the rotational speed approaches the critical limit, i.e., the dynamical fates of the regolith material shed from the surface will be highly weighted to being transported farther away, which is probably occurring with the Didymos primary considering that κ* is close to 1.

With the dimensionless model developed above, the dynamical flow around the target body is completely determined by its shape and the parameter κ. Considering the asymmetric shape, a great diversity is found in the dynamics of lofted particles, contrasting with the simple orbital motion around an ideal spherical body. This section explores the orbital dynamics of lofted particles near the surface of the progenitor asteroid, and we start with a theoretical survey of the synchronous orbits, changing as functions of κ. The synchronous orbits are isolated in the asymmetric case, rather than continuously distributed over a certain radius. These synchronous orbits are referred to as equilibrium points of Equation (2) in the body-fixed frame, which satisfy

$$\begin{eqnarray}&&{\rm{\nabla }}V=0.\end{eqnarray} \tag{ 27 }$$

Assuming that the shape of the body is constant, we can independently observe the qualitative changes of the equilibrium points as a function of κ. In particular, when κ is approaching the critical spin limit, the equilibrium points are close to the surface of the target body and it will be relatively easy for the lofted particles to reach up to the altitudes of these points. As was shown in Section 2.3, this change in mechanical environment exerts a significant effect on the secular evolution of the lofted particles. Scheeres et al. (2006) identified that the primary of 1999 KW4 has four equilibrium points outside the body, one of which lies just at the surface. In comparison, the primary of Didymos has a poor roundness at the equator, and hence exhibits a more complex geopotential structure around the surface. Figure 5 shows the counter lines of the geopotential on the equatorial plane (gray solid lines), which are plotted for κ = 11 (2.43 hr in period). Eight equilibrium points are found in total in this case, and according to Equation (27), their positions can be determined by searching for the self-crossing points of the countour lines. As shown in Figure 5, all the eight equilibrium points (marked with solid dots) are close to the equator of the asteroid (shadowed area).

Zoom In Zoom Out Reset image size

It has been recognized that equilibrium points always exist outside a uniformly rotating body of arbitrary mass distribution, as long as the rotation is substantially slow. In theory, the position of a given equilibrium point depends on the rotation period, or in the normalized form, on κ. Hence as κ continuously changes, the location of the equilibrium point varies and its locus extends to an arc, which is hereafter named the equilibrium branch. As illustrated in Figure 5, all the eight exterior equilibrium points extend to corresponding branches (dashed lines). These branches extend outward from the surface as κ increases, except when two of them collide and vanish. The positions where such collisions occur are called bifurcation points, marked by circles in Figure 5. In our system, the bifurcations of equilibrium branches essentially define the κ values at which the topology of geopotential structure changes, which implies a direct impact on the motion of the lofted particles (see Section 3.2 for more details).

Jiang et al. (2015) studied both the exterior and interior equilibria about the asteroid 216 Kleopatra, and found that their total number is always odd and changes in pairs as the spin rate varies. We focus on the exterior equilibrium points for practical reasons, and find that the number of such points presents a $4\to 6\to 8$ change for Didymos as κ decreases, which still follows the rule "appear and vanish in pairs." Figure 6 shows the profiles of the eight equilibrium branches on the κ versus $J/{J}_{q}$ diagram, in which J indicates the Jacobi integral of the equilibrium point. The axis $J/{J}_{q}$ defines the relative bias in Jacobi integral with respect to the result for an ideal spherical model. The reason we adopt the scaled index $J/{J}_{q}$ instead of J is that it enables a clear definition of all eight branches in one single diagram, which gives a broad view of how the equilibrium points shift as κ changes. The diagram reveals that A–A', B–B', C–C', and D–D' are the trivial branches that remain as κ increases (the prime notation is used to distinguish between the starting and ending points). Furthermore, these trivial branches converge asymptotically on the benchmark value $J/{J}_{q}=1$ (no bias) when $\kappa \to \infty$. This is because when the equilibrium point moves away from the body, the configuration of the nearby gravitational field can be approximated to that induced by a point of the same mass, or an ideal sphere, which means that the Jacobi integral J is approximately J_q. While the target body is spinning up, equilibrium points H' (or I') first emerge from a bifurcation at κ = 2.44 (3.52 hr period), and split into branches H–H' and I–I' as κ decreases further, bringing the number of equilibria to six. Likewise, the branches F–F' and G–G' appear at κ = 1.24 (2.53 hr period), and bring the number of equilibria to eight as κ further approaches the critical spin limit. Note that not all these branches remain until the critical limit, i.e., some will end up at the asteroid surface before κ reaches 1. In fact, only one exterior equilibrium point remains on branch D–D' when κ = κ*, and its position is almost stuck to the surface of the asteroid shape.

Zoom In Zoom Out Reset image size

Figure 6 also shows the stabilities of the equilibrium branches. Four subbranches are identified to be stable, including C–C' (1.48 < κ < 4.24), D–D' (1.80 < κ < 4.24), G–G' (1.12 < κ < 1.24), and I–I' (1.38 < κ < 2.44). The definition of linearized stability is adopted here, i.e., the stability of an equilibrium point is determined by the local linearized system, which is guaranteed by the center manifold theorem (Shilnikov et al. 1998). In the system described by Equation (2), there are six eigenvalues of the linearized system λ_i (i = 1, ..., 6), which are included in one of the following cases (Yu 2016): I. Two opposite real numbers; II. four opposite conjugate complex numbers; and III. two opposite purely imaginary numbers. Defining a topology index

$$\begin{eqnarray}&&({N}_{r},{N}_{c},{N}_{i}),\end{eqnarray} \tag{ 28 }$$

where N_r, N_c, and N_i are the number of real, complex, and purely imaginary eigenvalues, respectively, there are only six combinations of the three types, as listed in Table 1.

The topological types described in Table 1 represent different manifold structures at the equilibrium point that determine the general behaviors of trajectories in its neighborhood. Shilnikov et al. (1998) discussed in detail the orbital schemes around equilibria of the above six types. We apply the methods of qualitative theory to asteroid Didymos, and found that only three types appear in the eight branches: the (2, 0, 4) saddle-focus, the (0, 4, 2) center-focus, and the (0, 0, 6) center point. Of these three types, only the center point is stable in the sense that particles nearby tend to cycle around it following the local central manifolds. The saddle-focus and center-focus topologies are unstable, and the motion of particles around them is dominated by the unstable local manifolds, i.e., particles depart from them in exponential/spiral ways. Table 2 presents the topological changes of the eight branches around Didymos (primary) as a function of κ. The results essentially reflect in geometry the evolution of the dynamical flow near the equilibria. Within the κ range examined, the topology of branches A–A', B–B', F–F', and H–H' sticks to the saddle-focus type and remains unstable, while the topology of branches C–C', D–D', G–G', and I–I' transfers from the center-focus type to center point type; i.e., the stability accordingly switches from "unstable" to "stable" in these branches. The latter is a rare case in the equilibria bifurcations of Hamiltonian systems because it requires the four imaginary eigenvalues at the bifurcation point to be degenerate, i.e., the two pairs of imaginary eigenvalues coincide on the imaginary axis. Before the stability transition occurs, the equilibrium branch is in complex instability, which exhibits low-speed divergence; after that, the branch is in center-type stability, which is critical and non-asymptotical. Table 2 employs a stability indicator (Shilnikov et al. 1998):

$$\begin{eqnarray}&&{\rm{\Xi }}=\mathop{\max }\limits_{1\leqslant i\leqslant 6}\{{\rm{Re}}({\lambda }_{i})\},\end{eqnarray} \tag{ 29 }$$

which is the maximum of the real parts of the six eigenvalues. For a Hamiltonian system defined by Equation (2), Liouville's theorem ensures that Ξ is always non-negative, therefore an equilibrium point is stable only if Ξ = 0, and is unstable if Ξ > 0. The magnitude of Ξ (>0) reflects the degree of instability of the equilibrium point, i.e., a larger Ξ corresponds to a higher degree of instability. As illustrated in the table, Ξ shows a descending trend as κ increases, except for a few local anomalies. In effect, this suggests that the equilibria on all the eight branches become increasingly unstable as the asteroid is spinning up.

Table 2.  Some Metrics on the Topology and Stability of the Eight Equilibrium Branches about Didymos

Branch	κ	Topology	Ξ	Stability
A–A'	$1.14\to 4.24$	(2, 0, 4)	$7.93\to 1.28$	U
B–B'	$1.11\to 4.24$	(2, 0, 4)	$3.76\to 2.52\to 2.58\to 0.97$	U
C–C'	$1.06\to 1.48$	(0, 4, 2)	$2.70\to 0.00$	U
	$1.48\to 4.24$	(0,0,6)	0	S
D–D'	$0.98\to 1.80$	(0, 4, 2)	$2.68\to 0.00$	U
	$1.80\to 4.24$	(0,0,6)	0	S
F–F'	$1.05\to 1.24$	(2, 0, 4)	$1.40\to 1.80\to 0.44$	U
G–G'	$1.00\to 1.12$	(0, 4, 2)	$1.39\to 0.00$	U
	$1.12\to 1.24$	(0,0,6)	0	S
H–H'	$1.09\to 2.44$	(2, 0, 4)	$4.55\to 0.46$	U
I–I'	$1.03\to 1.38$	(0, 4, 2)	$2.09\to 0.00$	U
	$1.38\to 2.44$	(0, 0, 6)	0	S

Note. The topological type, stability indicator Ξ, and stability type are listed as a function of κ. An arrow is used to indicate a monotonous variation from tail to head. "U" is short for "unstable," and "S" for "stable."

Download table as:  ASCIITypeset image

We considered the implications of these topological transitions in the realistic evolution of asteroids. The lofted particles are allowed to gather and wander about a stable equilibrium point, i.e., to move in groups along the synchronous orbit with respect to the inertia frame. For an equilibrium point with sufficient stable margin, such orbits of particles clustered in its neighborhood may last for a relatively long time (we confirm that the interval can last up to at least several months in the case of Didymos; see Section 3.2). Growth of solid debris probably proceeds via mutual collisions across the region, depending on both the particle sizes and their orbital concentration (particles in this case have low relative speeds). As the debris growth essentially relies on the particle size, these processes (if occurring) can be relatively stochastic (Pater & de Lissauer 2015). Specific to the collisional growth of big particles, e.g., particles >1 cm in diameter, the stable equilibria might provide a mechanism that enables the dispersed particles to gather and aggregate by increasing their collisional probability. Another significant effect is caused by the topological transition, when an equilibrium branch switches from "stable" to "unstable." Particles originally trapped around the equilibrium point will experience a dramatic change of dynamical environment, and would be removed in a short interval (see Section 3.2 for details). This orbital transfer follows the unstable manifolds at the equilibrium point, and may eventually lead to an impact with the surface of the asteroid, to an escape from the system, or to injection into a new temporary cycling orbit. If the amount of transferred material is sufficiently large, each of the three outcomes will be capable of causing a sudden change in the angular momentum of the target body.

We consider a system including only the asteroid and the lofted particles, which enables us to examine the influence of the non-spherical gravity separately from other possible external perturbations. The asteroid is assumed to be rotating about the maximum principal axis of inertia, and its vicinal gravity field is calculated using the polyhedral method. The mass of lofted particles is assumed to be negligible compared to that of the asteroid. For each simulation, tracer particles are first sampled around the neighborhood of the asteroid and assigned initial positions and velocities. Then individual simulations are performed with this setup to track the month-long evolution of these particles. Note that this model applies primarily to centimeter-sized or larger particles because the motion of such particles is dominated by gravity.

For subcentimeter-sized particles from asteroids, the solar radiation pressure can play a significant role in the orbital motion, while for micron- and submicron-sized asteroidal particles, the orbital motion would also be strongly influenced by the Lorentz force if particles become charged (Liu & Schmidt 2018). Poynting–Robertson drag and solar wind drag are important only for the long-term dynamics of asteroidal particles (Yu et al. 2017), so these drag forces are also neglected over the timescales of our simulations (see below). The orbits of tracer particles are solved by integrating Equation (2), which defines Hamiltonian flows exterior to the central body (Yu 2016). Collisions with the surface have to be treated while considering the motion close to the asteroid. Since we focus on the on-orbit behaviors of shed regolith material, a particle is assumed to be stuck on the surface at the impact site when a collision occurs. The Laplacian of the polyhedral potential provides a simple criterion for fast detection of such collisions (Werner & Scheeres 1997), which was adopted in our numerical model.

Our first group of simulations aimed at a broad understanding of the dependency of the orbital behaviors on the system parameter κ, without considering the mutual collisions between the particles. The dynamical fates of lofted particles essentially depend on their initial conditions. We intentionally assigned an initial set of particles that are distributed uniformly over an annular disk lying on the equatorial plane, and the annular disk is originally centered on the asteroid mass center. After removing the particles lying inside the asteroid body, the remaining particles constitute a dense ring around the target body, with the inner edge close to its surface. After the initial positions of the particles are chosen, the corresponding velocities are calculated by assuming that they are on equivalent Keplerian circular orbits, as described in Section 2.3. This choice is reasonable, because we focus on the fates of the temporarily cycling orbits, rather than on those that are briefly lofted and fall back within a short interval. Equation (22) provides a rule to constrain the size of the ring composed of cycling particles. We continue to employ the scaled shape model of the primary of Didymos, and the observed parameter value κ*. Then assuming a large upper limit to a (scaled) velocity change of η = 0.1, the ring radius at the exterior rim lies at 1.2566 in scaled length units, which is ∼486.90 m in terms of the real dimensions of Didymos. Following the methodology stated above, we generated a ring system consisting of 46,160 particles. Simulations were performed to sweep a wide range of κ, i.e., κ* ≤ κ ≤ 4.24.

In this study, the choice of simulated time was made through observing the orbital distribution on the parameter plane of semimajor axis versus eccentricity. The orbital characteristics were calculated according to the approximated Keplerian orbits of the particles cycling around the asteroid (the instantaneous Keplerian orbits as ellipse, parabola, and hyperbola characterized by six elements). We note in all three cases that the trajectories of the particles entered a steady state evolution after 10–15 days, i.e., most tracer particles have either impacted on the surface of the asteroid, or have entered a period of relatively stable motion when the shapes of the approximated Keplerian orbits show slow variation and the instantaneous elements reflect the distribution of the orbits of the surviving particles (see the supplementary materials for details). On this basis, the simulated time was taken as 30 days, which we have determined to be sufficiently long and to enable a complete observation of the orbital evolution of the lofted particles.

Figure 7 shows three representative cases of the evolution of the ring system: κ = κ* (a), κ = 1.70 (b), and κ = 3.71 (c). In each subplot of Figure 7, snapshots were created at early, middle, and late time of the simulation, respectively, showing different configurations of the ring in these stages; the histograms at the bottom record the distributions of the orbits of cycling particles corresponding to the sampled frames. Here the orbits of particles were defined to be the instantaneous Keplerian orbits at a specified moment, which gives a measurable approximation to the trajectory of a lofted particle. Two parameters in Figure 7 are used to measure the shape of a particle orbit: the semimajor axis a describes the dimension of the orbit on average, and the periapse distance r_p describes the closest radial distance from the orbit to the asteroid mass center. Unscaled times are shown at the top left of each snapshot in order to facilitate comparison. Fictitious large particles are drawn in these frames for visual enhancement, and the green particles mark those transported back to the surface (stuck to the impact locations as defined in our model).

Figure 7 demonstrates the diverged dynamical behaviors of the ring corresponding to different rotation speeds. The spin–orbit resonance is found to be a key mechanism that plays an important role in all examined cases; i.e., when the spin period of the asteroid and the orbital periods of the lofted particles are commensurate, a periodic gravitational influence acts on the orbiting particles. The result of these perturbations may drastically modify the orbits of the particles, which will either lead to a collision with the asteroid or transfer the particles into more sustainable orbits.

The Roche lobe is beneath the asteroid surface at Figure 7(a) κ = κ*, thus all the particles are initially above it and have the opportunity to be transported to a region far away from the surface. Numerical simulations confirm that a period of one month is a reasonable interval to enable the ring particles to be scattered. The histograms show that the distribution of the semimajor axis spreads rapidly at the early stage: its upper limit grows from 1.26 to 6.37 in 10 days. At the same time, the range of the periapse distance does not change much, remaining in the 1.25–1.57 range. That is, most particles gained extra energy from orbital resonances with the asymmetric gravitational field, and their orbits underwent mainly eccentric drifting that locked the periapsis to the vicinity of the asteroid and transported the apoapsis far away. These particles in scattered orbits carried more total energy and angular momentum than in their original orbits, and in extreme cases, theory predicts that some particles will gain positive Keplerian energy and enter escaping orbits through a "slingshot" effect, which, however, was never seen in our simulations. This is because the largely round shape of the asteroid leads to relatively weak spin–orbit resonances.

When the rotation speed slows down, Section 3.1 shows that a stable equilibrium point first emerges on branch G–G' at κ = 1.12, which is close to the surface of the target body (see Figure 5). At the value examined in Figure 7(b) κ = 1.70, two stable equilibrium points in total are seen, lying on branch C–C' and branch I–I'. The snapshots demonstrate the dynamical influence of the equilibrium points on the orbital motion of shedding particles, i.e., both give rise to clustering of the orbiting particles that prevents particles near the equilibria from scattering. The histograms in Figure 7 show that in our simulation, over 40% of the particles remain in the clustering structures for at least 30 days. Figure 6 suggests that there can be up to three such clusters when 1.80 < κ < 2.44. Particles above the synchronous orbit have the same behavior as those in Figure 7(a), i.e., the orbits experience eccentric drifting so that the distribution of semimajor axis extends outward rapidly, while the periapse distance remains in the 1.25–1.68 range. On the other hand, the particles originally beneath the synchronous orbit are theoretically bound to the surface, and according to our model, a majority of these particles collide with the surface before the end of the simulation.

A qualitative difference occurs when the whole ring is initially inside the Roche lobe, as shown in Figure 7(c) κ = 3.71: except for the particles that collide on the surface, the other particles are permanently constrained within a narrow range defined by the altitude of the synchronous orbit. The histograms show that no orbital diffusion occurs in this case, i.e., both the semimajor axis and periapse distance preserve their distributions over the 1.25–1.37 range.

The second group of simulations focuses on the motion of the clustered particles near the equilibria, as illustrated in the snapshots of Figure 7(b). Theoretically, the clusters should only emerge around stable equilibrium points, as a result of the orbital transport of shed materials. Catastrophic changes may occur to such clustered particles if the system parameter κ shifts and leads to a qualitative modification of the flow around the equilibria. Section 3.1 described the bifurcations and stability changes of all eight branches of equilibrium points. Numerically, we can further explore the behaviors of the clustering particles, and show when and how often the catastrophic changes will occur in comparison with the YORP timescale.

Examples are provided to represent the particles moving around the equilibrium points on branch D–D'. The initial positions and velocities of particles are randomly chosen subject to a given upper limit, ΔJ, which is chosen to limit the maximum deviation of the Jacobi constant of these test particles from that of the corresponding equilibrium point. In our simulations, we examined sample clusters consisting of 1000 particles in two cases with different spin rates: I. κ = 2.55, which is away from the stability transition point (for branch D–D', the transition point is κ = 1.80; see Table 2); and II. κ = 1.80, which is slightly above the transition point and remains on the stable branch. Simulations of the sample clusters are performed under the original spin rates (i.e., κ values defined above) to create benchmarking results for comparisons and to numerically demonstrate the stability of the cluster during the simulated time (30 days, which is confirmed to be sufficiently long to account for the stability change we try to observe here). Assuming that the asteroid experiences a spin-up process due to the YORP effect, we simulate the evolution of the same clusters under the changed spin rates. We calculate the κ change based on the YORP spin-up/spin-down timescale for kilometer-sized asteroids, which is estimated to range from 10⁴ to 10⁶ year (Rubincam 2000; Ćuk 2007). As an example of a strong YORP effect, i.e., spin-up period $| T/({dT}/{dt})| \,\approx {10}^{4}\,\mathrm{years}$, we calculate that the dimensionless factor κ drops from 2.55 to 2.44 within ∼130 years for case I; or within the same time, drops from 1.80 to 1.74 for case II. Then controlled simulations of the sample clusters are performed under the accelerated spin rates for cases I and II. Figure 8 presents a zoomed view of the position distributions of the particles that belong to the cluster at the original and accelerated spin rates for both cases. For each specified κ value (2.55 or 1.80), the initial configurations of the clusters are colored in red, their configurations after a 30-day-long simulation under the original spin rate are colored in blue, and the configurations after a 30-day-long simulation under the accelerated spin rate are colored in black.

Zoom In Zoom Out Reset image size

Figure 8 shows how the fate of a cluster depends on the stability of the affiliated equilibrium point. In case I, the equilibrium point is far from the unstable branch of D–D', and thus has a relatively wide margin in stability, i.e., the system is structurally stable so that after the asteroid experiences a spin-up, the particles of the cluster will stay around the shifted equilibrium point. In case II, however, the cluster before/after the spin-up exhibits distinct outcomes; i.e., the particles under the original spin rate remain clustered about the equilibrium point after one month, while the particles under the shifted spin rate experience catastrophic disruption within the same period, and a majority of these particles finally collide with the surface of the asteroid. Figure 9 illustrates the disruption process of the simulated cluster that corresponds to case II (shifted equilibrium point). Snapshots are shown at times 0 minute, 8 hr 39 minutes, 15 hr 51 minutes, 1 day 0 hr 30 minutes, 1 day 23 hr 34 minutes, 6 days 20 hr 19 minutes, 13 days 0 hr 48 minutes, 13 days 0 hr 48 minutes, and 30 days, showing the representative configurations of the cluster.

Zoom In Zoom Out Reset image size

Figure 9 demonstrates that the cluster starts to spread out of the neighborhood of the equilibrium point within hours after the initial time. The disruption begins just hours after the first expansion, during which the cluster is stretched and distorted rapidly. Most particles of the cluster are dispelled from the equilibrium point in one week, and land on the close side of the surface. This numerically confirms that when the system is evolving into the critical spin limit, the equilibrium point is no longer structurally stable and a slight change in the spin rate may lead to catastrophic disruption of the equilibrium cluster. Such a stability transition of the equilibrium point provides a mechanism to rapidly remove the lofted particles that are previously transported to the near-synchronous orbit. Because of this effect, it is hard to accrete a satellite around a fast-rotating body (close to the critical spin limit) or to capture it into an orbit close to its surface, i.e., the lofted regolith will be transported to a large spatial region and may finally lead to "mass leaking" from the system.

We further analyze the accretion of particles that are clustered around an equilibrium point. It is known that the collisional growth of the grains essentially depends on the particle size, the mean relative speeds, and the volume number density of particles in the cluster. The collision rate increases with increasing volume density of particles in the cluster. The settling speed (the relative speed that enables two particles to agglomerate rather than bounce off after the collision) is correlated with the particle size in a complicated way depending on the scale of particle size, and it may be influenced by cohesion forces. As a general discussion regarding particles clustered around the equilibrium point, the volume density highly depends on the topology of the equilibrium point, because the motion of the particles, except when there are collisions, is governed by Equation (2) (if no extra perturbations are considered). Statistically, this means the spatial region occupied by the cluster evolves under the control of the dynamical flow. The geometry of the dynamical flow near the equilibrium point was discussed in Section 3.1. Here we consider a small volume originating near the equilibrium point, w; the evolution of w is approximated using the linearized system of Equation (2) (see a detailed discussion in Yu & Baoyin 2012). The linearized neighborhood is transported via the Jacobian matrix, which is a 6 × 6 matrix defined as

$$\begin{eqnarray}{\{{J}_{{ij}}\}}^{t}=\left[\begin{array}{cc}\displaystyle \frac{\partial {{\boldsymbol{r}}}^{t}}{\partial {{\boldsymbol{r}}}^{0}} & \displaystyle \frac{\partial {{\boldsymbol{r}}}^{t}}{\partial {{\boldsymbol{v}}}^{0}}\\ \displaystyle \frac{\partial {{\boldsymbol{v}}}^{t}}{\partial {{\boldsymbol{r}}}^{0}} & \displaystyle \frac{\partial {{\boldsymbol{v}}}^{t}}{\partial {{\boldsymbol{v}}}^{0}}\end{array}\right].\end{eqnarray} \tag{ 30 }$$

The superscripts 0, t indicate the time when a variable is accessed. ${\{{J}_{{ij}}\}}^{t}$ is evaluated on a trajectory segment that starts at state $({{\boldsymbol{r}}}^{0},{{\boldsymbol{v}}}^{0})$ and ends at state $({{\boldsymbol{r}}}^{t},{{\boldsymbol{v}}}^{t})$. If the trajectory segment is specific to the equilibrium point, the Jacobian matrix takes the explicit form

$$\begin{eqnarray}&&{\{{J}_{{ij}}\}}^{t}={\rm{\exp }}(\{{a}_{{ij}}\}t),\end{eqnarray} \tag{ 31 }$$

in which $\{{a}_{{ij}}\}$ is the linearized matrix of Equation (2) at the equilibrium point, i.e.,

$$\begin{eqnarray}\{{a}_{{ij}}\}=\left[\begin{array}{cc}{\rm{O}} & {\rm{I}}\\ -{\rm{\nabla }}{\rm{\nabla }}V & -2\tilde{\omega }\end{array}\right].\end{eqnarray} \tag{ 32 }$$

In Equation (32), ∇∇V indicates the gradient matrix of the potential V (a detailed expression of the matrix can be found in Yu 2016), and $\tilde{\omega }$ indicates the skew-symmetric matrix of angular velocity ${\boldsymbol{\omega }}$. ${\rm{O}}$ and I indicate the 3 × 3 zero matrix and the 3 × 3 identity matrix, respectively. Note that in our normalized formulation, ${\boldsymbol{\omega }}$ is a constant vector, thus Equations (31) and (32) show that the Jacobian matrix is totally determined by the gravitational potential of the asteroid, which has a complex geometry due to the asymmetric mass distribution. Furthermore, the variation of the test volume at time t depends on a 3 × 3 block of the Jacobian matrix:

$$\begin{eqnarray}&&\mathrm{volume}({w}^{t})={\int }_{{w}^{t}}d{{\boldsymbol{r}}}^{t}={\int }_{w}\left|{\rm{\det }}\displaystyle \frac{\partial {{\boldsymbol{r}}}^{t}}{\partial {{\boldsymbol{r}}}^{0}}\right|d{{\boldsymbol{r}}}^{0}.\end{eqnarray} \tag{ 33 }$$

Combining Equation (30), the compression/expansion rate of the test volume w yields

$$\begin{eqnarray}&&\chi =\displaystyle \frac{\mathrm{volume}({w}^{t})}{\mathrm{volume}(w)}=| {\rm{\det }}({\{{J}_{{ij}}\}}_{i=1..3,j=1..3}^{t})| .\end{eqnarray} \tag{ 34 }$$

χ is a function of time that provides a measurement to the phase flow close to the equilibrium point. It also gives a mathematical explanation to the results of the second group of simulations. We calculated the compression/expansion rates at the equilibrium point on branch D–D' for cases I and II and compared the curves of χ(t) before and after the shifts of the spin rates (Figure 10). The curves are consistent with the simulation results: the volume near a stable equilibrium point (κ = 2.55, 2.44 and 1.80) oscillates between 0 (fully compressed) and the maximum value (fully expanded), and the maximum value shows a positive correlation with the stability indicator Ξ; i.e., the amplitude of χ decreases as the stability of the equilibrium point increases. After κ drops below the transition point (κ = 1.74 in case II), the amplitude of the oscillation becomes divergent, and the divergence rate increases rapidly with time.

Zoom In Zoom Out Reset image size

The quasi-periodic oscillation of the volume occupied by the cluster hints at the fact that particles moving around a stable equilibrium point experience frequent extrusion, which accordingly increases the collision probability among the particles when the volume is compressed. It therefore leads to an effect that the collision growth of grains near such points may be efficient. If the shed regolith occasionally contributes to a high volume number density of particles orbiting near the stable equilibrium point, we conjecture that the nearby dynamical environment favors the formation of temporary satellites. This is a qualitative understanding based on our analysis and simulations using our "tracer" model. A study of the growth of satellites due to mass shedding would rely on understanding the collision physics, which is crucial for the quantitative understanding of satellite formation. Future work will be performed to explore the role of collisions and its influence on the characteristics of a debris cloud.