The Formation of Bilobate Comet Shapes through Sublimative Torques

Our dynamical modeling effort numerically integrated the evolution of orbits in the TNO population over the age of the solar system, due to gravitational perturbations of the giant planets.  This model neglects sublimative, nongravitational forces acting on Centaurs, which may additionally perturb these orbits. Rigorous testing has verified our TNO population model's accuracy against both the observed population of TNOs and reservoir formation theories (Nesvorný & Morbidelli 2012; Nesvorný et al. 2016, 2017). Our numerical integrator is the swift_rmvs4 code from the Swift N-body integration package (Levison & Duncan 1994), which has been modified to account for the radial migration and damping of planetary orbits in the early solar system via the exponential e-folding timescales (Nesvorný & Morbidelli 2012). Further details of our dynamical model, including a detailed description of the integration method, planet migration, initial orbital distribution of disk planetesimals, and comparison of results with the orbital structure of the TNO population can be found in Nesvorný et al. (2017).

Our numerical integrator tracks the orbital evolution of 10⁶ outer disk planetesimals from the onset of Neptune's migration to the present time. To improve statistics for modern day Centaurs, we produced 100 clones of objects that entered the Centaur region during the past 1 Gyr and repeated their dynamical evolution through the Centaur region. The cloning was accomplished by introducing a small (random) change of the velocity vector of a particle when it first evolves to an orbit with semimajor axis a < 30 au (i.e., enters the Centaur region). We recorded the cloned orbits every 100 years, until the object was either ejected from the solar system or evolved into a JFC orbit. In this manner, we produced a total of ∼10⁷ Centaur orbits, resulting in a statistically significant sample of 346 representative dynamical pathways for JFCs through the Centaur region.

When the object is between 5 and 30 au (i.e., in the Centaur region), we record the semimajor axis, eccentricity, and current heliocentric distance in 100 yr intervals for each dynamical pathway (this interval is ∼4–5 orders of magnitude shorter than the timescale of dynamical evolution through the Centaur region; e.g., di Sisto & Brunini 2007). This step unfortunately introduces an ambiguity into our simulations, as we do not record whether the particle is approaching aphelion or perihelion (or, equivalently, the mean, true, or eccentric anomalies). We therefore consider both cases for each dynamical pathway, which we ultimately find introduces a very small uncertainty into our constraints of rotational disruption.

Sublimative torques arise from the asymmetric emission of sublimating volatiles from the irregular shapes of icy bodies. To compute their effects on the rotation state of Centaurs, we use the parameterized Sublimative Yarkovsky–O'Keefe–Radzievskii–Paddack (SYORP) sublimative torque model (Steckloff & Jacobson 2016). The SYORP model uses the Yarkovsky–O'Keefe–Radzievskii–Paddack (YORP) effect of radiative torques as a template to compute the torques resulting from sublimation. Formulations of the YORP effect compute the magnitude of torques from a shape-dependent parameter and the magnitude of the radiation pressure at zero solar phase angle (Scheeres 2007; Rozitis & Green 2013). Similarly, the SYORP model uses first principles to solve for the mass flux and thermal velocity of sublimating volatiles leaving the surface, and thus the dynamic sublimative pressure that these gases exert on the surface. The SYORP model treats Centaurs as spheres of a single, pure volatile ice with negligibly low albedos (consistent with the observed low bond albedos of JFC nuclei), resulting in a pair of equations that can be solved simultaneously:

$$\begin{eqnarray}&&\left(1-A\right)\displaystyle \frac{{L}_{\mathrm{solar}}}{4\pi {r}_{h}^{2}\lambda }\cos \phi ={\alpha }_{\left(T\right)}\sqrt{\displaystyle \frac{{m}_{\mathrm{mol}}}{2\pi {RT}}}{P}_{{ref}}{e}^{\tfrac{\lambda }{R}\left(\tfrac{1}{{T}_{{ref}}}-\tfrac{1}{T}\right)},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{P}_{\mathrm{sub}}=\displaystyle \frac{2}{3}\left(1-A\right)\displaystyle \frac{{L}_{\mathrm{solar}}}{4\pi {r}_{h}^{2}\lambda }{\left(\displaystyle \frac{8{RT}}{\pi {m}_{\mathrm{mol}}}\right)}^{\tfrac{1}{2}}\cos \phi ,\end{eqnarray} \tag{ 2 }$$

