The Formation of Bilobate Comet Shapes through Sublimative Torques

Recent spacecraft and radar observations have found that ~70 percent of short-period comet nuclei, mostly Jupiter-family comets (JFCs), have bilobate shapes (two masses connected by a narrow neck). This is in stark contrast to the shapes of asteroids of similar sizes, of which ~14% are bilobate. This suggests that a process or mechanism unique to comets is producing these shapes. Here we show that the bilobate shapes of JFC nuclei are a natural byproduct of sublimative activity during their dynamical migration from their trans-Neptunian reservoir, through the Centaur population, and into the Jupiter family. We model the torques resulting from volatile sublimation during this dynamical migration and find that they tend to spin up these nuclei to disruption. Once disrupted, the rubble pile-like material properties of comet nuclei (tensile strengths of ~1-10 Pa and internal friction angles of ~35$^\circ$) cause them to reform as bilobate objects. We find that JFCs likely experienced rotational disruption events prior to entering the Jupiter family, which could explain the prevalence of bilobate shapes. These results suggest that the bilobate shapes of observed comets developed recently in their history (within the past ~1-10 Myr), rather than during solar system formation or collisions during planet migration and residency in the trans-Neptunian population.

suggesting that either material properties and/or physical processes unique to comets are favoring bilobate shape formation. Comets 8P/Tuttle and 1P/Halley are not classified as JFCs, but rather as Halley-family comets: highly thermally and dynamically evolved isotropic comets. 8P/Tuttle was observed as bilobate by the Arecibo radio telescope; the rest of these nuclei were observed by spacecraft. While detailed three-dimensional shape models do not exist for 19P/Borrelly and 1P/Halley, spacecraft observations suggest that they have bilobate shapes. Comets 9P/Tempel 1 and 81P/Wild are JFCs that are not observed to be bilobate (each has an oblate nucleus that is not considered elongated).
Recent work has shown that re-accretion following either catastrophic collisional (Schwartz et al. 2018) (Jutzi and Benz 2017) or rotational disruption (Sánchez and Scheeres 2016; can produce bilobate objects. The collisional pathway to bilobate shapes is well studied for low velocity (up to a few m/s) (Jutzi and Asphaug 2015), sub-catastrophic (~100s m/s) (Jutzi and Benz 2017), and catastrophic impacts (Schwartz et al. 2018). Remarkably, catastrophic collisions only result in minimal thermal processing of cometary volatiles (Schwartz et al. 2018), allowing impacts to produce these objects under the right conditions. Nevertheless, comets with diameters smaller than 4 km should experience a catastrophic impact (Morbidelli & Nesvorny 2020) with typical impact speeds of 2-4 km/s (Morbidelli & Rickman 2015). However, such high impact speeds struggle to create the larger or more highly elongated bilobate comet nuclei observed by spacecraft (Schwartz et al. 2018), such as 19P/Borrelly or 103P/Hartley 2.
Alternatively, studies have found that the cold classical Kuiper Belt Object (KBO) (486958) Arrokoth obtained its shape through the collapse of a binary system, which merged to form one object (McKinnon et al. 2020). However, the orbital characteristics of the JFC reservoir population are markedly different from those of the cold classical KBOs, which have not experienced planetary encounters (Delsanti & Jewitt 2006). In contrast, the JFCs are thought to originate from the Scattered Disc population (Duncan & Levison 1997;Sarid et al. 2019) and the "stirred" parts of the classical Kuiper Belt (Sarid et al. 2019), which have experienced significant collisional evolution and encounters with the giant planets (Morbidelli & Nesvorny 2020) that can destroy such binary systems. As a result, the Scattered Disk population (the JFC reservoir population), is less likely to retain many similar binary systems, suggesting that some other, enigmatic process must be fundamentally altering the structure, shape, and surfaces of Jupiter-family comets.
Sánchez and Scheeres found that, after rotationally disrupting, small bodies with comet-like strength (Steckloff et al. 2015;Attree et al. 2018) and internal friction angles (Groussin et al. 2015;Steckloff and Samarasinha 2018) tend to reform into bilobate shapes (Sánchez and Scheeres 2016;. Furthermore, once a comet becomes bilobate, further rotational spin up merely separates the lobes, which eventually reaccrete, preserving the bilobate shape (Hirabayashi et al. 2016). Finally, rotational disruption produces fragments with relative velocities small enough that the overwhelming majority of the disrupted body reaccretes into a single body (Sánchez and Scheeres 2016;. This prevents rotational disruptions from producing significant numbers of fragments, consistent with observed JFC and TNO size-frequency distributions. The different dependences of rotational inertia and sublimative torque on the size of the body makes rotational disruptions highly size-dependent (Steckloff and Jacobson 2016). Thus, smaller bodies spin up more rapidly, while larger bodies spin up more slowly. This leads to a critical radius, above which objects are simply too large to spin up to disruption over the dynamical lifetime of a Centaur.
Additionally, not all rotational disruptions necessarily produce bilobate shapes. Bilobate shapeformation requires that cohesion influence an object's internal stress state and disruption dynamics (Sánchez and Scheeres 2016;. However, gravity dominates the cohesive forces of nuclei with radii greater than ~10 km (Sánchez and Scheeres 2014;Kokotanekova et al. 2017) favoring single masses rather than bilobate shapes. Nevertheless, the interior nucleus stress state does not scale analytically in the ~1-10 km transitional size regime (between cohesion-dominated and gravity-dominated dynamics; (Sánchez and Scheeres 2020).

Methods
To understand the limits and constraints of sublimative torques on inducing Centaurs to disrupt (which may then reaccrete/deform into bilobate shapes), we studied how sublimative torques affect these objects during their dynamical migration into the JFC population. We first conducted numerical simulations of Centaur dynamics to obtain a representative sample of the dynamical pathways that JFCs take through the Centaur region. We then computed the sublimative torques experienced by Centaurs along these dynamical pathways, to understand the effects of rotational disruption statistically on the Centaur population as a whole.

Dynamical Evolution Modeling
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, non-gravitational 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ý and Morbidelli 2012;Nesvorný, Vokrouhlický, and Roig 2016;Nesvorný et al. 2017). Our numerical integrator is the swift_rmvs4 code from the Swift N-body integration package (Levison and Duncan 1994), that 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ý and 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 6 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 7 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-year 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 Torque Modeling
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 SYORP sublimative torque model (Steckloff and Jacobson 2016). The SYORP model uses the 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 and 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: (1 − ) where is the Bond albedo of the material, 1678$ is the Sun's luminosity, 9 is the heliocentric distance, )67 is the molar mass of the sublimating volatile, is the ideal gas constant, is the position of each sublimating area relative to the object's subsolar point, and $,-and $,-are a laboratory-determine pressure and temperature on the sublimative phase-change curve (Steckloff et al. 2015;Steckloff and Jacobson 2016). We chose CO and CO2 as representative volatiles which undergo sublimation in the Centaur region. CO sublimation is active throughout the Centaur region, while CO2 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 H2O 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 H2O from the remainder of our analyses.
The SYORP model combines this sublimation pressure with a parameter : (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 and Jacobson 2016;Steckloff and Samarasinha 2018). Thus, the angular accelerations resulting from sublimative torques are described by: where is the bulk density of the object (assumed to be uniform), is the effective radius of the object, and : 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 -4 to 10 -3 (Steckloff and Jacobson 2016), more recent work computed the SYORP coefficients of JFC nuclei to lie in the range ~10 -5 to 10 -4 (Steckloff et al. 2020). Nevertheless, we consider SYORP coefficients that span the entire range (10 -5 to 10 -3 ) 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-year 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-year 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. For our simulations, we used carbon monoxide and carbon dioxide as volatiles being sublimated from the surface of Centaurs. To first order, a molecule's volatility is determined by and inversely related to its latent heat of sublimation. Carbon monoxide is around 3 times more volatile than carbon dioxide by this metric; carbon dioxide is twice as volatile as water. These are the most relevant volatiles sublimating in the Centaur region. is the volatile's latent heat of sublimation and !"# is its molar mass. $%& is an experimentally determined reference vapor pressure at reference temperature $%& .
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 of 100 years, by the angular acceleration at that distance): 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 non-rotating 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 accrued angular velocity for this average orbit is shown as a function of object radius in Figure 1.
We compared this angular velocity to critical angular velocities for breakup in both strength-and gravity-dominated regimes. In the strength-dominated regime, the critical angular velocity C$A= is given by: where is the radius of the object, = 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 and Scheeres 2014). In the gravity-dominated regime, the critical angular velocity C$A= is given by: where is the universal gravitational constant (Pravec and 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 : . 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 : .

Results
We find that sublimation-driven spin-up is capable of disrupting typical JFCs into bilobate shapes during their migration through the Centaur region. We fitted the distribution of critical sizes to Gaussian distributions, to obtain a statistical description of sizes below which objects rotationally disrupt during transit through the Centaur region (see Figure 4 for results for bulk densities of = 500 kg/m 3 ). Mean critical radii based on the fitted Gaussian for each histogram range from 7.2 km for CO2 with : =10 -5 to 134.9 km for CO with : =10 -3 . As expected, the high relative volatility of CO results in larger critical radii. Large critical size outliers appear to skew the distribution a bit. Nevertheless, we find that the distribution of critical sizes is typically much larger than the effective radii of JFCs with known shapes.
We also used our critical size computations to compute the probability of disruption in the Centaur region for each combination of parameters as a function of effective nucleus size, and found that objects with sizes typical for JFCs (on the order of ~1 km in radius) have high probabilities of disrupting in the Centaur region (Figure 3). Even for the weakest combination of parameters (CO2driven activity, : = 10 /F , high density of = 700 kg/m 3 ), the probability of disruption is above 75% for all JFC nuclei with known shapes. For more potent combinations of parameters (COdriven activity, : = 10 /" , low density of = 300 kg/m 3 ), this rises to above 95% for such nuclei. Since such spin up tends to cause objects with JFC-like strength and internal friction angles to disrupt/deform into bilobate shapes (Sánchez and Scheeres 2016), these results suggest that rotational disruption of nuclei in the Centaur region may explain the prevalence of bilobate JFCs. For CO, we find a critical radius between 10 km (for ' = 10 () and body bulk density 700 kg/m 3 ) and 190.2 km (for ' = 10 (* and body bulk density 300 kg/m 3 ). For CO, we find a critical radius between 5.4 km (for ' = 10 () and body bulk density 700 kg/m 3 ) and 100.7 km (for ' = 10 (* and body bulk density 300 kg/m 3 ). Each area corresponds to the upper limit on radius that can be disrupted for a range of densities and the indicated ' : the low edge corresponds to a body bulk density of 700 kg/m 3 , the line through the area corresponds to a body bulk density of 500 kg/m 3 , and the upper edge corresponds to a body bulk density of 300 kg/m 3 . Objects with radii smaller than these values may spin up during their dynamical evolution to angular velocities large enough to cause breakup; after breakup, conditions are favorable for such bodies to reform in bilobate shapes.

Figure 3: Probability of disruption as a function of object radius for CO/CO2 sublimation.
The probability of disruption is calculated statistically from our data set of Centaur orbits by comparing accrued angular velocity for each with the critical values as defined above and in Figure 1. We plot known effective radii of cometary nuclei on this result to show that our computed probabilities of disruption are high for even the weakest combinations of sublimating volatile and SYORP coefficient ' . CO sublimates over the full Centaur region, while CO2 sublimation is turned on within around 13 AU. Figure 2a shows probabilities for CO, while Figure 2b shows probabilities for CO2. Each area corresponds to the probability of disruption as a function of radius for the indicated ' and a range of densities: the low edge corresponds to a body bulk density of 700 kg/m 3 , the line through the area corresponds to a body bulk density of 500 kg/m 3 , and the upper edge corresponds to a body bulk density of 300 kg/m 3 .

Discussion
Our results suggest that the TNO reservoir of JFC progenitors may have a shape distribution that differs significantly from the JFCs. Small non-bilobate TNOs are likely to experience changes that reshape them into bilobate objects in the Centaur region. Nevertheless, were an SDO to already have a bilobate shape, such a shape would likely be preserved during rotational evolution (Hirabayashi et al. 2016), and would survive into the JFC population. Thus, we would only expect the JFC population to have significantly more bilobate objects than their trans-Neptunian reservoir if the population of small objects in this reservoir were not already dominated by bilobate shapes. Conversely, our results suggest that large objects are unlikely to change shape in the Centaur region. Thus, the shapes of the largest Centaurs are likely representative of the shape distribution  Figures 1 and 2). Critical radius is defined as the largest body radius for which the accrued angular velocity over an evolution is sufficient to spin the body to disruption (assuming gravitationally bound bodies). Mean critical radii range from 6.1 km for CO2 with ' = 10 -5 to 114.1 km for CO with ' = 10 -3 . Each distribution is calculated from the same set of 692 orbital pathways, as described above; the distributions are therefore not independent but scaled corresponding to the combination of volatile and ' . of large objects in the JFC reservoir population. However, known JFCs have radii of a few km (Snodgrass et al. 2011); no centaurs or TNOs are known with such small sizes due to observational limitations. Therefore, the effect of this process, and the connection between the TNO, Centaur, and JFC size-frequency distributions is presently unknown.
Previous studies have found that close encounters with planets (such as with giant planets during migration through the Centaur region) could produce tidal forces that elongate objects, potentially into binaries or bilobate shapes (Richardson et al. 2002;Walsh 2018), or even disrupt them into fragments as happened with comet Shoemaker-Levy 9 (Asphaug and Benz 1996). Such dramatic tidal evolution requires close encounters on order the Roche radius or less (a few planetary radii), which happen but are extremely rare for any Centaur. Famously, such an encounter occurred between comet Shoemaker-Levy 9 and Jupiter in 1992 (Asphaug and Benz 1996). To explore the importance of this potential alternative mechanism to forming bilobate shapes in the Centaur region, we used the SyMBA symplectic N-body integrator (Duncan et al. 1998), modified to handle the typical close encounter distances between Centaurs and giant planets. We evolved the orbits of 144 known Centaurs for 10 Myr with a 0.5yr time step and found that Centaurs seldom encounter giant planets as distances less than ~100 planet radii ( Figure 5). Thus, tidal disruptions contribute negligibly to the creation of bilobate comet nuclei.
Furthermore, our results suggest that any typical JFC has a very high likelihood of disrupting into a bilobate nucleus; this mechanism may even work too well, given that not all JFCs (such as Figure 5: Cumulative giant planet encounter distance frequency distribution. We plot the close encounter distance of known centaurs with the giant planets over 10 Myr. Even encounters closer than ~100 planet radii, which is much greater than the encounter distances required to create bilobate shapes through tidal effects (a few planet radii), are rare. 9P/Tempel 1 and 81P/Wild 2) are bilobate. This begs the question: why did these nuclei fail to become bilobate?

9P/Tempel 1, 81P/Wild 2, and Non-bilobate Nuclei
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 and 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: For a tensile strength on the order of ~10 Pa and a density of ~500 kg/m 3 , 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;Richardson et al. 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 non-zero 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 and 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, non-bilobate 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 and 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 3 . We found that such a typical Centaur would lose 340 km or 7100 km to CO2 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).
We can estimate what these actual size limits are using a simple calculation. In our model, angular acceleration depends on 1 * ⁄ (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 : = 10 /" and : = 10 /F would need to have radius smaller than ~1 km and ~100 m, respectively, to be disrupted, while an object sublimating CO2 would need to have radius smaller than ~200 m and ~20 m (respectively for the above : values). This suggests that CO2 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 : values may also produce an overabundance of bilobate nuclei. The range of empirical : values (10 -5 -10 -4 ) was computed from observations of rotation state changes in JFC nuclei (Steckloff et al. 2020), either from lightcurve 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 : parameter, with higher levels of activity favoring greater : values, potentially biasing this sample toward larger : values. In such a case, : values for a representative population may plausibly be lower than 10 -5 -10 -4 . As angular acceleration from sublimation scales linearly with : , a lower average : would reduce the critical radius limit for disruption that we calculated. The probabilities of disruption for non-bilobate JFCs could in reality be sufficiently low to explain their existence.
Additionally, our : values may overestimate the concentration of supervolatile species across the nucleus surface (the volatile abundance/distribution component of : ). Presently, all empirically derived : values are for H2O sublimation in the inner solar system, with the notable exception of 103P/Hartley 2, whose activity is dominated by CO2 sublimation (Steckloff et al. 2020). Although volatile abundances within Centaurs are poorly constrained, CO and CO2 are nevertheless typically significantly less abundant than H2O 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 CO2 production rates, which may be as much as an order of magnitude smaller than that of H2O. 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 : .
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 cancelations (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.22x10 -5 (Steckloff et al. 2020), 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. 2020); 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.

Split Comets
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 and 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 and 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 Comets
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 1 . 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 and 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 4 -10 5 years (Steinberg and Sari 2011). However, sublimative forces are typically ~10 5 times stronger than radiative forces when sublimation is vigorous (Steckloff and Jacobson 2016), causing vigorously sublimating systems to evolve ~10 3 -10 4 times faster, once accounting for the relative weakness of SYORP coefficients relative to their YORP coefficient analogues (Steckloff et al. 2020) 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-year 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.

Conclusions
Unlike the near-Earth asteroid population, the majority of short period comet nuclei have bilobate shapes, suggesting that material properties and/or physical processes unique to comets are favoring bilobate shape formation. Previously proposed mechanisms to explain this discrepancy involve either improbably low-velocity collisions between similarly sized objects, or high-velocity impacts that are inconsistent with the size-frequency distribution of the JFC reservoir population. We explore an alternative mechanism, in which sublimative torques rotationally disrupt JFC nuclei through the Centaur region from their trans-Neptunian reservoir population. Given the material properties of comet nuclei, such disruption events are likely to result in nuclei reforming or deforming into bilobate shapes. We find that this mechanism is extremely efficient at producing bilobate shapes among the population of JFCs with known sizes. Even our worst-case scenario produces a probability of disruption for the largest JFCs with known shapes of greater than 75%, with our best-case scenarios disrupting all JFCs with near certainty. Nevertheless, there is an unknown cutoff point in body size beyond which rotational disruption is unlikely to result in bilobate shapes; this size is presently unknown.