where A is the Bond albedo of the material, L_solar is the Sun's luminosity, r_h is the heliocentric distance, m_mol is the molar mass of the sublimating volatile, R is the ideal gas constant, ϕ is the position of each sublimating area relative to the object's subsolar point, and P_ref and T_ref are a laboratory-determine pressure and temperature on the sublimative phase-change curve (Steckloff et al. 2015; Steckloff & Jacobson 2016). We chose CO and CO₂ as representative volatiles which undergo sublimation in the Centaur region. CO sublimation is active throughout the Centaur region, while CO₂ sublimation is active at heliocentric distances within approximately 13 au. We chose an appropriate set of properties for each volatile used in our model, which are shown in Table 1 below. We also considered the effect of H₂O sublimation but found that it is ultimately too refractory in the Centaur region to contribute to any significant changes in rotational states; we therefore neglect H₂O from the remainder of our analyses.

The SYORP model combines this sublimation pressure with a parameter C_S (the "SYORP coefficient"), which encodes the effects of shape, spin pole orientation/obliquity, depth to sublimation fronts, and volatile-distribution on the magnitude of sublimative torques (Steckloff & Jacobson 2016; Steckloff & Samarasinha 2018). Thus, the angular accelerations resulting from sublimative torques are described by

$$\begin{eqnarray}&&\displaystyle \frac{d\omega }{{dt}}=\displaystyle \frac{3{P}_{\mathrm{sub}}{C}_{S}}{4\pi \rho {r}^{2}},\end{eqnarray} \tag{ 3 }$$

where ρ is the bulk density of the object (assumed to be uniform), r is the effective radius of the object, and C_S is the shape- and distribution-dependent SYORP coefficient. Although small, pristine icy bodies with surface-located sublimation fronts may have SYORP coefficients on the order of ∼10⁻⁴ to 10⁻³ (Steckloff & Jacobson 2016), more recent work computed the SYORP coefficients of JFC nuclei to lie in the range ∼10⁻⁵ to 10⁻⁴ (Steckloff et al. 2021). Nevertheless, we consider SYORP coefficients that span the entire range (10⁻⁵ to 10⁻³) for completeness.

To compute the effects of sublimative torques during Centaur dynamical evolution, we applied this SYORP model to the computed Centaur dynamical pathways as they evolve inward. We used the Centaurs' orbital elements in each 100 yr interval to compute the amount of time that it spent in 1 au heliocentric distance bins in this range using the given orbital elements integrated forward for each 100 yr timestep. For each body, we were able to compute a "best-case" (most spin-up) and "worst-case" (least spin-up) evolution given the uncertainty in orbital elements; we defined "best-case" as approaching perihelion at the start of the integration, and therefore spending more time near the Sun; "worst-case" was likewise defined as moving away from perihelion at the start of the integration and therefore spending less time near the Sun. In this way, we calculated 692 possible dynamical pathways from the 346 simulated pathways provided. Because sublimative changes in a Centaur's rotation state occur over secular (rather than orbital) timescales, this method is equivalent to computing the angular acceleration resulting from sublimative torques (in 1 au intervals) and integrating over an entire orbit to compute the change in angular speed over a single orbit.

We calculated the total accrued angular velocity over an orbital evolution by summing the angular velocities accrued in each 1 au heliocentric distance bin (given by multiplying the time spent at that distance, with minimum dt of 100 years, by the angular acceleration at that distance):

$$\begin{eqnarray}&&{\omega }_{\mathrm{final}}=\sum \displaystyle \frac{d\omega }{{dt}}{dt}.\end{eqnarray} \tag{ 4 }$$

We assumed that rotation state changes are solely the result of sublimative torques, ignoring the possibility of collisions or any other mechanism, and that such torques monotonically increased the angular velocity of each body from a nonrotating starting state. We also assumed that spin state changes were cumulative over the dynamical evolution of a Centaur, and therefore we were interested only in time spent at each heliocentric distance. With this method, we calculated the total change in rotation state (angular speed) for each of the 692 orbital cases as a function of object radius, as well as for a computed average orbital evolution (found by taking the average of the time spent in each 1 au heliocentric distance bin for each of our orbits). The final spin period for this average orbit is shown as a function of object radius in Figure 2.

We compared this angular velocity with critical angular velocities for breakup in both strength- and gravity-dominated regimes. In the strength-dominated regime, the critical angular velocity ω_crit is given by

$$\begin{eqnarray}&&{\omega }_{\mathrm{crit}}=\sqrt{\displaystyle \frac{2{\sigma }_{t}}{\rho {r}^{2}}},\end{eqnarray} \tag{ 5 }$$

where r is the radius of the object, σ_t is the tensile material strength and ρ is the material density. For cometary nuclei, we assumed a tensile strength on the order of a few Pa, consistent with a rubble pile bound by Van der Waals forces (Sánchez & Scheeres 2014). In the gravity-dominated regime, the critical angular velocity ω_crit is given by

$$\begin{eqnarray}&&{\omega }_{\mathrm{crit}}=\sqrt{\displaystyle \frac{4\pi \rho G}{3}}\end{eqnarray} \tag{ 6 }$$

where G is the universal gravitational constant (Pravec & Harris 2000). Note that, unlike the strength-dominated critical angular velocity, the gravity-dominated critical velocity does not depend on the size of the object. Based on these critical angular velocities, we can compute a critical radius along every dynamical pathway for each combination of sublimating volatile and SYORP coefficient C_S . Below this critical radius, an object will disrupt during its migration. We used our model to compute the critical radius of disruption for each orbital pathway in our evolution set, and then computed a probability of disruption and critical size distribution (Figure 3) as a function of size, sublimating volatile, bulk density, and C_S .

Rotational disruption produces bilobate shapes (rather than losses of surface materials) when cohesive forces are strong enough to induce structural failure in the interior of the body, rather than at the surface (Sánchez & Scheeres 2016). As objects increase in size, gravitational forces play an increasing role in binding the object together. By setting Equations (5) and (6) equal to each other we can find the size at which nuclei transition from being predominantly bound by cohesion to predominantly bound by gravity:

$$\begin{eqnarray}&&r=\displaystyle \frac{1}{\rho }\sqrt{\displaystyle \frac{3{\sigma }_{t}}{2\pi G}}.\end{eqnarray} \tag{ 7 }$$

For a tensile strength on the order of ∼10 Pa and a density of ∼500 kg m⁻³, we obtain a radius of ∼500 m. Above this radius, cohesion still helps to bind the nucleus, but plays a decreasing role as nuclei increase in radius. Above some size (perhaps an order of magnitude larger than this transitional size), gravitational forces sufficiently dominate cohesive forces to prevent cohesion from influencing the outcome of rotational disruption. 9P/Tempel 1 and 81P/Wild 2 may be in this category of objects small enough to rotationally disrupt, yet too large to reform/deform into a bilobate shape. Indeed, Tempel 1 seems to exhibit large scale layering throughout its nucleus (Thomas et al. 2007), which could be the result of its disruption.

To form our results, we assumed that all comets that spin up to disruption reform into bilobate shapes. However, previous numerical studies have shown that rotational spin-up of rubble piles does not necessarily produce bilobate shapes (Richardson et al. 2005, 2009). These studies found that a rotationally disrupted aggregate's fate is sensitive to material cohesion: in strengthless cases, objects tend to elongate (Richardson et al. 2005), and in high-strength cases, objects tend to shed material (Richardson et al. 2009). However, in cases of low but nonzero cohesive strength (on the order of ∼100 Pa), this work found that nucleus deformation can occur before failure, which could lead to the formation of binary or bilobate objects (Richardson et al. 2009). Sánchez & Scheeres (2016) found that bilobate formation also depends on internal friction angle and found that the cohesive strengths and internal friction angles that result in bilobate formation are ∼0.6–600 Pa and ∼25°–35°, respectively. Nevertheless, it is plausible that some Centaurs could have material properties that lie outside of this range, which would preclude bilobate shape formation. Therefore, our results may overestimate the probability of bilobate formation.

Alternatively, 9P/Tempel 1 and 81P/Wild 2 could have become bilobate but lost their smaller lobes. Once a nucleus becomes bilobate, rotational spin-up will cause the two lobes to separate (Hirabayashi et al. 2016). The fate of such lobes depends on the relative masses of the two lobes; if their mass ratio is greater than 0.2 (similar-size lobes), the two lobes remain gravitationally bound, and will reaccrete into a bilobate shape. However, if the lobes' mass ratio is less than 0.2, the lobes tend to be gravitationally unbound; the lobes separate permanently and the comet splits (Hirabayashi et al. 2016). If either 9P/Tempel 1 or 81P/Wild 2 had a sufficiently small lobe, it may have split from the larger lobe, leaving a single, nonbilobate lobe comprising the nucleus. However, the nucleus could then spin up again to bilobate formation; a process that could plausibly repeat until a stable bilobate nucleus with lobes comparable in size (mass ratio greater than 0.2) formed. Nuclei close to the critical disruption radius would be less likely to have sufficient time to experience multiple disruption cycles; if such nuclei only produced small lobes less than 20% the mass of the larger lobe and continued to escape, then the prevalence of bilobate shapes should generally decrease with increasing size near the disruption limit. Nevertheless, both of these comets are much smaller than the critical disruption limit, and previous modeling work tends to show that roughly equal lobes would be created (Sánchez & Scheeres 2016). We therefore consider this scenario to be less likely.

Assumptions made in our model may produce errors in the calculated sizes at which bilobate objects can form. For example, real bodies are finite in size and may be depleted in volatiles before they are able to disrupt. However, our model assumes that objects have an infinite supply of volatiles. To estimate this effect, we calculated the depth of ice that would be lost from a body undergoing an average Centaur evolution (as we used to create Figure 2), assuming a density of 500 kg m⁻³. We found that such a typical Centaur would lose 340 km or 7100 km to CO₂ or CO sublimation, respectively, before entering the Jupiter family. No known Centaur is this large, and bodies with this amount of ice would not be able to spin up to disruption, suggesting that the actual size limits for bilobate formation are smaller than our calculated limits as given in Figure 4 (see the mean of each distribution).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can estimate what these actual size limits are using a simple calculation. In our model, angular acceleration depends on 1/r² (as in Equation (3)), so a body with half of the radius of another would require a fourth of the amount of ice to spin to disruption. Applying this same scaling to our modeled size limits, an object sublimating CO with SYORP coefficients C_S  = 10⁻⁴ and C_S  = 10⁻⁵ would need to have radius smaller than ∼1 km and ∼100 m, respectively, to be disrupted, while an object sublimating CO₂ would need to have radius smaller than ∼200 m and ∼20 m (respectively for the above C_S values). This suggests that CO₂ is much more important than CO in spinning up JFCs (at the observed size range of ∼1 km radius). Additionally, the distribution of bilobate JFCs is likely sensitive to the distribution of SYORP coefficients among Centaurs.

Errors and biases in estimates of C_S values may also produce an overabundance of bilobate nuclei. The range of empirical C_S values (10⁻⁵–10⁻⁴) was computed from observations of rotation state changes in JFC nuclei (Steckloff et al. 2021), either from light-curve studies or spacecraft observations. However, such observations are inherently biased toward more active comets, which are more visible from Earth. However, the level and distribution of sublimative activity across the nucleus surface is folded into the C_S parameter, with higher levels of activity favoring greater C_S values, potentially biasing this sample toward larger C_S values. In such a case, C_S values for a representative population may plausibly be lower than 10⁻⁵–10⁻⁴. As angular acceleration from sublimation scales linearly with C_S , a lower average C_S would reduce the critical radius limit for disruption that we calculated. The probabilities of disruption for nonbilobate JFCs could in reality be sufficiently low to explain their existence.

Additionally, our C_S values may overestimate the concentration of supervolatile species across the nucleus surface (the volatile abundance/distribution component of C_S ). Presently, all empirically derived C_S values are for H₂O sublimation in the inner solar system, with the notable exception of 103P/Hartley 2, whose activity is dominated by CO₂ sublimation (Steckloff et al. 2021). Although volatile abundances within Centaurs are poorly constrained, CO and CO₂ are nevertheless typically significantly less abundant than H₂O in cometary comae (Bockelée-Morvan et al. 2004), resulting in lower volatile production rates which drive sublimative torques. If we assume that these coma abundances reflect surface compositions, our model may overestimate CO or CO₂ production rates, which may be as much as an order of magnitude smaller than that of H₂O. This would ultimately result in smaller critical radii for disruption than computed. However, these reduced abundances would have a surprisingly small effect on the net sublimative torque experienced by comets. Most comet outgassing produces torques that stochastically cancel out; reducing volatile production rates produces a corresponding change in the measured net torque of only ∼3%–20% (Samarasinha & Mueller 2013; Mueller & Samarasinha 2018). Thus, this error is unlikely to produce significant errors in critical radii for disruption, which are already dominated by uncertainties in the SYORP parameter C_S .

Finally, our assumption of monotonic changes in spin rate may be overestimating actual spin rate changes. The analogous radiative YORP torques are thought to be stochastic or varying in magnitude and direction over time (Statler 2009; Cotto-Figueroa et al. 2015), causing such monotonic assumptions to overestimate maximum spin rates (specifically, angular speeds) by a small factor on the order of unity (Cotto-Figueroa et al. 2015). Sublimative torques likely follow a similar pattern, as rotation changes may activate new areas of sublimative activity (Steckloff & Samarasinha 2018) that change the magnitude and direction of sublimative torques. These torque changes may be somewhat self-limiting due to stochastic torque cancellations (Samarasinha & Mueller 2013). Nevertheless, if sublimative torque changes behave similarly to radiative torques, the critical sizes for each dynamical pathway (and thus the size of disrupted nuclei) may be overestimated by a small factor (on the order of unity). Given that 9P/Tempel 1 has an SYORP coefficient of only 1.22 × 10⁻⁵ (Steckloff et al. 2021), such a shift could make it unlikely that Tempel 1 would have ever spun up to disruption. Similarly, comet 19P/Borrelly has a large SYORP coefficient (Steckloff et al. 2021); if 81P/Wild 2 has a smaller SYORP coefficient, then this time varying effect could explain both why Borrelly has a bilobate shape and Wild 2 does not.

Rotational disruption, perhaps driven by sublimative torques, has long been proposed as a mechanism for comet splitting (e.g., Boehnhardt 2000). For example, comet 3D/Biela split into two distinct objects in 1846 before disappearing, presumed disintegrated (Jenniskens & Vaubaillon 2007). If Biela was originally a bilobate comet with lobes having a mass ratio of less than 0.2, rotational disruption could have split Biela into a gravitationally unbound system, causing its lobes to separate into two visually distinct objects. A few such systems exist; recent observations have revealed another split comet candidate in comets 252P and BA14 (Li et al. 2017). More generally, Chen and Jewitt predict a splitting rate of cometary nuclei of 0.01 per comet per year (Chen & Jewitt 1994). Sublimation-driven rotational spin-up may explain this rate of comet bursting; however, such an investigation is outside the scope of this work.

Binary asteroids are thought to form through rotational disruption, in which an asteroid rotates rapidly enough for material to be rotationally shed from the equator before collecting into a small moon. The formation or rotational evolution of comet nuclei should similarly produce binary systems, if only temporarily. However, no comets or Centaurs are known to be binary systems. ⁷ One reason for this may be the differences in the strengths of the forces that evolve these systems.

Small asteroid systems evolve primarily though radiative forces (e.g., YORP and binary-YORP or "BYORP"), which can cause the secondary to either spiral in or leave the system (Cuk & Burns 2005). A typical 1 km radius asteroid with a 0.2 km radius secondary ("moon") at 1au would experience such evolution over a timescale of ∼10⁴–10⁵ yr (Steinberg & Sari 2011). However, sublimative forces are typically ∼10⁵ times stronger than radiative forces when sublimation is vigorous (Steckloff & Jacobson 2016), causing vigorously sublimating systems to evolve ∼10³–10⁴ times faster, once accounting for the relative weakness of SYORP coefficients relative to their YORP coefficient analogs (Steckloff et al. 2021) and the relatively low density of cometary bodies relative to asteroids. Thus, sublimative forces would cause a comparable binary cometary system with a 1 km primary at 0.2 km secondary at 1 au to evolve to either reaccrete into a binary system or separate over a ∼1–100 yr timescale. Such timescales are incredibly short, resulting in a significantly lower probability that a binary comet system would be observed relative to a binary asteroid system. Furthermore, the comae of vigorously sublimating objects could further obscure any binary system, rendering it unlikely that any binary cometary system would have ever been detected.