[orcid=0000-0002-0140-4475]

\cormark

[1]

\fnmark

[1]

1]organization=Dept. of Astronomy, University of Michigan, city=Ann Arbor, postcode=48109, state=MI, country=USA

[orcid=0000-0002-1717-2226] 2]organization=The Planetary Science Institute, city=Tucson, state=AZ, country=USA

3]organization=Dept. of Aerospace Engineering and Engineering Mechanics, The University of Texas at Austin, city=Austin, state=TX, country=USA

[orcid=0000-0002-0726-6480] \fnmark[2] 4]organization=Dept. of Astronomy and Carl Sagan Institute, Cornell University, addressline=122 Sciences Drive, city=Ithaca, postcode=14853, state=NY, country=USA

[orcid=0000-0003-0774-884X] 5]organization=Jet Propulsion Laboratory, California Institute of Technology, addressline=4800 Oak Grove Dr., city=Pasadena, postcode=91109, state=CA, country=USA

[orcid=0000-0001-8573-7412] 6]organization=Southwest Research Institute, addressline=1050 Walnut Street, Suite 300, city=Boulder, postcode=80302, state=CO, country=USA

[orcid=0000-0002-6034-5452] 7]organization=Institute of Astronomy, Charles University, addressline=V Holešovičkách 2, CZ-18000 , city=Prague, citysep=, postcode=8, country=Czech Republic

[orcid=0000-0002-4547-4301]

[orcid=0000-0001-7895-8209] 8]organization=ESA NEO Coordination Centre, addressline=Largo Galileo Galilei 1, I-00044 Frascati (RM), country=Italy

\cortext

[cor1]Corresponding author

\fntext

[fn1]Fannie and John Hertz Foundation Fellow \fntext[fn2]NSF Astronomy and Astrophysics Postdoctoral Fellow

The Dynamical Origins of the Dark Comets and a Proposed Evolutionary Track

Aster G. Taylor agtaylor@umich.edu [    Jordan K. Steckloff [ [    Darryl Z. Seligman [    Davide Farnocchia [    Luke Dones [    David Vokrouhlický [    David Nesvorný    Marco Micheli [
Abstract

So-called ‘dark comets’ are small, morphologically inactive near-Earth objects (NEOs) that exhibit nongravitational accelerations inconsistent with radiative effects. These objects exhibit short rotational periods (minutes to hours), where measured. We find that the strengths required to prevent catastrophic disintegration are consistent with those measured in cometary nuclei and expected in rubble pile objects. We hypothesize that these dark comets are the end result of a rotational fragmentation cascade, which is consistent with their measured physical properties. We calculate the predicted size-frequency distribution for objects evolving under this model. Using dynamical simulations, we further demonstrate that the majority of these bodies originated from the ν6subscript𝜈6\nu_{6}italic_ν start_POSTSUBSCRIPT 6 end_POSTSUBSCRIPT resonance, implying the existence of volatiles in the current inner main belt. Moreover, one of the dark comets, (523599) 2003 RM, likely originated from the outer main belt, although a JFC origin is also plausible. These results provide strong evidence that volatiles from a reservoir in the inner main belt are present in the near-Earth environment.

keywords:
Asteroids (72) \sepCelestial mechanics (211) \sepComets (280) \sepNear-Earth Objects (1092) \sepOrbital evolution (1178)

1 Introduction

Classically, small solar-system bodies are divided into two broad populations — asteroids and comets. In general, asteroids are on short-period orbits and lack visibly sublimating volatiles. The previously-measured nongravitational accelerations on asteroids are generally attributed to radiation pressure (Vokrouhlický & Milani, 2000) or, in most cases, the Yarkovsky effect (Farnocchia et al., 2013; Vokrouhlický et al., 2015; Greenberg et al., 2020). Most asteroids have typical densities of ρ10003000similar-to𝜌10003000\rho\sim 1000-3000italic_ρ ∼ 1000 - 3000 kg/m3 (Carry, 2012) and reside in the main belt between Mars and Jupiter or in the Trojan populations co-orbiting with Jupiter. Due to the Yarkovsky effect, some asteroids drift into resonances such as the ν6subscript𝜈6\nu_{6}italic_ν start_POSTSUBSCRIPT 6 end_POSTSUBSCRIPT secular resonance and the 3:1 mean-motion resonance with Jupiter (Gladman et al., 1997; Migliorini et al., 1998; Granvik et al., 2017). These resonances excite the asteroids’ eccentricities, causing many of them to dynamically evolve into near-Earth objects (NEOs) with perihelia of q1.3𝑞1.3q\leq 1.3italic_q ≤ 1.3 au (Nesvorný et al., 2017; Seligman et al., 2021; Nesvorný et al., 2023). However, the near-Earth environment is dynamically unstable and NEOs are continuously removed via scattering into the Sun, scattering onto Jupiter-crossing orbits, and planetary collisions (Farinella et al., 1994; Nesvorný et al., 2023).

On the other hand, most comets are thought to originate from either the Oort Cloud or the Trans-Neptunian Object (TNO) population, where low ambient temperatures allow volatiles to exist for extended periods (Volk & Malhotra, 2008; Nesvorný, 2018; Lisse et al., 2022). The presence of relatively light volatiles and high nuclear porosity in comets results in consistently lower bulk densities than asteroids, typically on the order of ρ500similar-to𝜌500\rho\sim 500italic_ρ ∼ 500 kg/m3 (Knight et al., 2023). Comets typically exhibit their eponymous comae and nongravitational accelerations, which arise from the outgassing of sublimating volatiles upon entering the inner solar system (Whipple, 1950, 1951; Yeomans et al., 2004). Recoil forces from outgassing tend to spin up these nuclei, resulting in rotational fragmentation as centrifugal forces overcome the binding forces (Steckloff & Jacobson, 2016; Jewitt, 2021, 2022). This effect, in combination with mass loss due to outgassing, generally leads to “cometary fading”, the destruction of comets on relatively short timescales (Wiegert & Tremaine, 1999; Brasser & Wang, 2015).

It is expected that most short-period comets, particularly Jupiter Family and Encke-type comets, originate from the TNO reservoir (Levison & Duncan, 1997; Nesvorný, 2018). Analogously to the main belt, interactions with the giant planets (particularly Neptune) can destabilize and scatter the orbits of Kuiper Belt objects (KBOs) inwards (Nesvorný et al., 2017). These KBOs then transition to the giant planet region between Jupiter and Neptune. At this point, these objects are classified as Centaurs, named after the first object recognized to be in this class, (2060) Chiron (Kowal et al., 1979). While definitions vary, Centaurs are generally understood to be objects with semimajor axes of a5𝑎5a\geq 5italic_a ≥ 5 au and perihelia of q30𝑞30q\leq 30italic_q ≤ 30 au (Jewitt, 2009).

Centaur orbits are unstable and chaotically evolve due to their proximity to the giant planets. About one-third of Centaurs are scattered inwards via interactions with Jupiter and become JFCs after 110similar-toabsent110\sim 1-10∼ 1 - 10 Myr as a Centaur (Tiscareno & Malhotra, 2003; Di Sisto & Brunini, 2007; Seligman et al., 2021). JFCs are then destroyed by repeated perihelion passages (Nesvorný, 2018), ejected from the solar system, or dynamically detached from Jupiter via migration into the near-Earth environment (Grav et al., 2023; Nesvorný et al., 2023).

Recent developments have demonstrated that a continuum of objects exists between the two extremes of asteroids and comets, such as the active asteroids and inactive comets. Active asteroids are objects on asteroid-like orbits with observable cometary activity (Jewitt, 2012; Hsieh, 2017). A particular example is the Main Belt Comets (MBCs), which reside within the main belt but exhibit comae (Hsieh & Jewitt, 2006). It is worth mentioning that the MBCs are a subset of active asteroids. MBC activity is due to sublimation of ice, while active asteroids may be active due to rotational instability or impacts (Jewitt, 2012). Only a handful of objects have been identified as MBCs, and targeted searches suggest occurrence rates of <<<1/500 and similar-to\sim1/300 in the outer main belt (Sonnett et al., 2011; Bertini, 2011; Ferellec et al., 2023). See Jewitt & Hsieh (2022) for a recent review.

Inactive comets are on cometary orbits but display no detectable activity. While there are several categories of inactive comets, classified based on their orbits, nearly all are likely extinct or dormant comet nuclei that were not destroyed by outgassing or outgassing-driven processes (Asher et al., 1994; Jewitt, 2005; Licandro et al., 2018). There are also a few known weakly active, tailless “Manx” comets (Meech et al., 2016; Piro et al., 2021; Kwon et al., 2022). It has been suggested that the Manx comets are asteroids that were scattered into the Oort cloud and developed a volatile-rich frost before traveling back toward the Sun (Weissman & Levison, 1997; Shannon et al., 2015; Meech et al., 2016).

Recently, Farnocchia et al. (2023) and Seligman et al. (2023) (building on the work of Chesley et al., 2016) defined “dark comets” as a population of active asteroids that are distinct from these previously identified categories. The dark comets are on near-Earth orbits and exhibit no detectable comae (hence the “dark” moniker), yet experience nongravitational accelerations that are inconsistent with radiative effects. Farnocchia et al. (2023) and Seligman et al. (2023) demonstrated that outgassing is consistent with both the reported accelerations and the nondetections of comae. In addition to their accelerations and absence of comae, dark comets exhibit several other notable properties — they are generally (i) small (on the order of 101001010010-10010 - 100 m), (ii) rapidly rotating (0.0461.990.0461.990.046-1.990.046 - 1.99 h, Seligman et al., 2023), and (iii) accelerating in nonradial directions.

Taylor et al. (2024) demonstrated that if rapid rotation causes diurnal outgassing to be negligible, differential heating across the hemispheres of these objects due to seasonal effects could generate the observed out-of-plane acceleration. This mechanism is consistent with the accelerations measured on all the currently known dark comets except for 2003 RM, which also differs from the other objects in size and orbital parameters. However, this mechanism relies on the rotation rate, which is unmeasured on a majority of the dark comets.

Here, we further investigate the properties and origins of the dark comets. This paper is organized as follows: in Sec. 2, we calculate a preliminary estimate of the population fraction in the current near-Earth environment. In Sec. 3, we calculate the minimum material strength required to hold the dark comets stable against rotational disruption, and compare with the estimated strengths of known comet nuclei. In Sec. 4, we introduce a proposed evolutionary track for these objects, in which rotational splitting and fragmentation results in nuclei consistent with dark comets. In Sec. 5, we calculate the size-frequency distribution that our model predicts for dark comets. In Sec. 6, we discuss the dynamical transfer of objects to the dark comet orbits. Finally, in Sec. 7, we discuss the implications of our results.

2 Approximate Occurrence Rate

Given the long temporal baselines of astrometric data required to measure nongravitational accelerations and the small sizes of the dark comets, it is possible that more as-yet undetected dark comets exist in the inner solar system. In this section, we provide an order-of-magnitude estimate of the population of dark comets, albeit one subject to large uncertainties due to small-number statistics, unknown detection biases, and poorly constrained evolutionary histories.

To estimate the fraction of NEOs that are dark comets, we first estimate the fraction of NEOs that can be ruled out as dark comets. Specifically, we calculate a best-fit nongravitational acceleration for each NEO, including uncertainty. We then find the number of objects with sufficiently small uncertainties in hypothetical nongravitational accelerations to rule out dark comet–like values. The population-level occurrence rate is then approximately the number of known dark comets divided by the total number of objects within this subset. As only 7 dark comets are currently known, this is a first-order estimate.

We estimated nongravitational accelerations for 10 553 NEOs using the least-squares fitting method of Farnocchia et al. (2015). Notably, the same methodology was used to identify the dark comets, thereby cancelling any biases in the estimation of the nongravitational accelerations. The optical astrometric data are obtained from the Minor Planet Center111https://www.minorplanetcenter.net/ and radar measurements from the Solar System Dynamics Group222https://ssd.jpl.nasa.gov/sb/radar.html at JPL. The nongravitational accelerations are modeled as Ai(1 au/r)2subscript𝐴𝑖superscript1 au𝑟2A_{i}\,(\text{1 au}/r)^{2}italic_A start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ( 1 au / italic_r ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT in each of the orbital coordinate directions (radial, transverse, out-of-plane), following the model of Marsden (1969).

Many of these objects will have nongravitational accelerations from radiation pressure or the Yarkovsky effect, independent of any dark comet–like behavior. To account for this, we obtain best-fit lines for the acceleration magnitude of radiation pressure and the Yarkovsky effect, with object data given by the JPL Small-Body Database333https://ssd.jpl.nasa.gov/tools/sbdb_query.html. We assume that radiation pressure and the Yarkovsky effect both exhibit a 1/RN1subscript𝑅N1/R_{\rm N}1 / italic_R start_POSTSUBSCRIPT roman_N end_POSTSUBSCRIPT dependence, where RNsubscript𝑅NR_{\rm N}italic_R start_POSTSUBSCRIPT roman_N end_POSTSUBSCRIPT is the radius of the nucleus and is computed from the absolute magnitude H𝐻Hitalic_H by

RN=(1329 km2p)100.2H.subscript𝑅N1329 km2𝑝superscript100.2𝐻R_{\rm N}=\bigg{(}\frac{1329\text{ km}}{2\sqrt{p}}\bigg{)}10^{-0.2H}\,.italic_R start_POSTSUBSCRIPT roman_N end_POSTSUBSCRIPT = ( divide start_ARG 1329 km end_ARG start_ARG 2 square-root start_ARG italic_p end_ARG end_ARG ) 10 start_POSTSUPERSCRIPT - 0.2 italic_H end_POSTSUPERSCRIPT . (1)

Eq. (1) is given by Pravec & Harris (2007). The expected radiative nongravitational accelerations are given by

Ai,est=100.2H+bi au/yr2.subscript𝐴𝑖estsuperscript100.2𝐻subscript𝑏𝑖 au/yr2A_{i,{\rm est}}=10^{0.2H+b_{i}}\text{ au/yr${}^{2}$}\,.italic_A start_POSTSUBSCRIPT italic_i , roman_est end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT 0.2 italic_H + italic_b start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT end_POSTSUPERSCRIPT au/yr . (2)

In Eq. (2), bisubscript𝑏𝑖b_{i}italic_b start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT is a constant scaling factor, which is 10.734±0.075plus-or-minus10.7340.075-10.734\pm 0.075- 10.734 ± 0.075 for radiation pressure (i=1𝑖1i=1italic_i = 1) and 12.315±0.014plus-or-minus12.3150.014-12.315\pm 0.014- 12.315 ± 0.014 for the Yarkovsky effect (i=2𝑖2i=2italic_i = 2; see Seligman et al., 2023).444These uncertainties are small compared to the parameter values and are a minimal source of error. Note that this equation depends only on the absolute magnitude H𝐻Hitalic_H and the empirically-determined bisubscript𝑏𝑖b_{i}italic_b start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT and has no dependence on the albedo. Therefore, this estimate of the occurrence rate holds for any subpopulation of asteroids with a representative albedo. These values are empirically obtained as the best-fit line to object data. Although there is no known out-of-plane radiative effect, we assume that the expected out-of-plane nongravitational acceleration (i=3𝑖3i=3italic_i = 3) is equivalent to the Yarkovsky effect as a proxy.

Table 1: Dark Comet Nongravitational Accelerations. The nongravitational accelerations’ magnitudes, significance, and the maximum ratio of measured nongravitational acceleration to expected nongravitational acceleration for all of the dark comets. Nongravitational accelerations and uncertainties are given by Seligman et al. (2023) and the references therein.
Object A1subscript𝐴1A_{1}italic_A start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPTa Signif.b A2subscript𝐴2A_{2}italic_A start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPTc Signif.d A3subscript𝐴3A_{3}italic_A start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPTe Signif.f max(Cisubscript𝐶𝑖C_{i}italic_C start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT)g
[105superscript10510^{-5}10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT au yr-2] [σ𝜎\sigmaitalic_σ] [105superscript10510^{-5}10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT au yr-2] [σ𝜎\sigmaitalic_σ] [105superscript10510^{-5}10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT au yr-2] [σ𝜎\sigmaitalic_σ]
1998 KY26 2.31±plus-or-minus\pm± 1.21 2 -0.00168±plus-or-minus\pm± 0.00081 2 0.426±plus-or-minus\pm± 0.153 3 66.86
2005 VL1 -8.87±plus-or-minus\pm± 10.7 <<<1 -0.00947±plus-or-minus\pm± 0.00789 1 -0.320±plus-or-minus\pm± 0.0546 6 33.90
2016 NJ33 12.4±plus-or-minus\pm± 3.94 3 -0.00754±plus-or-minus\pm± 0.00257 3 1.13±plus-or-minus\pm± 0.217 5 186.38
2010 VL65 8.75±plus-or-minus\pm± 17.3 <<<1 -0.00195±plus-or-minus\pm± 0.00711 <<<1 -1.22±plus-or-minus\pm± 0.173 7 36.01
2010 RF12 0.650±plus-or-minus\pm± 0.795 <<<1 -0.00181±plus-or-minus\pm± 0.00381 <<<1 -0.224±plus-or-minus\pm± 0.028 8 9.58
2006 RH120 1.84±plus-or-minus\pm± 0.107 18 -0.676±plus-or-minus\pm± 0.0849 8 -0.173±plus-or-minus\pm± 0.0426 4 17.58
2003 RM -1.39±plus-or-minus\pm± 1.62 <<<1 0.0286±plus-or-minus\pm± 0.0005 56 0.0200±plus-or-minus\pm± 0.0723 <<<1 679.90
  • a

    Radial nongravitational acceleration.

  • b

    Significance of A1subscript𝐴1A_{1}italic_A start_POSTSUBSCRIPT 1 end_POSTSUBSCRIPT.

  • c

    Transverse nongravitational acceleration.

  • d

    Significance of A2subscript𝐴2A_{2}italic_A start_POSTSUBSCRIPT 2 end_POSTSUBSCRIPT.

  • e

    Out-of-plane nongravitational acceleration.

  • f

    Significance of A3subscript𝐴3A_{3}italic_A start_POSTSUBSCRIPT 3 end_POSTSUBSCRIPT.

  • g

    Maximum value of the ratio of the measured to expected nongravitational acceleration across all three components.

Refer to caption
Figure 1: Cumulative fraction/number of objects with nongravitational accelerations smaller than a given cutoff, which is a factor of C𝐶Citalic_C larger than the estimated nongravitational acceleration from radiation pressure (i=1𝑖1i=1italic_i = 1) or the Yarkovsky effect (i=2𝑖2i=2italic_i = 2, and i=3𝑖3i=3italic_i = 3 by proxy). All three components are shown individually. The number of objects that satisfy Eq. (3) for all 3 components are also plotted. The known dark comets have accelerations similar-to\sim10 – 680×\times× the expected radiative values, which are shown as a shaded region bounded by vertical dashed lines.

We define a cutoff value CAi,est𝐶subscript𝐴𝑖estCA_{i,{\rm est}}italic_C italic_A start_POSTSUBSCRIPT italic_i , roman_est end_POSTSUBSCRIPT, signifying the value where nongravitational accelerations are considered to be significant in comparison to typical values from radiation-driven effects. For example, an object with a radius of R100similar-to-or-equals𝑅100R\simeq 100italic_R ≃ 100 m and an albedo of p=0.1𝑝0.1p=0.1italic_p = 0.1 is expected to have a radiation-pressure acceleration of A1,est4×107similar-tosubscript𝐴1est4superscript107A_{1,{\rm est}}\sim 4\times 10^{-7}italic_A start_POSTSUBSCRIPT 1 , roman_est end_POSTSUBSCRIPT ∼ 4 × 10 start_POSTSUPERSCRIPT - 7 end_POSTSUPERSCRIPT au yr-2. For a given C𝐶Citalic_C, we count the number of objects which have nongravitational acceleration magnitudes smaller than this cutoff by >3σabsent3𝜎>3\sigma> 3 italic_σ. Therefore, the estimated fraction of objects with nongravitational accelerations definitively lower than CAi,est𝐶subscript𝐴𝑖estCA_{i,{\rm est}}italic_C italic_A start_POSTSUBSCRIPT italic_i , roman_est end_POSTSUBSCRIPT is given by the fraction of objects that satisfy

|Ai|+3σiCAi,est.subscript𝐴𝑖3subscript𝜎𝑖𝐶subscript𝐴𝑖est|A_{i}|+3\sigma_{i}\leq CA_{i,{\rm est}}\,.| italic_A start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | + 3 italic_σ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ≤ italic_C italic_A start_POSTSUBSCRIPT italic_i , roman_est end_POSTSUBSCRIPT . (3)

In Eq. (3), σisubscript𝜎𝑖\sigma_{i}italic_σ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT is the uncertainty in the i𝑖iitalic_ith nongravitational acceleration component. In Fig. 1, we plot the number/fraction555The total number is the same for all three components of the nongravitational acceleration by construction of our sample. of objects that satisfy this condition for a range of values of C𝐶Citalic_C. The dark comet accelerations are similar-to\sim10 – 680×\times× the expected radiative values (Table 1), providing an upper and lower bound of the estimated fraction of objects that can be ruled out as dark comets. These boundaries depend on the typical magnitude of a dark comet nongravitational acceleration, and will be further constrained with future discoveries. We find that for C=10𝐶10C=10italic_C = 10680680680680, a range of 5 – 1561 objects satisfy Eq. (3) for all 3 acceleration components simultaneously.666Note that since dark comets are identified by a detection in at least one of the acceleration components, ruling out an object requires a nondetection in all three acceleration components. Therefore, this population represents the subset that can be ruled out as dark comets via this method. Since 7 dark comets are currently known, we estimate that 0.5% – 60% of all sufficiently-constrained NEOs are dark comets. It is important to note that this estimate is subject to small-number statistics and heretofore uncharacterized observational biases. As a result, it is not currently possible to accurately estimate the occurrence rate of dark comets from given data. Future data are necessary to refine this estimate.

If this occurrence rate is representative of the entire NEO population, it is possible that a significant fraction of all NEOs, or tens to thousands of objects, could be dark or dead comets. However, our current lack of understanding of the history, origins, and observational constraints on the dark comets makes an exact calculation impossible. The small sizes of the dark comets also restrict observations to objects with close approaches to Earth and long data arcs, which may not reflect the overall NEO population. As data arcs extend and new dark comets are discovered, the known population will increase and population estimates will become more robust.

3 Material Strengths

In this section, we calculate the minimum tensile strength required to stabilize dark comets against rotational disruption (a basic property of geologic materials). We construct a simple three-force model comparing the centripetal force with two binding forces (material strength and self-gravity) required to hold the rotating object together. We approximate the nucleus as a cube of side length 2s2𝑠2s2 italic_s, with the rotation axis passing through the center of two opposing faces. With this approximation, the cube can be divided into two rectangular prisms joined along a plane passing through the rotation axis, each with dimensions of s𝑠sitalic_s:2s2𝑠2s2 italic_s:2s2𝑠2s2 italic_s (Fig. 2).

While previous work has calculated the internal stress tensors and critical rotation rates for triaxial ellipsoids (Davidsson, 1999, 2001; Breiter et al., 2012), this model allows for an intuitive understanding of the relevant forces, while being sufficient to obtain the tensile strength with errors of order unity. In addition, observational uncertainties dominate the error that arises from our model choice, so improving this model would not necessarily produce a more reliable strength measurement. Our cubic nucleus approximation also has a smaller aspect ratio than many comet nuclei (Knight et al., 2023). The resulting tensile strengths that we calculate are therefore biased downwards, so our calculations provide a minimum strength value. However, we can find a similar constraint by modeling the nucleus as a 2:1 rectangular prism comprised of two cubes, which produces nearly-identical estimates of the material strength. The aspect ratio is therefore not a significant factor in the material strength of these objects. Note that only 1998 KY26 has a measured shape model, which is roughly spherical (Ostro et al., 1999). The rapid rotation rates of 2016 NJ33 and 2006 RH120 imply that these objects likely do not have extreme aspect ratios that make them more susceptible to rotational disruption.

Tensile forces on a homogeneous rotating object of uniform cross-section are maximized along a plane passing through the rotation axis. In the case of fragmentation, we assume the object materially fails and produces the two rectangular prisms as child nuclei.777Although this binary division is simplistic, it simplifies the force calculations while introducing only minor errors. The strength of the gravitational binding force is the force resulting from two identical rectangular prisms of side length s𝑠sitalic_s, 2s2𝑠2s2 italic_s, and 2s2𝑠2s2 italic_s and density ρ𝜌\rhoitalic_ρ in contact with one another such that their centers are separated by a distance s𝑠sitalic_s. The physical binding force Fssubscript𝐹sF_{\rm s}italic_F start_POSTSUBSCRIPT roman_s end_POSTSUBSCRIPT of the nucleus is the product of the contact area of the two cubes and their tensile strength σtsubscript𝜎𝑡\sigma_{t}italic_σ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT. The gravitational force is Gm2/s2𝐺superscript𝑚2superscript𝑠2Gm^{2}/s^{2}italic_G italic_m start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT / italic_s start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT (G𝐺Gitalic_G is the gravitational constant), the tensile strength force is σt(2s)2subscript𝜎𝑡superscript2𝑠2\sigma_{t}(2s)^{2}italic_σ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ( 2 italic_s ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT, and the mass of each prism is m=4ρs3𝑚4𝜌superscript𝑠3m=4\rho s^{3}italic_m = 4 italic_ρ italic_s start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT. The binding force is therefore given by

Fbind=16Gρ2s4+4σts2.subscript𝐹bind16𝐺superscript𝜌2superscript𝑠44subscript𝜎𝑡superscript𝑠2F_{\rm bind}=16G\rho^{2}s^{4}+4\sigma_{t}s^{2}\,.italic_F start_POSTSUBSCRIPT roman_bind end_POSTSUBSCRIPT = 16 italic_G italic_ρ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_s start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT + 4 italic_σ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT italic_s start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT . (4)

For a given rotation rate, the centripetal force needed to keep a mass of m𝑚mitalic_m in a circle at a distance of s/2𝑠2s/2italic_s / 2 from the axis (see Fig. 2) is given by

Fr=2ρs4ω2.subscript𝐹r2𝜌superscript𝑠4superscript𝜔2F_{\rm r}=2\rho s^{4}\omega^{2}.italic_F start_POSTSUBSCRIPT roman_r end_POSTSUBSCRIPT = 2 italic_ρ italic_s start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT . (5)

In Eq. (5), ω𝜔\omegaitalic_ω is the angular velocity, which is related to the rotation period P𝑃Pitalic_P by P=2π/ω𝑃2𝜋𝜔P=2\pi/\omegaitalic_P = 2 italic_π / italic_ω. These objects are minimally stable against rotational disruption when this centripetal force balances the binding forces. The force required at this minimum strength is therefore a lower limit on the bulk tensile strength of these dark comets. We set Eqs. (4) and (5) equal and rearrange to find that the minimum strength is given by

σt=ρs2(12ω24Gρ).subscript𝜎𝑡𝜌superscript𝑠212superscript𝜔24𝐺𝜌\sigma_{t}=\rho s^{2}\left(\frac{1}{2}\omega^{2}-4G\rho\right)\,.italic_σ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT = italic_ρ italic_s start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT ( divide start_ARG 1 end_ARG start_ARG 2 end_ARG italic_ω start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT - 4 italic_G italic_ρ ) . (6)

As most of the dark comets have been observed for a decade or more, we assume they are currently stable against rotational disruption.

Refer to caption
Figure 2: A diagram of our model of a cometary nucleus. We represent a simplified nucleus as a cube with a side length of 2s2𝑠2s2 italic_s. The cube is composed of two identical rectangular prisms with sides of length s𝑠sitalic_s, 2s2𝑠2s2 italic_s, 2s2𝑠2s2 italic_s. The rotation axis passes through the center of this structure. The relevant destabilizing and binding forces are shown in orange and blue, respectively.
Table 2: Dark Comets. Orbital parameters, absolute magnitudes, sizes, rotation periods, and estimated limits on material strengths of the dark comets. Parameters are given by the JPL Horizons database (https://ssd.jpl.nasa.gov/horizons.cgi). Sizes are given by Seligman et al. (2023) and the references therein. The rotational period of 1998 KY26 was reported by Ostro et al. (1999), the period of 2016 NJ33 was reported by Seligman et al. (2023), and the period of 2006 RH120 was reported by Kwiatkowski et al. (2009). The remaining objects do not have reported measurements of rotational period.
Object a𝑎aitalic_aa e𝑒eitalic_eb i𝑖iitalic_ic Hd R𝑅Ritalic_Re P𝑃Pitalic_Pf σt,c>subscript𝜎𝑡𝑐absent\sigma_{t,c}>italic_σ start_POSTSUBSCRIPT italic_t , italic_c end_POSTSUBSCRIPT >g σt,a>subscript𝜎𝑡𝑎absent\sigma_{t,a}>italic_σ start_POSTSUBSCRIPT italic_t , italic_a end_POSTSUBSCRIPT >h
[au] [] [mag] [m] [h] [Pa] [Pa]
1998 KY26 1.23 0.20 1.48 25.60 15 0.178 5.3 26
2005 VL1 0.89 0.23 0.25 26.45 11
2016 NJ33 1.31 0.21 6.64 25.49 16 0.410.410.410.41-1.991.991.991.99 0.00.00.00.0-1.11.11.11.1 0.00.00.00.0-5.65.65.65.6
2010 VL65 1.07 0.14 4.41 29.22 3
2010 RF12 1.06 0.19 0.88 28.42 4
2006 RH120 1.00 0.04 0.31 29.50 2-7 0.046 1.4-18 7.5-91
2003 RM 2.92 0.60 10.86 19.70 230
  • a

    Semimajor axis.

  • b

    Eccentricity.

  • c

    Inclination.

  • d

    Absolute magnitude.

  • e

    Nuclear radius.

  • f

    Rotation period.

  • g

    Minimum material strength for a cometary density of ρ500similar-to-or-equals𝜌500\rho\simeq 500italic_ρ ≃ 500 kg/m3.

  • h

    Minimum material strength for an asteroidal density of ρ2600similar-to-or-equals𝜌2600\rho\simeq 2600italic_ρ ≃ 2600 kg/m3.

Using Eq. (6), we compute this lower limit for the dark comet tensile strengths from their measured rotational periods (currently measured for three dark comets). We assume that the length scale s𝑠sitalic_s is equal to the estimated radius for each. For each dark comet with a known rotational period (1998 KY26, 2016 NJ33, and 2006 RH120), we provide a range of minimum strength values, based on uncertainty in the periods and sizes. We calculate these tensile strengths assuming both (i) a typical comet nucleus density (500 kg/m3, Richardson et al., 2007; Pätzold et al., 2016) and (ii) a typical chondrite meteorite density (2600 kg/m3, Yeomans et al., 2000). These results are presented in Table 2 as σt,c>subscript𝜎𝑡𝑐absent\sigma_{t,c}>italic_σ start_POSTSUBSCRIPT italic_t , italic_c end_POSTSUBSCRIPT > and σt,a>subscript𝜎𝑡𝑎absent\sigma_{t,a}>italic_σ start_POSTSUBSCRIPT italic_t , italic_a end_POSTSUBSCRIPT > respectively.

These computed minimum strengths are on the order of σt,c0.01greater-than-or-equivalent-tosubscript𝜎𝑡𝑐0.01\sigma_{t,c}\gtrsim 0.01italic_σ start_POSTSUBSCRIPT italic_t , italic_c end_POSTSUBSCRIPT ≳ 0.0118181818 Pa or σt,a0.0greater-than-or-equivalent-tosubscript𝜎𝑡𝑎0.0\sigma_{t,a}\gtrsim 0.0italic_σ start_POSTSUBSCRIPT italic_t , italic_a end_POSTSUBSCRIPT ≳ 0.091919191 Pa, assuming the objects have a typical comet density or typical asteroid density, respectively. These strengths are consistent with those of comet nuclei, which are on the order of σt1similar-tosubscript𝜎𝑡1\sigma_{t}\sim 1italic_σ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ∼ 110101010 Pa (Sekanina & Yeomans, 1985; Asphaug & Benz, 1996; Steckloff et al., 2015; Attree et al., 2018), with other estimates reaching as high as σt10200similar-tosubscript𝜎𝑡10200\sigma_{t}\sim 10-200italic_σ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ∼ 10 - 200 Pa for 67P (Hirabayashi et al., 2016). Even for rubble pile bodies held together purely by van der Waals interactions between grains, tensile strengths are on the order of a few tens of Pa (Sánchez & Scheeres, 2014; Rozitis et al., 2014), such as for the DART target’s primary Didymos (Zhang et al., 2017, 2021). Therefore, these strengths are entirely consistent with rubble-pile bodies, and no additional cohesion is required.

However, dark comets are much smaller than typical cometary nuclei. Because the binding force terms have different dependencies on the size s𝑠sitalic_s of the object, dark comets may exhibit fundamentally different evolutionary behavior if they are predominantly bound by a different force than the larger cometary nuclei. This is evidenced in the “spin barrier” for larger comets and asteroids, in which rubble piles (gravity-dominated objects) break apart at rotation periods shorter than 2.2similar-toabsent2.2\sim 2.2∼ 2.2 hours for asteroids (Pravec et al., 2002; Warner et al., 2009; Sánchez & Scheeres, 2014) or 6similar-toabsent6\sim 6∼ 6 hours for comets (Kokotanekova et al., 2017; Safrit et al., 2021). However, smaller, strength-dominated objects can rotate more rapidly if they have significant material strength.

We estimate the approximate size at which objects transition between strength- and gravity-dominated regimes. To do this we calculate the angular velocity at which the rotational forces overcome the tensile strength and gravitational forces and set these velocities equal to each other.

The rotational failure criterion for a strength-dominated object is given by

ωcritstr=2σtρs2.superscriptsubscript𝜔critstr2subscript𝜎𝑡𝜌superscript𝑠2\omega_{\rm crit}^{\rm str}=\sqrt{\frac{2\sigma_{t}}{\rho s^{2}}}\,.italic_ω start_POSTSUBSCRIPT roman_crit end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_str end_POSTSUPERSCRIPT = square-root start_ARG divide start_ARG 2 italic_σ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT end_ARG start_ARG italic_ρ italic_s start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG end_ARG . (7)

For a derivation of this relationship, see Steckloff & Jacobson (2016). Similarly, the failure criterion for a gravity-dominated object is

ωcritgrav=4πρG3.superscriptsubscript𝜔critgrav4𝜋𝜌𝐺3\omega_{\rm crit}^{\rm grav}=\sqrt{\frac{4\pi\rho G}{3}}\,.italic_ω start_POSTSUBSCRIPT roman_crit end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_grav end_POSTSUPERSCRIPT = square-root start_ARG divide start_ARG 4 italic_π italic_ρ italic_G end_ARG start_ARG 3 end_ARG end_ARG . (8)

Notably, this is independent of object radius (Pravec & Harris, 2000; Safrit et al., 2021), although it will depend on the object’s aspect ratio. Setting these two equations equal to one another provides an expression for the strength-to-gravity transition size — R200similar-to𝑅200R\sim 200italic_R ∼ 200 m for asteroids and R1similar-to𝑅1R\sim 1italic_R ∼ 1 km for comets (assuming tensile strength of 10 Pa and densities of 2600 kg/m3 and 500 kg/m3 respectively). This transition size can be scaled to different densities and tensile strengths, such that

stransast=206 m (σt10 Pa)1/2(ρ2600 kg m3)1.superscriptsubscript𝑠transast206 m superscriptsubscript𝜎𝑡10 Pa12superscript𝜌2600superscript kg m31s_{\rm trans}^{\rm ast}=206\text{ m }\left(\frac{\sigma_{t}}{10\text{ Pa}}% \right)^{1/2}\left(\frac{\rho}{2600\text{ kg m}^{-3}}\right)^{-1}\,.italic_s start_POSTSUBSCRIPT roman_trans end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_ast end_POSTSUPERSCRIPT = 206 m ( divide start_ARG italic_σ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT end_ARG start_ARG 10 Pa end_ARG ) start_POSTSUPERSCRIPT 1 / 2 end_POSTSUPERSCRIPT ( divide start_ARG italic_ρ end_ARG start_ARG 2600 kg m start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT end_ARG ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT . (9)

Similarly, the transition size for comets is

stranscomet=1070 m (σt10 Pa)1/2(ρ500 kg m3)1.superscriptsubscript𝑠transcomet1070 m superscriptsubscript𝜎𝑡10 Pa12superscript𝜌500superscript kg m31s_{\rm trans}^{\rm comet}=1070\text{ m }\left(\frac{\sigma_{t}}{10\text{ Pa}}% \right)^{1/2}\left(\frac{\rho}{500\text{ kg m}^{-3}}\right)^{-1}\,.italic_s start_POSTSUBSCRIPT roman_trans end_POSTSUBSCRIPT start_POSTSUPERSCRIPT roman_comet end_POSTSUPERSCRIPT = 1070 m ( divide start_ARG italic_σ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT end_ARG start_ARG 10 Pa end_ARG ) start_POSTSUPERSCRIPT 1 / 2 end_POSTSUPERSCRIPT ( divide start_ARG italic_ρ end_ARG start_ARG 500 kg m start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT end_ARG ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT . (10)

This asteroid transition radius is consistent with that calculated by Pravec et al. (2002). This suggests that all the dark comets are small enough to lie in the strength-dominated regime (Fig. 3), in contrast to the comet nuclei visited by spacecraft, all of which are large enough to be gravity-dominated: 9P/Tempel 1 (A’Hearn et al., 2005), 19P/Borrelly (Soderblom et al., 2002), 67P/Churyumov-Gerasimenko (Jorda et al., 2016), and 81P/Wild 2 (Brownlee et al., 2006), with 103P/Hartley 2 (s600similar-to𝑠600s\sim 600italic_s ∼ 600 m; A’Hearn et al., 2011) possibly representing a transition object.

Refer to caption
Figure 3: Minimum stable rotation period (fastest stable rotation rate) as a function of size and tensile strength. Here we use our simple cubic model to compute the binding forces (self-gravity and strength) to determine the size-dependence of the minimum stable rotation period for various tensile strengths. Note that a larger aspect ratio would increase the tensile strength required to maintain the cohesion of the body. We find that the strength of a typical comet nucleus (similar-to\sim1–100 Pa; Sekanina & Yeomans, 1985; Asphaug & Benz, 1996; Bowling et al., 2014; Steckloff et al., 2015; Hirabayashi et al., 2016; Attree et al., 2018) can hold small comet nuclei together at rotation rates faster than those of the dark comets. Therefore, dark comets do not require atypical strengths, in spite of their short rotation periods/fast rotation rates. The three dark comets with measured rotation periods are shown as white diamonds. JFCs from Kokotanekova et al. (2017) are shown as dark gray squares, and MBCs from Jewitt & Hsieh (2022) are shown as light gray circles. As the radius gets larger, gravity begins to dominate over strength, so the minimum stable period converges.
Refer to caption
Figure 4: The proposed evolutionary track for the dark comets — a comet is injected into the near-Earth environment, where outgassing causes spinup. The spinup causes a rotational disintegration cascade and devolatilization, leading to smaller, less active objects (i.e., dark comets). Extinction is also possible at this point in the evolutionary track, well before becoming an NEO (4b). If the disintegration cascade operates more rapidly than devolatilization (Sec. 4.2), then the object will eventually be destroyed completely. However, if volatiles are removed before this can take place, then the end state will be a small, rapidly-rotating, inactive asteroid in the near-Earth environment. Example objects are also given, which represent analogous physical (not dynamical) stages of this model.

4 Evolution of the Dark Comets

In this section, we discuss the evolutionary pathway of the dark comets suggested by the limits on material strengths calculated in Sec. 3.

4.1 Production via Rotational Fragmentation Cascade

A rotational fragmentation cascade may explain the small sizes and rapid rotations of the dark comets. The weak nongravitational accelerations of these objects suggest that they are highly thermally evolved and largely devolatilized or covered with refractory mantles, if they are cometary in origin. As such objects devolatilize, sublimative torques operating on their nuclei would produce spinup, potentially to their rotational disruption limits (Jewitt, 1997; Steckloff & Jacobson, 2016; Hirabayashi et al., 2016; Jewitt, 2021; Safrit et al., 2021; Jewitt, 2022). Such rotational disruption could spall off pieces of material from the nucleus. So long as these ejected pieces have less than similar-to\sim20% of the mass of the original nucleus, the ejected piece will gravitationally escape the parent nucleus (Hirabayashi et al., 2016), and may continue to rotationally disrupt.

This rotational cascade could produce small fragments that reside in or near the strength-dominated regime. In this case, breakup would produce even smaller fragments that are more resistant to further disruption, because smaller objects require even faster rotation to disrupt (see Eq. 6). However, continual outgassing torque would continue to spin fragments up to their disruption limits. Eventually, these objects would be sufficiently small and devolatilized to remain stable against further disruption. As a result, one would expect rotationally-disrupted objects to continually spin up until all that remains are small fragments that are stable to further disruption. This final state is consistent with the observed properties of dark comets, which are stable against rotational disruption if they have typical cometary bulk tensile strengths. Therefore, dark comets may represent one of the final endpoints of comet nucleus evolution — the small (s110similar-to𝑠110s\sim 1-10italic_s ∼ 1 - 10 m), largely devolatilized remnants of their parent nuclei.

This implies a possible evolutionary sequence for dark comets — they begin as part of a volatile-rich body that fragments due to rotational disruption. The resulting fragments may then undergo a rotational spin-up cascade until they reach small sizes that are stable against disruption, where they remain until they become devolatilized. At some point during this evolutionary sequence, these objects would dynamically transfer from their original environment to an orbit with a semimajor axis a1similar-to-or-equals𝑎1a\simeq 1italic_a ≃ 1 au. We show this hypothetical evolutionary sequence in Fig. 4.

Comet 252P/LINEAR may be undergoing the first steps toward evolving into a dark comet. 252P is thought to be a fragment of comet 460P/2016 BA14 Pan-STARRS, due to their similar orbits (Li et al., 2017). Their relative sizes (300 m for 252P and 1 km for 460P; Li et al., 2017) are consistent with this origin, according to the conditions found by Hirabayashi et al. (2016). 252P is also a fast rotator for a comet nucleus, with a light curve–derived rotation period of only 5.41 h (Li et al., 2017). However, it is worth noting that 460P may have a rotation period as slow as 40 h (Naidu et al., 2016). Because the size ratio of 252P and 460P is 1/3similar-toabsent13\sim 1/3∼ 1 / 3, 252P may carry away enough angular momentum to slow the rotation of 460P to this rate. Therefore, this slow rotation rate is compatible with this fragmentation mechanism.

As the sizes of the fragments decrease, they may continue to spin up. Some of these fragments may undergo a sublimation-driven rotational cascade that disrupts them into fine debris, potentially producing dust trails that resemble isolated striae (Steckloff & Jacobson, 2016). Alternatively, these fragments may mostly devolatilize prior to the completion of the rotational fragmentation cascade. In such cases, the fragments would reach small sizes with fast rotation periods and minimal outgassing — the dark comets.888Note that the small number of measured rotation rates on the dark comets makes any comparison highly uncertain.

This suggests that the near-Earth space may be littered with numerous dead (fully devolatilized) comet fragments. Nevertheless, they may be exceedingly difficult to detect, or identify if already detected. In addition to being small and correspondingly dim, the weak sublimation-driven nongravitational accelerations that are the currently-identified hallmark of dark comets (Farnocchia et al., 2023; Seligman et al., 2023) will fade until their nuclei are completely devolatilized. Thereafter, dark comets would appear only as small, inactive asteroids/meteoroids (“dead comets”). The detected dark comets may therefore represent objects finishing their transition into this ultimate end state.

This is a significantly different evolutionary path than that taken by larger comets that are predominantly held together by self-gravity (comet nuclei 500greater-than-or-equivalent-toabsent500\gtrsim 500≳ 500 m in radius, based on Safrit et al., 2021).999Bottke et al. (2023) found that the strength-gravity transition was at R20similar-to𝑅20R\sim 20italic_R ∼ 20 m, using a model of collisional evolution of Jupiter Trojans and outer solar system populations, combined with laboratory data. This value is consistent with our results for σt0.1similar-to-or-equalssubscript𝜎𝑡0.1\sigma_{t}\simeq 0.1italic_σ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ≃ 0.1 Pa (see Fig. 3). Such large comet nuclei predominantly evolve via mass-wasting events, which flatten the surface of the nucleus over time, expose buried volatiles, and maintain sublimative activity (Vincent et al., 2015; Steckloff & Samarasinha, 2018). As these nuclei flatten, the frequency of outbursts decreases, leading to episodic activity that can reactivate comet nuclei (Steckloff & Samarasinha, 2018). Eventually, these nuclei can become sufficiently gravitationally flat to preclude further mass-wasting reactivation, and lapse into large dead comets (Steckloff & Samarasinha, 2018). In addition, the weak relative tensile strength of comets at gravity-dominated scales causes comets to be destroyed in the near-Earth environment, preventing gravity-dominated objects from remaining as NEOs for long periods of time, since close encounters will tend to destroy them (Levison & Duncan, 1997; Fernández et al., 2002; Di Sisto et al., 2009; Granvik & Walsh, 2024; Nesvorný et al., 2024a).

However, throughout this process, gravity-dominated nuclei generally do not change size. Sublimative spin-up would tend to deform the nucleus (Safrit et al., 2021) or fission off gravitationally bound fragments that can reaccrete (Hirabayashi et al., 2016). This is fundamentally different from the small strength-bound dark comets, which can more easily spin up and completely divide. Although both of these tracks end with a largely devolatilized nucleus, their sizes and final spin rates are likely to diverge based on the track followed. Under this hypothesis, the evolutionary path of a comet is strongly dependent on the size of its nucleus (and its dynamical history), which is consistent with the JFCs, MBCs, and dark comets shown in Fig. 3.

4.2 Devolatilization and Fission Timescales

The timescale over which these objects would devolatilize τdevsubscript𝜏dev\tau_{\rm dev}italic_τ start_POSTSUBSCRIPT roman_dev end_POSTSUBSCRIPT is set, at minimum, by the time required for a thermal heat pulse to reach the centers of these objects. Safrit et al. (2021) found that

τdev=R2α.subscript𝜏devsuperscript𝑅2𝛼\tau_{\rm dev}=\frac{R^{2}}{\alpha}\,.italic_τ start_POSTSUBSCRIPT roman_dev end_POSTSUBSCRIPT = divide start_ARG italic_R start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT end_ARG start_ARG italic_α end_ARG . (11)

In Eq. (11), R𝑅Ritalic_R is the radius of the dark comet nucleus and α𝛼\alphaitalic_α is the thermal diffusivity of the material, typically on the order of α108107similar-to𝛼superscript108superscript107\alpha\sim 10^{-8}-10^{-7}italic_α ∼ 10 start_POSTSUPERSCRIPT - 8 end_POSTSUPERSCRIPT - 10 start_POSTSUPERSCRIPT - 7 end_POSTSUPERSCRIPT m2 s-1 (Gundlach & Blum, 2012; Jewitt et al., 2017; Groussin et al., 2019; Steckloff et al., 2021). Thus, for dark comets with radii between 2 and 15 m, τdev11000similar-tosubscript𝜏dev11000\tau_{\rm dev}\sim 1-1000italic_τ start_POSTSUBSCRIPT roman_dev end_POSTSUBSCRIPT ∼ 1 - 1000 years.

We can compare this to the timescale required for sublimative torques to spin up an object, τtorsubscript𝜏tor\tau_{\rm tor}italic_τ start_POSTSUBSCRIPT roman_tor end_POSTSUBSCRIPT. Using the sublimative YORP (SYORP) formalism from Steckloff & Jacobson (2016), we write that

τtor=R32ρσtPSYS.subscript𝜏tor𝑅32𝜌subscript𝜎𝑡subscript𝑃𝑆subscript𝑌𝑆\tau_{\rm tor}=\frac{R\sqrt{32\rho\sigma_{t}}}{P_{S}Y_{S}}\,.italic_τ start_POSTSUBSCRIPT roman_tor end_POSTSUBSCRIPT = divide start_ARG italic_R square-root start_ARG 32 italic_ρ italic_σ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT end_ARG end_ARG start_ARG italic_P start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT italic_Y start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT end_ARG . (12)

In Eq. (12), PSsubscript𝑃𝑆P_{S}italic_P start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT is the material sublimation pressure (Steckloff et al., 2015; Steckloff & Jacobson, 2016) and YSsubscript𝑌𝑆Y_{S}italic_Y start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT is the SYORP parameter, which measures the fraction of the sublimative momentum flux that must be directed tangentially to produce an equivalent torque. In the near-Earth environment, the sublimation pressure PSsubscript𝑃𝑆P_{S}italic_P start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT for water ice is on the order of PS0.01similar-tosubscript𝑃𝑆0.01P_{S}\sim 0.01italic_P start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT ∼ 0.010.10.10.10.1 Pa (Steckloff & Jacobson, 2016). SYORP parameters for objects of this size are thought to be on the order of YS104similar-tosubscript𝑌𝑆superscript104Y_{S}\sim 10^{-4}italic_Y start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT ∼ 10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT103superscript10310^{-3}10 start_POSTSUPERSCRIPT - 3 end_POSTSUPERSCRIPT (Steckloff & Jacobson, 2016), but have been measured in larger, strength-dominated comet nuclei to be an order of magnitude lower, at YS105similar-tosubscript𝑌𝑆superscript105Y_{S}\sim 10^{-5}italic_Y start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT ∼ 10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT104superscript10410^{-4}10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT (Steckloff et al., 2021). This results in SYORP timescales (time required to monotonically spin up from rest to disruption) of 2200similar-toabsent2200\sim 2-200∼ 2 - 200 years for the dark comets with known rotation periods. However, this may be up to an order of magnitude longer if SYORP parameters are weaker. Furthermore, the stochastic nature of torque directions on comets over secular timescales (Hirabayashi et al., 2016) may lengthen this by an additional order of magnitude.

Regardless, this spin-up timescale is comparable to the devolatilization timescale. By dividing Eqs. (11) and (12), and plugging in typical values, we find that

τdevτtor=1.25(R5 m)(PS0.1 Pa)(YS104)(α107 m2/s)1×(ρ500 kg/m3)1/2(σt10 Pa)1/2.subscript𝜏devsubscript𝜏tor1.25𝑅5 msubscript𝑃𝑆0.1 Pasubscript𝑌𝑆superscript104superscript𝛼superscript107 m2/s1superscript𝜌500 kg/m312superscriptsubscript𝜎𝑡10 Pa12\begin{split}\frac{\tau_{\rm dev}}{\tau_{\rm tor}}=1.25&\Big{(}\frac{R}{5\text% { m}}\Big{)}\Big{(}\frac{P_{S}}{0.1\text{ Pa}}\Big{)}\Big{(}\frac{Y_{S}}{10^{-% 4}}\Big{)}\Big{(}\frac{\alpha}{10^{-7}\text{ m${}^{2}$/s}}\Big{)}^{-1}\\ &\times\Big{(}\frac{\rho}{500\text{ kg/m${}^{3}$}}\Big{)}^{-1/2}\Big{(}\frac{% \sigma_{t}}{10\text{ Pa}}\Big{)}^{-1/2}\,.\end{split}start_ROW start_CELL divide start_ARG italic_τ start_POSTSUBSCRIPT roman_dev end_POSTSUBSCRIPT end_ARG start_ARG italic_τ start_POSTSUBSCRIPT roman_tor end_POSTSUBSCRIPT end_ARG = 1.25 end_CELL start_CELL ( divide start_ARG italic_R end_ARG start_ARG 5 m end_ARG ) ( divide start_ARG italic_P start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT end_ARG start_ARG 0.1 Pa end_ARG ) ( divide start_ARG italic_Y start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT end_ARG start_ARG 10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT end_ARG ) ( divide start_ARG italic_α end_ARG start_ARG 10 start_POSTSUPERSCRIPT - 7 end_POSTSUPERSCRIPT m /s end_ARG ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT end_CELL end_ROW start_ROW start_CELL end_CELL start_CELL × ( divide start_ARG italic_ρ end_ARG start_ARG 500 kg/m end_ARG ) start_POSTSUPERSCRIPT - 1 / 2 end_POSTSUPERSCRIPT ( divide start_ARG italic_σ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT end_ARG start_ARG 10 Pa end_ARG ) start_POSTSUPERSCRIPT - 1 / 2 end_POSTSUPERSCRIPT . end_CELL end_ROW (13)

It is then plausible that while many of these objects can experience rotational disintegration, many others will devolatilize prior to completing a fragmentation cascade, as suggested by Steckloff & Jacobson (2016). In this case, the cascade would terminate with a small, rapidly-rotating dead comet — a dark comet without any nongravitational acceleration.

Refer to caption
Figure 5: The normalized size frequency distribution (SFD) of the objects evolving under this model, for a range of removal rate parameters q0subscript𝑞0q_{0}italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT. We show q0[1014,106]subscript𝑞0superscript1014superscript106q_{0}\in[10^{-14},10^{-6}]italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ∈ [ 10 start_POSTSUPERSCRIPT - 14 end_POSTSUPERSCRIPT , 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT ] s-1. We expect q0=1014subscript𝑞0superscript1014q_{0}=10^{-14}italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT - 14 end_POSTSUPERSCRIPT s-1 to be the true removal rate, and q0=106subscript𝑞0superscript106q_{0}=10^{-6}italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT - 6 end_POSTSUPERSCRIPT s-1 is chosen because the distribution is flat at this point. The model SFD f(R)𝑓𝑅f(R)italic_f ( italic_R ) is shown in terms of the steady-state divisionless SFD, p(R)/q0𝑝𝑅subscript𝑞0p(R)/q_{0}italic_p ( italic_R ) / italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT. We assume that the injection SFD power-law slope is β=1𝛽1\beta=1italic_β = 1, but this only affects the exact values rather than the slope. The predicted small-R limit of f(R)Rp(R)proportional-to𝑓𝑅𝑅𝑝𝑅f(R)\propto R\,p(R)italic_f ( italic_R ) ∝ italic_R italic_p ( italic_R ) is shown as a red dashed line.

5 Size Frequency Distribution

In this section, we estimate a size frequency distribution (SFD) for objects undergoing the rotational fragmentation cascade in the NEO environment discussed in Sec. 4.

We consider an idealized case where the rotational cascade only splits objects in half. Specifically, each splitting event replaces an object of radius R𝑅Ritalic_R with two objects of radius R/(21/3)𝑅superscript213R/(2^{1/3})italic_R / ( 2 start_POSTSUPERSCRIPT 1 / 3 end_POSTSUPERSCRIPT ). We ignore collisional fragmentation in these calculations because the collisional timescale (Bottke et al., 2005) is

τcol=1 Myr (R10 m)1/2.subscript𝜏col1 Myr superscript𝑅10 m12\tau_{\rm col}=1\text{ Myr }\left(\frac{R}{10\text{ m}}\right)^{1/2}\,.italic_τ start_POSTSUBSCRIPT roman_col end_POSTSUBSCRIPT = 1 Myr ( divide start_ARG italic_R end_ARG start_ARG 10 m end_ARG ) start_POSTSUPERSCRIPT 1 / 2 end_POSTSUPERSCRIPT . (14)

Since τdivτcolmuch-less-thansubscript𝜏divsubscript𝜏col\tau_{\rm div}\ll\tau_{\rm col}italic_τ start_POSTSUBSCRIPT roman_div end_POSTSUBSCRIPT ≪ italic_τ start_POSTSUBSCRIPT roman_col end_POSTSUBSCRIPT for the relevant sizes, collisions are negligible. We then define the differential object fraction per radius R𝑅Ritalic_R as

f(R)dndR.𝑓𝑅d𝑛d𝑅f(R)\equiv\frac{{\rm d}n}{{\rm d}R}\,.italic_f ( italic_R ) ≡ divide start_ARG roman_d italic_n end_ARG start_ARG roman_d italic_R end_ARG . (15)

This is normalized such that the total number of objects N=0f(R)dR𝑁superscriptsubscript0𝑓𝑅differential-d𝑅N=\int_{0}^{\infty}\!\!f(R)\,{\rm d}Ritalic_N = ∫ start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT italic_f ( italic_R ) roman_d italic_R.

Next, we define a function p(R)𝑝𝑅p(R)italic_p ( italic_R ), which is the number of objects of a size R𝑅Ritalic_R that are injected into the population over a time step dtd𝑡{\rm d}troman_d italic_t. This function is set by the input from the source population into the NEO environment. Note that this analysis is agnostic to the source population and is applicable to a JFC or MBC source. We further define a constant q0subscript𝑞0q_{0}italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT, which is the fraction of objects that are removed from the population in a given time dtd𝑡{\rm d}troman_d italic_t by orbital instabilities. This is assumed to be constant across all object sizes.

Finally, the population of objects of a given size is determined by the rotational cascade. We define a function q(R)𝑞𝑅q(R)italic_q ( italic_R ), which is the fraction of objects of size R𝑅Ritalic_R that divide in a time step dtd𝑡{\rm d}troman_d italic_t. Combining these definitions, we find that the change in the population is given by

dfdt(R)=p(R)+2q(ϵR)f(ϵR)q(R)f(R)q0f(R).d𝑓d𝑡𝑅𝑝𝑅2𝑞italic-ϵ𝑅𝑓italic-ϵ𝑅𝑞𝑅𝑓𝑅subscript𝑞0𝑓𝑅\frac{{\rm d}f}{{\rm d}t}(R)=p(R)+2q(\epsilon R)f(\epsilon R)-q(R)f(R)-q_{0}f(% R)\,.divide start_ARG roman_d italic_f end_ARG start_ARG roman_d italic_t end_ARG ( italic_R ) = italic_p ( italic_R ) + 2 italic_q ( italic_ϵ italic_R ) italic_f ( italic_ϵ italic_R ) - italic_q ( italic_R ) italic_f ( italic_R ) - italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT italic_f ( italic_R ) . (16)

The first term on the right-hand side is the number of objects injected into the population, the second term is the addition of two objects via splitting of an object with size ϵRitalic-ϵ𝑅\epsilon Ritalic_ϵ italic_R, the third term is the removal of objects via fragmentation events, and the final term is the removal of objects via orbital instabilities. In this case, ϵ=21/3italic-ϵsuperscript213\epsilon=2^{1/3}italic_ϵ = 2 start_POSTSUPERSCRIPT 1 / 3 end_POSTSUPERSCRIPT. In order to find the steady-state SFD, we set df/dt=0d𝑓d𝑡0{\rm d}f/{\rm d}t=0roman_d italic_f / roman_d italic_t = 0.

We assume that p(R)𝑝𝑅p(R)italic_p ( italic_R ) is constant in time, so that the age distribution of objects is uniform. Therefore, if objects live in the NEO environment for 106superscript10610^{6}10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT yr (Gladman et al., 2000; Nesvorný et al., 2023), then q01014similar-to-or-equalssubscript𝑞0superscript1014q_{0}\simeq 10^{-14}italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ≃ 10 start_POSTSUPERSCRIPT - 14 end_POSTSUPERSCRIPT s-1.

There are two competing factors to determine q(R)𝑞𝑅q(R)italic_q ( italic_R ): spin-up and devolatilization. The timescales for these processes are generally comparable (Sec. 4.2, Eq. (13)), but fragmentation does not occur when the devolatilization timescale is shorter than the spinup timescale. However, if the spinup timescale is shorter than the devolatilization timescale, then f(R)dt/τdiv(R)𝑓𝑅d𝑡subscript𝜏div𝑅f(R){\rm d}t/\tau_{\rm div}(R)italic_f ( italic_R ) roman_d italic_t / italic_τ start_POSTSUBSCRIPT roman_div end_POSTSUBSCRIPT ( italic_R ) objects are removed, where τdivsubscript𝜏div\tau_{\rm div}italic_τ start_POSTSUBSCRIPT roman_div end_POSTSUBSCRIPT is the spinup timescale. We assume that all objects begin at the critical rotation rate of the progenitor object, even if they were recently inserted into the population. Therefore, τdiv=τtor(11/ϵ)subscript𝜏divsubscript𝜏tor11italic-ϵ\tau_{\rm div}=\tau_{\rm tor}(1-1/\epsilon)italic_τ start_POSTSUBSCRIPT roman_div end_POSTSUBSCRIPT = italic_τ start_POSTSUBSCRIPT roman_tor end_POSTSUBSCRIPT ( 1 - 1 / italic_ϵ ), where τtorsubscript𝜏tor\tau_{\rm tor}italic_τ start_POSTSUBSCRIPT roman_tor end_POSTSUBSCRIPT is given by Eq. (12), and the factor of (121/3)1superscript213(1-2^{-1/3})( 1 - 2 start_POSTSUPERSCRIPT - 1 / 3 end_POSTSUPERSCRIPT ) is introduced by beginning at ω=ωcrit(ϵR)𝜔subscript𝜔crititalic-ϵ𝑅\omega=\omega_{\rm crit}(\epsilon R)italic_ω = italic_ω start_POSTSUBSCRIPT roman_crit end_POSTSUBSCRIPT ( italic_ϵ italic_R ) rather than ω=0𝜔0\omega=0italic_ω = 0. With the devolatilization timescale given by Eq. (11), q(R)𝑞𝑅q(R)italic_q ( italic_R ) is given by

q(R)={PSYS32ρσt1R1(11/ϵ)τdivτdev0otherwise.𝑞𝑅casessubscript𝑃𝑆subscript𝑌𝑆32𝜌subscript𝜎𝑡1𝑅111italic-ϵsubscript𝜏divsubscript𝜏dev0otherwiseq(R)=\begin{dcases}\frac{P_{S}Y_{S}}{\sqrt{32\rho\sigma_{t}}}\frac{1}{R}\frac{% 1}{(1-1/\epsilon)}&\tau_{\rm div}\leq\tau_{\rm dev}\\ 0&\text{otherwise}\,.\end{dcases}italic_q ( italic_R ) = { start_ROW start_CELL divide start_ARG italic_P start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT italic_Y start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT end_ARG start_ARG square-root start_ARG 32 italic_ρ italic_σ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT end_ARG end_ARG divide start_ARG 1 end_ARG start_ARG italic_R end_ARG divide start_ARG 1 end_ARG start_ARG ( 1 - 1 / italic_ϵ ) end_ARG end_CELL start_CELL italic_τ start_POSTSUBSCRIPT roman_div end_POSTSUBSCRIPT ≤ italic_τ start_POSTSUBSCRIPT roman_dev end_POSTSUBSCRIPT end_CELL end_ROW start_ROW start_CELL 0 end_CELL start_CELL otherwise . end_CELL end_ROW (17)

Solving Eq. (16) for f(R)𝑓𝑅f(R)italic_f ( italic_R ) is complicated by the fact that both f(R)𝑓𝑅f(R)italic_f ( italic_R ) and f(ϵR)𝑓italic-ϵ𝑅f(\epsilon R)italic_f ( italic_ϵ italic_R ) appear in the equation. Moreover, a single solution does not exist because q(R)𝑞𝑅q(R)italic_q ( italic_R ) is a piecewise function of radius R𝑅Ritalic_R. We therefore consider three domains of the function.

Case 1. First we consider the case where R𝑅Ritalic_R is small and q(R)=q(ϵR)=0𝑞𝑅𝑞italic-ϵ𝑅0q(R)=q(\epsilon R)=0italic_q ( italic_R ) = italic_q ( italic_ϵ italic_R ) = 0. In this instance, Eq. (16) can be easily solved to show that

f(R)=p(R)q0.𝑓𝑅𝑝𝑅subscript𝑞0f(R)=\frac{p(R)}{q_{0}}\,.italic_f ( italic_R ) = divide start_ARG italic_p ( italic_R ) end_ARG start_ARG italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG . (18)

In this limit, f(R)𝑓𝑅f(R)italic_f ( italic_R ) has identical behavior to the input function, albeit modified by a factor that accounts for the constant removal by orbital instability.

Case 2. Next, we consider the case where R𝑅Ritalic_R is large and q(R)0𝑞𝑅0q(R)\neq 0italic_q ( italic_R ) ≠ 0, at the extreme end from the previous case. We can rewrite Eq. (16), solving for f(r)𝑓𝑟f(r)italic_f ( italic_r ) in terms of f(ϵR)𝑓italic-ϵ𝑅f(\epsilon R)italic_f ( italic_ϵ italic_R ) to obtain

f(R)=p(R)+2q(ϵR)f(ϵR)q0+q(R).𝑓𝑅𝑝𝑅2𝑞italic-ϵ𝑅𝑓italic-ϵ𝑅subscript𝑞0𝑞𝑅f(R)=\frac{p(R)+2q(\epsilon R)f(\epsilon R)}{q_{0}+q(R)}\,.italic_f ( italic_R ) = divide start_ARG italic_p ( italic_R ) + 2 italic_q ( italic_ϵ italic_R ) italic_f ( italic_ϵ italic_R ) end_ARG start_ARG italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + italic_q ( italic_R ) end_ARG . (19)

We continue this series, writing f(ϵR)𝑓italic-ϵ𝑅f(\epsilon R)italic_f ( italic_ϵ italic_R ) in terms of f(ϵ2R)𝑓superscriptitalic-ϵ2𝑅f(\epsilon^{2}R)italic_f ( italic_ϵ start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT italic_R ) up to infinity. From Eq. (17), q(R)=A/R𝑞𝑅𝐴𝑅q(R)=A/Ritalic_q ( italic_R ) = italic_A / italic_R where A=(121/3)1PSYS/32ρσt𝐴superscript1superscript2131subscript𝑃𝑆subscript𝑌𝑆32𝜌subscript𝜎𝑡A=(1-2^{-1/3})^{-1}P_{S}Y_{S}/\sqrt{32\rho\sigma_{t}}italic_A = ( 1 - 2 start_POSTSUPERSCRIPT - 1 / 3 end_POSTSUPERSCRIPT ) start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT italic_P start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT italic_Y start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT / square-root start_ARG 32 italic_ρ italic_σ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT end_ARG is a constant. Therefore, q(ϵkR)=ϵkq(R)𝑞superscriptitalic-ϵ𝑘𝑅superscriptitalic-ϵ𝑘𝑞𝑅q(\epsilon^{k}R)=\epsilon^{-k}q(R)italic_q ( italic_ϵ start_POSTSUPERSCRIPT italic_k end_POSTSUPERSCRIPT italic_R ) = italic_ϵ start_POSTSUPERSCRIPT - italic_k end_POSTSUPERSCRIPT italic_q ( italic_R ). We also assume that p(R)𝑝𝑅p(R)italic_p ( italic_R ) is a power law such that p(R)Rβproportional-to𝑝𝑅superscript𝑅𝛽p(R)\propto R^{-\beta}italic_p ( italic_R ) ∝ italic_R start_POSTSUPERSCRIPT - italic_β end_POSTSUPERSCRIPT. Continuing out the series and simplifying, we find that

f(R)=p(R)n=02nϵn(2β+n+1)/2q(R)n×k=0n[q0+q(ϵkR)]1.𝑓𝑅𝑝𝑅superscriptsubscript𝑛0superscript2𝑛superscriptitalic-ϵ𝑛2𝛽𝑛12𝑞superscript𝑅𝑛superscriptsubscriptproduct𝑘0𝑛superscriptdelimited-[]subscript𝑞0𝑞superscriptitalic-ϵ𝑘𝑅1\begin{split}f(R)=p(R)\sum_{n=0}^{\infty}&2^{n}\epsilon^{-n(2\beta+n+1)/2}q(R)% ^{n}\\ \times&\prod_{k=0}^{n}\Big{[}q_{0}+q(\epsilon^{k}R)\Big{]}^{-1}\,.\end{split}start_ROW start_CELL italic_f ( italic_R ) = italic_p ( italic_R ) ∑ start_POSTSUBSCRIPT italic_n = 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT ∞ end_POSTSUPERSCRIPT end_CELL start_CELL 2 start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT italic_ϵ start_POSTSUPERSCRIPT - italic_n ( 2 italic_β + italic_n + 1 ) / 2 end_POSTSUPERSCRIPT italic_q ( italic_R ) start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT end_CELL end_ROW start_ROW start_CELL × end_CELL start_CELL ∏ start_POSTSUBSCRIPT italic_k = 0 end_POSTSUBSCRIPT start_POSTSUPERSCRIPT italic_n end_POSTSUPERSCRIPT [ italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT + italic_q ( italic_ϵ start_POSTSUPERSCRIPT italic_k end_POSTSUPERSCRIPT italic_R ) ] start_POSTSUPERSCRIPT - 1 end_POSTSUPERSCRIPT . end_CELL end_ROW (20)

In general, Eq. (20) can be used to compute f(R)𝑓𝑅f(R)italic_f ( italic_R ) for the domain where q(R)0𝑞𝑅0q(R)\neq 0italic_q ( italic_R ) ≠ 0, so long as the series converges.

Case 3. In the final case, q(R)=0𝑞𝑅0q(R)=0italic_q ( italic_R ) = 0 and q(ϵR)0𝑞italic-ϵ𝑅0q(\epsilon R)\neq 0italic_q ( italic_ϵ italic_R ) ≠ 0. This domain is transitional between the other two domains considered, and occurs for middling values of R𝑅Ritalic_R. In this instance, Eq. (19) can be used to refer f(R)𝑓𝑅f(R)italic_f ( italic_R ) to f(ϵR)𝑓italic-ϵ𝑅f(\epsilon R)italic_f ( italic_ϵ italic_R ). By construction, f(ϵR)𝑓italic-ϵ𝑅f(\epsilon R)italic_f ( italic_ϵ italic_R ) falls into the large-R𝑅Ritalic_R regime and can be solved by evaluating Eq. (20).

We show f(R)𝑓𝑅f(R)italic_f ( italic_R ) for a range of q0subscript𝑞0q_{0}italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT values in Fig. 5. We normalize the SFD to the small-R𝑅Ritalic_R limit p(R)/q0𝑝𝑅subscript𝑞0p(R)/q_{0}italic_p ( italic_R ) / italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT. We assume that β=1𝛽1\beta=1italic_β = 1, PS=0.1subscript𝑃𝑆0.1P_{S}=0.1italic_P start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT = 0.1 Pa, YS=104subscript𝑌𝑆superscript104Y_{S}=10^{-4}italic_Y start_POSTSUBSCRIPT italic_S end_POSTSUBSCRIPT = 10 start_POSTSUPERSCRIPT - 4 end_POSTSUPERSCRIPT, σt=50subscript𝜎𝑡50\sigma_{t}=50italic_σ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT = 50 Pa, ρ=500𝜌500\rho=500italic_ρ = 500 kg/m3, α=107𝛼superscript107\alpha=10^{-7}italic_α = 10 start_POSTSUPERSCRIPT - 7 end_POSTSUPERSCRIPT m2 s-1, and ϵ=21/3italic-ϵsuperscript213\epsilon=2^{1/3}italic_ϵ = 2 start_POSTSUPERSCRIPT 1 / 3 end_POSTSUPERSCRIPT, identically to Eq. (13). For small R𝑅Ritalic_R, the form of the SFD is unchanged from the input function (Eq. 18). For large R𝑅Ritalic_R, the relationship between Eq. (20) and p(R)𝑝𝑅p(R)italic_p ( italic_R ) depends on the removal parameter q0subscript𝑞0q_{0}italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT. If q0=0subscript𝑞00q_{0}=0italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT = 0, then f(R)Rp(R)proportional-to𝑓𝑅𝑅𝑝𝑅f(R)\propto Rp(R)italic_f ( italic_R ) ∝ italic_R italic_p ( italic_R ) and the power-law index is increased by a factor of order unity. Similarly, if q00subscript𝑞00q_{0}\neq 0italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ≠ 0 then f(R)p(R)proportional-to𝑓𝑅𝑝𝑅f(R)\propto p(R)italic_f ( italic_R ) ∝ italic_p ( italic_R ). As q00subscript𝑞00q_{0}\rightarrow 0italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT → 0, the small-n𝑛nitalic_n terms will be dominated by Rp(R)𝑅𝑝𝑅Rp(R)italic_R italic_p ( italic_R ), before eventually converging to p(R)𝑝𝑅p(R)italic_p ( italic_R ) at large R𝑅Ritalic_R. Therefore, the value of q0subscript𝑞0q_{0}italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT determines the slope of the power law for large objects. There are significantly more objects undergoing evolution than are injected, due to the relatively low values of the removal rate q0subscript𝑞0q_{0}italic_q start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT. The spikes at R1similar-to𝑅1R\sim 1italic_R ∼ 1 m correspond to the cutoff where objects are too small to be removed by division. At even smaller sizes there is another cutoff as the population is no longer fed by the division of larger objects.

Table 3: Dark Comet Source Probabilities. The probability that each dark comet originates from a given source based on our dynamical models. These results should be considered order of magnitude estimates, since (i) the source populations are normalized to the NEOMOD dataset and are not adjusted to our model’s modified orbital distributions and (ii) these probabilities do not account for the composition of the dark comets and the source populations, which may significantly modify estimations of their true origins. These source probabilities are generally consistent with the NEOMOD model (Nesvorný et al., 2024b) and the model of Granvik et al. (2018).
Object JFCs ν6subscript𝜈6\nu_{6}italic_ν start_POSTSUBSCRIPT 6 end_POSTSUBSCRIPT 3:1 5:2 7:3 8:3 9:4 11:5 2:1 Forest High i𝑖iitalic_i
1998 KY26 0 0.893 0.088 0.002 0 0.001 0 0.001 0 0.012 0.004
2005 VL1 0 0.850 0.125 0.010 0 0.001 0 0.001 0 0.012 0.001
2016 NJ33 0 0.831 0.096 0.049 0 0.002 0 0 0 0.020 0.003
2010 VL65 0 0.797 0.198 0 0 0 0 0 0 0.004 0.001
2010 RF12 0 0.948 0.026 0.010 0 0.002 0 0 0 0.010 0.004
2006 RH120 0 0.389 0.505 0 0 0 0 0 0 0.106 0
2003 RM 0.011 0.001 0.040 0.911 0.005 0.010 0.003 0.019 0 0 0
Refer to caption
Figure 6: The bin probabilities for JFCs and main belt objects (from the ν6subscript𝜈6\nu_{6}italic_ν start_POSTSUBSCRIPT 6 end_POSTSUBSCRIPT and 5:2 resonances) entering the near-Earth environment. The populations are initialized from the NEOMOD model (Nesvorný et al., 2023) and then integrated further with nongravitational accelerations. A white background indicates that no object reached that point in the simulations. Only samples with an inclination between 0superscript00^{\circ}0 start_POSTSUPERSCRIPT ∘ end_POSTSUPERSCRIPT and 15superscript1515^{\circ}15 start_POSTSUPERSCRIPT ∘ end_POSTSUPERSCRIPT are shown. The large white circle is 2003 RM, which has an orbit consistent with a JFC or 5:2 resonance origin. The white diamonds are the other, shorter-period dark comets, which are the most similar to a ν6subscript𝜈6\nu_{6}italic_ν start_POSTSUBSCRIPT 6 end_POSTSUBSCRIPT origin.

6 Dynamical Origins

In this section, we use the NEOMOD model (Nesvorný et al., 2023) to investigate the dynamical origins of the dark comets. While NEOMOD makes predictions for the source population for a given set of orbital parameters, this model may not be fully accurate for the dark comets. Principally, NEOMOD does not include nongravitational accelerations, which are the defining feature of dark comets.

As a result, we run dynamical simulations with nongravitational accelerations, using initial conditions from the NEOMOD model. We consider 11 of the 12 initial populations in NEOMOD — the Hungarias and Phocaeas are combined into a single high-inclination source. We therefore use the JFCs, the high-inclination sources, the ν6subscript𝜈6\nu_{6}italic_ν start_POSTSUBSCRIPT 6 end_POSTSUBSCRIPT resonance, the 3:1 resonance, the 5:2 resonance, the 7:3 resonance, the 8:3 resonance, the 9:4 resonance, the 11:5 resonance, the 2:1 resonance, and the forest of weak resonances in the inner main belt as individual sources.

For each source, we generate a representative initial sample of 104superscript10410^{4}10 start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT objects, which are chosen from the probability density function (PDF) given by the NEOMOD data. We set each object to have a nongravitational acceleration of the form described in Appendix A with a magnitude of |A0|=105subscript𝐴0superscript105|A_{0}|=10^{-5}| italic_A start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT | = 10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT au yr-2, consistent with the nongravitational accelerations of the dark comets (see Table 2). Each source population is then numerically integrated for 105superscript10510^{5}10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT yr, consistent with the mass budget described in Appendix B. Models are constructed for each source following the methodology described in Appendix C.

From this dataset, we calculate a PDF for comparison with the dark comets. We then normalize each source population by multiplying each source’s PDF by the fraction of all NEOs in that source (according to NEOMOD). This then provides the probability for each dark comet to originate from a given source. A more detailed discussion of this methodology is given in Appendix D. These probabilities are given in Table 3, although they are limited to order-of-magnitude estimates by several factors. First, our methodology of renormalization may not be accurate, since our population samples are not identical to those implemented in NEOMOD. Second, these probabilities do not account for the different degrees of volatile enrichment in these populations but only reflect our dynamical calculations. Since the dark comets are most likely outgassing, the volatile enrichment of these populations is an important factor to consider in identifying their origins and interpretation of these probabilities must be done with care. Future observations of the dark comets’ surface features and outgassing properties (species and rates) will constrain the source populations of the dark comets beyond the dynamics alone.

In Fig. 6, we show the PDF101010Note that this is not the probability that an object is from a source given its bin, but is the probability of an object being found in the bin given its source. The difference is subtle, but important. See Appendix D for more details. for various source populations as a function of semimajor axis and eccentricity, with the dark comet locations indicated. Specifically, we show PDFs for the JFCS, the ν6subscript𝜈6\nu_{6}italic_ν start_POSTSUBSCRIPT 6 end_POSTSUBSCRIPT resonance, and the 5:2 resonance (see Appendix C for more complete source region calculations). We restrict the PDF to an inclination of i[0,15]𝑖superscript015i\in[0,15]^{\circ}italic_i ∈ [ 0 , 15 ] start_POSTSUPERSCRIPT ∘ end_POSTSUPERSCRIPT, which includes all of the dark comets (see Table 2).111111Note that the analogous calculations in Appendix C do not include an inclination cutoff. 2003 RM is most likely a main belt object, specifically from the nearby 5:2 resonance, but is also located where a JFC is most likely to occur. At this time, it appears that 2003 RM is volatile enriched and could plausibly originate in either the JFC region or the main belt, presumably via an MBC. Follow-up observations may help to constrain the origin of this object.

The remaining dark comet orbits are consistent with a main belt source, particularly the ν6subscript𝜈6\nu_{6}italic_ν start_POSTSUBSCRIPT 6 end_POSTSUBSCRIPT source in the inner main belt. However, the retention of volatiles in the inner main belt is more challenging, suggesting the presence of undiscovered MBCs or subsurface volatiles on asteroids in this region. If the dark comets are produced by a rotational fragmentation cascade, then a single object en route to a year-period orbit may produce several orders of magnitude more small objects, increasing the prevalence of this population and accounting for the number of dark comets.

These simulations show a significant mixing of the NEO source populations when nongravitational accelerations are incorporated. This makes it challenging to identify a dynamical origin based on orbital elements alone. Other characteristics — color, albedo, nongravitational accelerations, material strength, etc. — will be necessary to clarify any given object’s dynamical origin. Given these results, the inner dark comets likely originated from the main belt. While 2003 RM’s dynamics imply an origin in the outer main belt, a JFC origin is also plausible. Given the apparent degeneracies in source population probabilities, follow-up observations of physical properties may aid in the identification of the dark comets’ source populations.

Refer to caption
Figure 7: Sizes and semimajor axes of dark comets, the JFCs 252P/LINEAR and 460P/PANSTARRS, the MBCs 133P/Elst-Pizarro and 238P/Read, and the possible dead comets. The color indicates the eccentricity.

7 Discussion

In this paper, we investigate the physical and dynamical origins of the “dark comets” identified by Farnocchia et al. (2023) and Seligman et al. (2023). We propose that these objects begin as smaller, strength-dominated objects and evolve through gravitational effects into the NEO environment. During this transition, and especially once in the NEO environment, these progenitor objects undergo significant nongravitational accelerations that cause a rotational fragmentation cascade. The resulting fragments would continue to spin up and devolatilize. However, the rotational torque is reduced as the volatiles are depleted. Eventually, the torque will be insufficient to spin the object up beyond its critical rate and it will be stable against further rotational destruction. At this point, the fragments will be small, rapidly rotating, and mostly devolatilized.

While small-number statistics make comparison difficult, there are several lines of observational evidence that suggest that the dark comets may be the product of a rotational fragmentation cascade. These objects have small radii, weak nongravitational accelerations, and rapid rotation rates (where measured). Notably, these rotation periods are above the stability threshold for larger objects. Therefore, these objects may have been generated by the fragmentation of a larger object. We also calculate the SFD of the dark comets under the action of the rotational fragmentation cascade. The SFD has two observable signatures: (i) a modification to the power-law slope of the injected objects’ SFD and (ii) a sharp increase in the number of objects with sizes of R1similar-to𝑅1R\sim 1italic_R ∼ 110101010 m. The second of these signatures is consistent with the currently known population of dark comets, although limited by small-number statistics.

This model may also explain the origins of several NEOs with low densities and detected nongravitational acceleration due to radiation pressure — 2009 BD, 2012 LA, and 2011 MD. These objects are small (R10similar-to𝑅10R\sim 10italic_R ∼ 10 m) and have densities of ρ(500±300)similar-to-or-equals𝜌plus-or-minus500300\rho\simeq(500\pm 300)italic_ρ ≃ ( 500 ± 300 ) kg/m3, consistent with comets (Richardson et al., 2007; Micheli et al., 2012, 2013; Micheli, 2013; Micheli et al., 2014; Pätzold et al., 2016). Mommert et al. (2014a) and Mommert et al. (2014b) found larger densities consistent with asteroids for 2009 BD and 2011 MD, although 2011 MD’s density is still consistent with 500 kg/m3. The extinct devolatilized comet fragments predicted in this model should be meter-scale with relatively low densities, consistent with these NEOs. These bodies, as well as the possible precursor JFCs 252P/LINEAR and 460P/PANSTARRS and volatile-active MBCs 133P/Elst-Pizarro and 238P/Read, are shown in Fig. 7 along with the dark comets.

The rotational fragmentation model presented here both is self-consistent and explains the observed properties of the dark comets. However, there are significant sources of uncertainties in the current data and rotation periods are only measured on a few dark comets. As a result, alternative models may be similarly consistent with the limited data. For example, tidal disruption may also be responsible for the creation of the dark comets. A larger rubble-pile comet could undergo tidal disruption through a close encounter with a terrestrial planet, producing a family of smaller NEOs with similar sizes to the dark comets (Granvik & Walsh, 2024; Nesvorný et al., 2024a). These events will produce an overdensity of meter-scale objects, similar to our predicted rotational-fragmentation SFD (Sec. 5). Both tidal and rotational disruption may therefore be responsible for producing dark comets, depending on the relative size and binding force of the progenitor object.

The rotational fragmentation evolutionary track predicts a relatively large population of devolatilized fragments in the NEO environment, either with undetectably low nongravitational accelerations (dark comets) or fully extinct (dead comets). However, this calculation relies on the residence time and mass-loss time of these objects, which must be adjusted for the presence of nongravitational accelerations. Moreover, our estimates only apply to small NEOs that are capable of being produced by this model (Sec. 5). While this prediction is consistent with our estimate that 0.5 – 60% of NEOs could originate along this pathway, more work is needed to refine this value.

If these objects are generated by a rotational fragmentation cascade, few objects should be observed in the interim stage between the large, relatively young progenitors and the rapidly-rotating, small, thermally evolved inner dark comets. The rotational cascade timescale is 103similar-toabsentsuperscript103\sim 10^{3}∼ 10 start_POSTSUPERSCRIPT 3 end_POSTSUPERSCRIPT yr, significantly shorter than the residence time in the near-Earth environment of 106similar-toabsentsuperscript106\sim 10^{6}∼ 10 start_POSTSUPERSCRIPT 6 end_POSTSUPERSCRIPT yr (Gladman et al., 2000; Nesvorný et al., 2023). There should therefore be two distinct populations of dark comets — the progenitors and the fragments. The progenitors are relatively large (100 – 1000 m), new to the near-Earth environment, and may have only recently begun to sublimate, not yet triggering the rotational fragmentation cascade. Meanwhile, the small fragments have likely finished this cascade and are mostly devolatilized. Since these evolutionary stages have long lifetimes compared to the rotational fragmentation cascade, most of the dark comets will be observed in these two edge states.

This model also explains the outlier of 2003 RM, which is larger and has a primarily transverse nongravitational acceleration. In contrast, the smaller dark comets primarily have out-of-plane nongravitational components. This feature (although limited by small-number statistics) can be understood in the context of this model. Volatile-rich objects will likely have radial and transverse nongravitational accelerations when they first become active. These accelerations may spin up the objects and trigger a rotational fragmentation cascade. When the nuclei are rotating sufficiently rapidly to have longitudinally isothermal surfaces, the outgassing can align with the spin axis and produce out-of-plane accelerations (Taylor et al., 2024). The size of 2003 RM and the failure of the outgassing balancing mechanism to reproduce its acceleration are therefore consistent with this model, implying that 2003 RM is in the beginning stage of this evolutionary track and has not been spun-up to the isothermal limit.

We also investigate the dynamical origins of these objects. The orbit of 2003 RM implies that it likely originated in the outer main belt (see Table 3), although there are several caveats to this conclusion. First, the orbit of 2003 RM resides precisely at the maximum probability for a JFC origin (see Fig. 6). Second, the relatively low probability of a JFC origin compared to a main belt origin primarily reflects the small number of JFCs relative to main belt objects (Table 4). Finally, if the nongravitational acceleration of 2003 RM is driven by outgassing of volatiles, then the volatile-rich JFCs are a more likely a priori source. As a result, the JFCs are still a plausible source for 2003 RM.

Meanwhile, our numerical experiments demonstrate that the smaller dark comets likely originated in the inner main belt. This may have significant ramifications concerning the existence and abundance of volatiles in this region. While it is already understood that subsurface volatiles can exist in the inner main belt (Fanale & Salvail, 1989; Schörghofer, 2008) and survive into the near-Earth environment (Schörghofer et al., 2020), confirmation of this volatile reservoir remains elusive. The dark comets potentially provide evidence that a volatile reservoir exists in the inner main belt. Furthermore, the presence of these objects in the near-Earth environment potentially provides an additional pathway for terrestrial volatile delivery.

Ground- and space-based follow-up observations of dark comets may enable measurement of the outgassing rates and compositions, potentially constraining their dynamical origins. Spectral observations in particular may provide more information regarding their progenitor asteroid spectral type and the size of the main belt volatile reservoir. Future survey missions such as NEO Surveyor, the Rubin Observatory Legacy Survey of Space and Time (LSST), and the Gaia mission may also identify more dark comets than are currently known, which will refine our understanding of this evolutionary pathway and their source populations.

Finally, the Hayabusa2 Extended Mission (Hayabusa2#) will arrive at 1998 KY26 in 2031 (Hirabayashi et al., 2021; Kikuchi et al., 2023). If the rotational fragmentation model is correct, then Hayabusa2# should find an object with high porosity and a material strength of σt5.3subscript𝜎𝑡5.3\sigma_{t}\geq 5.3italic_σ start_POSTSUBSCRIPT italic_t end_POSTSUBSCRIPT ≥ 5.3 Pa, similar to a cometary body. In addition, the rotation period of 1998 KY26 was last measured in 1999 (Ostro et al., 1999). As a result, the SYORP effect may have further spun up 1998 KY26 in the intervening time, which would be detectable by Hayabusa2# and further radar observations. This will provide another constraint on the accuracy of the rotational fragmentation model. Given the dynamical origins and nongravitational acceleration of 1998 KY26, this mission will also provide an opportunity to further investigate and characterize the unique volatile reservoir in the inner main belt.

Acknowledgements

We thank the two anonymous reviewers for their helpful suggestions, which significantly improved the scientific content of this manuscript. We also thank the editor of this manuscript, Lori Feaga, for her assistance in the review process. We thank Fred Adams, Juliette Becker, Bill Chen, Fei Dai, Adina Feinstein, Thomas Kennedy, Garrett Levine, Nikole Lewis, Kevin Napier, Luis Salazar Manzano, Cindy Xiang, and Andrew Youdin for useful conversations and suggestions. A.G.T. thanks Fiona Corcoran for assistance with figures. A.G.T. acknowledges support from the Fannie and John Hertz Foundation and the University of Michigan’s Rackham Merit Fellowship Program. J.K.S. acknowledges support from NASA Grant No. 80NSSC19K1313 and NASA Grant No. 80NSSC22K1399. D.Z.S. is supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2202135. This research award is partially funded by a generous gift of Charles Simonyi to the NSF Division of Astronomical Sciences. The award is made in recognition of significant contributions to Rubin Observatory’s Legacy Survey of Space and Time. D.F. conducted this research at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). L.D. acknowledges support from an internal research grant at Southwest Research Institute.

This paper made use of the Julia programming language (Bezanson et al., 2017) and the plotting package Makie (Danisch & Krumbiegel, 2021). Simulations in this paper made use of the REBOUND N-body code (Rein & Liu, 2012) and the REBOUNDx library (Tamayo et al., 2020). The simulations were integrated using MERCURIUS (Rein et al., 2019).

References

  • A’Hearn et al. (2005) A’Hearn, M. F., Belton, M. J. S., Delamere, W. A., et al. 2005, Science, 310, 258, doi: 10.1126/science.1118923
  • A’Hearn et al. (2011) A’Hearn, M. F., Belton, M. J. S., Delamere, W. A., et al. 2011, Science, 332, 1396, doi: 10.1126/science.1204054
  • Asher et al. (1994) Asher, D. J., Bailey, M. E., Hahn, G., & Steel, D. I. 1994, MNRAS, 267, 26, doi: 10.1093/mnras/267.1.26
  • Asphaug & Benz (1996) Asphaug, E., & Benz, W. 1996, Icarus, 121, 225, doi: 10.1006/icar.1996.0083
  • Attree et al. (2018) Attree, N., Groussin, O., Jorda, L., et al. 2018, A&A, 611, A33, doi: 10.1051/0004-6361/201732155
  • Bertini (2011) Bertini, I. 2011, Planet. Space Sci., 59, 365, doi: 10.1016/j.pss.2011.01.014
  • Bezanson et al. (2017) Bezanson, J., Edelman, A., Karpinski, S., & Shah, V. B. 2017, SIAM review, 59, 65. https://doi.org/10.1137/141000671
  • Bottke et al. (2005) Bottke, W. F., Durda, D. D., Nesvorný, D., et al. 2005, Icarus, 179, 63, doi: 10.1016/j.icarus.2005.05.017
  • Bottke et al. (2023) Bottke, W. F., Vokrouhlický, D., Marshall, R., et al. 2023, \psj, 4, 168, doi: 10.3847/PSJ/ace7cd
  • Bowling et al. (2014) Bowling, T. J., Steckloff, J., Graves, K., & Melosh, H. J. 2014, in AAS/Division for Planetary Sciences Meeting Abstracts #46, 100.03
  • Brasser & Wang (2015) Brasser, R., & Wang, J. H. 2015, A&A, 573, A102, doi: 10.1051/0004-6361/201423687
  • Breiter et al. (2012) Breiter, S., RoŻek, A., & Vokrouhlický, D. 2012, MNRAS, 427, 755, doi: 10.1111/j.1365-2966.2012.21970.x
  • Brownlee et al. (2006) Brownlee, D., Tsou, P., Aléon, J., et al. 2006, Science, 314, 1711, doi: 10.1126/science.1135840
  • Carry (2012) Carry, B. 2012, Planet. Space Sci., 73, 98, doi: 10.1016/j.pss.2012.03.009
  • Chesley et al. (2016) Chesley, S. R., Farnocchia, D., Pravec, P., & Vokrouhlický, D. 2016, in Asteroids: New Observations, New Models, ed. S. R. Chesley, A. Morbidelli, R. Jedicke, & D. Farnocchia, Vol. 318, 250–258, doi: 10.1017/S1743921315008790
  • Danisch & Krumbiegel (2021) Danisch, S., & Krumbiegel, J. 2021, Journal of Open Source Software, 6, 3349, doi: 10.21105/joss.03349
  • Davidsson (1999) Davidsson, B. J. R. 1999, Icarus, 142, 525, doi: 10.1006/icar.1999.6214
  • Davidsson (2001) —. 2001, Icarus, 149, 375, doi: 10.1006/icar.2000.6540
  • Di Sisto & Brunini (2007) Di Sisto, R. P., & Brunini, A. 2007, Icarus, 190, 224, doi: 10.1016/j.icarus.2007.02.012
  • Di Sisto et al. (2009) Di Sisto, R. P., Fernández, J. A., & Brunini, A. 2009, Icarus, 203, 140, doi: 10.1016/j.icarus.2009.05.002
  • Fanale & Salvail (1989) Fanale, F. P., & Salvail, J. R. 1989, Icarus, 82, 97, doi: 10.1016/0019-1035(89)90026-2
  • Farinella et al. (1994) Farinella, P., Froeschlé, C., Froeschlé, C., et al. 1994, Nature, 371, 314, doi: 10.1038/371314a0
  • Farnocchia et al. (2015) Farnocchia, D., Chesley, S. R., Milani, A., Gronchi, G. F., & Chodas, P. W. 2015, in Asteroids IV (University of Arizona Press), 815–834, doi: 10.2458/azu_uapress_9780816532131-ch041
  • Farnocchia et al. (2013) Farnocchia, D., Chesley, S. R., Vokrouhlický, D., et al. 2013, Icarus, 224, 1, doi: 10.1016/j.icarus.2013.02.004
  • Farnocchia et al. (2023) Farnocchia, D., Seligman, D. Z., Granvik, M., et al. 2023, \psj, 4, 29, doi: 10.3847/PSJ/acb25b
  • Ferellec et al. (2023) Ferellec, L., Snodgrass, C., Fitzsimmons, A., et al. 2023, MNRAS, 518, 2373, doi: 10.1093/mnras/stac3199
  • Fernández et al. (2002) Fernández, J. A., Gallardo, T., & Brunini, A. 2002, Icarus, 159, 358, doi: 10.1006/icar.2002.6903
  • Gladman et al. (2000) Gladman, B., Michel, P., & Froeschlé, C. 2000, Icarus, 146, 176, doi: 10.1006/icar.2000.6391
  • Gladman et al. (1997) Gladman, B. J., Migliorini, F., Morbidelli, A., et al. 1997, Science, 277, 197, doi: 10.1126/science.277.5323.197
  • Granvik et al. (2017) Granvik, M., Morbidelli, A., Vokrouhlický, D., et al. 2017, A&A, 598, A52, doi: 10.1051/0004-6361/201629252
  • Granvik & Walsh (2024) Granvik, M., & Walsh, K. J. 2024, ApJ, 960, L9, doi: 10.3847/2041-8213/ad151b
  • Granvik et al. (2018) Granvik, M., Morbidelli, A., Jedicke, R., et al. 2018, Icarus, 312, 181, doi: 10.1016/j.icarus.2018.04.018
  • Grav et al. (2023) Grav, T., Mainzer, A. K., Masiero, J. R., et al. 2023, \psj, 4, 228, doi: 10.3847/PSJ/ad072e
  • Greenberg et al. (2020) Greenberg, A. H., Margot, J.-L., Verma, A. K., Taylor, P. A., & Hodge, S. E. 2020, AJ, 159, 92, doi: 10.3847/1538-3881/ab62a3
  • Groussin et al. (2019) Groussin, O., Attree, N., Brouet, Y., et al. 2019, Space Sci. Rev., 215, 29, doi: 10.1007/s11214-019-0594-x
  • Gundlach & Blum (2012) Gundlach, B., & Blum, J. 2012, Icarus, 219, 618, doi: https://doi.org/10.1016/j.icarus.2012.03.013
  • Hirabayashi et al. (2016) Hirabayashi, M., Scheeres, D. J., Chesley, S. R., et al. 2016, Nature, 534, 352, doi: 10.1038/nature17670
  • Hirabayashi et al. (2021) Hirabayashi, M., Mimasu, Y., Sakatani, N., et al. 2021, Advances in Space Research, 68, 1533, doi: 10.1016/j.asr.2021.03.030
  • Hsieh (2017) Hsieh, H. H. 2017, Philosophical Transactions of the Royal Society of London Series A, 375, 20160259, doi: 10.1098/rsta.2016.0259
  • Hsieh & Jewitt (2006) Hsieh, H. H., & Jewitt, D. 2006, Science, 312, 561, doi: 10.1126/science.1125150
  • Jewitt (1997) Jewitt, D. 1997, Earth Moon and Planets, 79, 35, doi: 10.1023/A:1006272914117
  • Jewitt (2005) —. 2005, AJ, 129, 530, doi: 10.1086/426328
  • Jewitt (2009) —. 2009, AJ, 137, 4296, doi: 10.1088/0004-6256/137/5/4296
  • Jewitt (2012) —. 2012, AJ, 143, 66, doi: 10.1088/0004-6256/143/3/66
  • Jewitt (2021) —. 2021, AJ, 161, 261, doi: 10.3847/1538-3881/abf09c
  • Jewitt (2022) —. 2022, AJ, 164, 158, doi: 10.3847/1538-3881/ac886d
  • Jewitt & Hsieh (2022) Jewitt, D., & Hsieh, H. H. 2022, arXiv e-prints, arXiv:2203.01397, doi: 10.48550/arXiv.2203.01397
  • Jewitt et al. (2017) Jewitt, D., Luu, J., Rajagopal, J., et al. 2017, ApJL, 850, L36, doi: 10.3847/2041-8213/aa9b2f
  • Jorda et al. (2016) Jorda, L., Gaskell, R., Capanna, C., et al. 2016, Icarus, 277, 257, doi: 10.1016/j.icarus.2016.05.002
  • Kikuchi et al. (2023) Kikuchi, S., Mimasu, Y., Takei, Y., et al. 2023, Acta Astronautica, 211, 295, doi: https://doi.org/10.1016/j.actaastro.2023.06.010
  • Knight et al. (2023) Knight, M. M., Kokotanekova, R., & Samarasinha, N. H. 2023, arXiv e-prints, arXiv:2304.09309, doi: 10.48550/arXiv.2304.09309
  • Kokotanekova et al. (2017) Kokotanekova, R., Snodgrass, C., Lacerda, P., et al. 2017, Monthly Notices of the Royal Astronomical Society, 471, 2974, doi: 10.1093/mnras/stx1716
  • Kowal et al. (1979) Kowal, C. T., Liller, W., & Marsden, B. G. 1979, in Dynamics of the Solar System, ed. R. L. Duncombe, Vol. 81, 245
  • Kwiatkowski et al. (2009) Kwiatkowski, T., Kryszczyńska, A., Polińska, M., et al. 2009, A&A, 495, 967, doi: 10.1051/0004-6361:200810965
  • Kwon et al. (2022) Kwon, Y. G., Masiero, J. R., & Markkanen, J. 2022, A&A, 668, A97, doi: 10.1051/0004-6361/202244853
  • Levison & Duncan (1997) Levison, H. F., & Duncan, M. J. 1997, Icarus, 127, 13, doi: 10.1006/icar.1996.5637
  • Li et al. (2017) Li, J.-Y., Kelley, M. S. P., Samarasinha, N. H., et al. 2017, AJ, 154, 136, doi: 10.3847/1538-3881/aa86ae
  • Licandro et al. (2018) Licandro, J., Popescu, M., de León, J., et al. 2018, A&A, 618, A170, doi: 10.1051/0004-6361/201832853
  • Lisse et al. (2022) Lisse, C. M., Gladstone, G. R., Young, L. A., et al. 2022, \psj, 3, 112, doi: 10.3847/PSJ/ac6097
  • Marsden (1969) Marsden, B. G. 1969, AJ, 74, 720, doi: 10.1086/110848
  • Meech et al. (2016) Meech, K. J., Yang, B., Kleyna, J., et al. 2016, Science Advances, 2, e1600038, doi: 10.1126/sciadv.1600038
  • Micheli (2013) Micheli, M. 2013, PhD thesis, University of Hawaii, Manoa
  • Micheli et al. (2012) Micheli, M., Tholen, D. J., & Elliott, G. T. 2012, New A, 17, 446, doi: 10.1016/j.newast.2011.11.008
  • Micheli et al. (2013) —. 2013, Icarus, 226, 251, doi: 10.1016/j.icarus.2013.05.032
  • Micheli et al. (2014) —. 2014, ApJ, 788, L1, doi: 10.1088/2041-8205/788/1/L1
  • Migliorini et al. (1998) Migliorini, F., Michel, P., Morbidelli, A., Nesvorný, D., & Zappala, V. 1998, Science, 281, 2022, doi: 10.1126/science.281.5385.2022
  • Mommert et al. (2014a) Mommert, M., Hora, J. L., Farnocchia, D., et al. 2014a, ApJ, 786, 148, doi: 10.1088/0004-637X/786/2/148
  • Mommert et al. (2014b) Mommert, M., Farnocchia, D., Hora, J. L., et al. 2014b, ApJ, 789, L22, doi: 10.1088/2041-8205/789/1/L22
  • Naidu et al. (2016) Naidu, S. P., Benner, L. A. M., Brozovic, M., et al. 2016, in AAS/Division for Planetary Sciences Meeting Abstracts #48, 219.05
  • Nesvorný (2018) Nesvorný, D. 2018, ARA&A, 56, 137, doi: 10.1146/annurev-astro-081817-052028
  • Nesvorný et al. (2017) Nesvorný, D., Vokrouhlický, D., Dones, L., et al. 2017, ApJ, 845, 27, doi: 10.3847/1538-4357/aa7cf6
  • Nesvorný et al. (2023) Nesvorný, D., Deienno, R., Bottke, W. F., et al. 2023, AJ, 166, 55, doi: 10.3847/1538-3881/ace040
  • Nesvorný et al. (2024a) Nesvorný, D., Vokrouhlický, D., Shelly, F., et al. 2024a, Icarus, 411, 115922, doi: 10.1016/j.icarus.2023.115922
  • Nesvorný et al. (2024b) —. 2024b, Icarus, 417, 116110, doi: 10.1016/j.icarus.2024.116110
  • Ostro et al. (1999) Ostro, S. J., Pravec, P., Benner, L. A. M., et al. 1999, Science, 285, 557, doi: 10.1126/science.285.5427.557
  • Pätzold et al. (2016) Pätzold, M., Andert, T., Hahn, M., et al. 2016, Nature, 530, 63, doi: 10.1038/nature16535
  • Piro et al. (2021) Piro, C., Meech, K. J., Bufanda, E., et al. 2021, \psj, 2, 33, doi: 10.3847/PSJ/abd552
  • Pravec & Harris (2000) Pravec, P., & Harris, A. W. 2000, Icarus, 148, 12, doi: 10.1006/icar.2000.6482
  • Pravec & Harris (2007) —. 2007, Icarus, 190, 250, doi: 10.1016/j.icarus.2007.02.023
  • Pravec et al. (2002) Pravec, P., Harris, A. W., & Michalowski, T. 2002, in Asteroids III (University of Arizona Press), 113–122
  • Rein & Liu (2012) Rein, H., & Liu, S. F. 2012, A&A, 537, A128, doi: 10.1051/0004-6361/201118085
  • Rein et al. (2019) Rein, H., Hernandez, D. M., Tamayo, D., et al. 2019, MNRAS, 485, 5490, doi: 10.1093/mnras/stz769
  • Richardson et al. (2007) Richardson, J. E., Melosh, H. J., Lisse, C. M., & Carcich, B. 2007, Icarus, 190, 357, doi: 10.1016/j.icarus.2007.08.001
  • Rozitis et al. (2014) Rozitis, B., Maclennan, E., & Emery, J. P. 2014, Nature, 512, 174, doi: 10.1038/nature13632
  • Safrit et al. (2021) Safrit, T. K., Steckloff, J. K., Bosh, A. S., et al. 2021, \psj, 2, 14, doi: 10.3847/PSJ/abc9c8
  • Sánchez & Scheeres (2014) Sánchez, P., & Scheeres, D. J. 2014, \maps, 49, 788, doi: 10.1111/maps.12293
  • Schörghofer (2008) Schörghofer, N. 2008, ApJ, 682, 697, doi: 10.1086/588633
  • Schörghofer et al. (2020) Schörghofer, N., Hsieh, H. H., Novaković, B., & Walsh, K. J. 2020, Icarus, 348, 113865, doi: 10.1016/j.icarus.2020.113865
  • Sekanina & Yeomans (1985) Sekanina, Z., & Yeomans, D. K. 1985, AJ, 90, 2335, doi: 10.1086/113939
  • Seligman et al. (2021) Seligman, D. Z., Kratter, K. M., Levine, W. G., & Jedicke, R. 2021, \psj, 2, 234, doi: 10.3847/PSJ/ac2dee
  • Seligman et al. (2023) Seligman, D. Z., Farnocchia, D., Micheli, M., et al. 2023, \psj, 4, 35, doi: 10.3847/PSJ/acb697
  • Shannon et al. (2015) Shannon, A., Jackson, A. P., Veras, D., & Wyatt, M. 2015, MNRAS, 446, 2059, doi: 10.1093/mnras/stu2267
  • Soderblom et al. (2002) Soderblom, L. A., Becker, T. L., Bennett, G., et al. 2002, Science, 296, 1087, doi: 10.1126/science.1069527
  • Sonnett et al. (2011) Sonnett, S., Kleyna, J., Jedicke, R., & Masiero, J. 2011, Icarus, 215, 534, doi: 10.1016/j.icarus.2011.08.001
  • Steckloff & Jacobson (2016) Steckloff, J. K., & Jacobson, S. A. 2016, Icarus, 264, 160, doi: 10.1016/j.icarus.2015.09.021
  • Steckloff et al. (2015) Steckloff, J. K., Johnson, B. C., Bowling, T., et al. 2015, Icarus, 258, 430, doi: https://doi.org/10.1016/j.icarus.2015.06.032
  • Steckloff et al. (2021) Steckloff, J. K., Lisse, C. M., Safrit, T. K., et al. 2021, Icarus, 356, 113998, doi: 10.1016/j.icarus.2020.113998
  • Steckloff & Samarasinha (2018) Steckloff, J. K., & Samarasinha, N. H. 2018, Icarus, 312, 172, doi: 10.1016/j.icarus.2018.04.031
  • Tamayo et al. (2020) Tamayo, D., Rein, H., Shi, P., & Hernandez, D. M. 2020, MNRAS, 491, 2885, doi: 10.1093/mnras/stz2870
  • Taylor et al. (2024) Taylor, A. G., Farnocchia, D., Vokrouhlický, D., et al. 2024, Icarus, 408, 115822, doi: 10.1016/j.icarus.2023.115822
  • Tiscareno & Malhotra (2003) Tiscareno, M. S., & Malhotra, R. 2003, AJ, 126, 3122, doi: 10.1086/379554
  • Vincent et al. (2015) Vincent, J.-B., Oklay, N., Marchi, S., Höfner, S., & Sierks, H. 2015, Planet. Space Sci., 107, 53, doi: 10.1016/j.pss.2014.06.008
  • Vokrouhlický et al. (2015) Vokrouhlický, D., Bottke, W. F., Chesley, S. R., Scheeres, D. J., & Statler, T. S. 2015, in Asteroids IV (University of Arizona Press), 509–531, doi: 10.2458/azu_uapress_9780816532131-ch027
  • Vokrouhlický & Milani (2000) Vokrouhlický, D., & Milani, A. 2000, A&A, 362, 746
  • Volk & Malhotra (2008) Volk, K., & Malhotra, R. 2008, ApJ, 687, 714, doi: 10.1086/591839
  • Warner et al. (2009) Warner, B. D., Harris, A. W., & Pravec, P. 2009, Icarus, 202, 134, doi: 10.1016/j.icarus.2009.02.003
  • Weissman & Levison (1997) Weissman, P. R., & Levison, H. F. 1997, ApJ, 488, L133, doi: 10.1086/310940
  • Whipple (1950) Whipple, F. L. 1950, ApJ, 111, 375, doi: 10.1086/145272
  • Whipple (1951) —. 1951, ApJ, 113, 464, doi: 10.1086/145416
  • Wiegert & Tremaine (1999) Wiegert, P., & Tremaine, S. 1999, Icarus, 137, 84, doi: 10.1006/icar.1998.6040
  • Yeomans et al. (2004) Yeomans, D. K., Chodas, P. W., Sitarski, G., Szutowicz, S., & Królikowska, M. 2004, in Comets II, ed. M. C. Festou, H. U. Keller, & H. A. Weaver (University of Arizona Press), 137
  • Yeomans et al. (2000) Yeomans, D. K., Antreasian, P. G., Barriot, J. P., et al. 2000, Science, 289, 2085, doi: 10.1126/science.289.5487.2085
  • Zhang et al. (2021) Zhang, Y., Michel, P., Richardson, D. C., et al. 2021, Icarus, 362, 114433, doi: https://doi.org/10.1016/j.icarus.2021.114433
  • Zhang et al. (2017) Zhang, Y., Richardson, D. C., Barnouin, O. S., et al. 2017, Icarus, 294, 98, doi: https://doi.org/10.1016/j.icarus.2017.04.027
\restartappendixnumbering

Appendix A Stochastic Nongravitational Acceleration

In this section, we introduce the model of nongravitational acceleration that we use in our dynamical simulations. Our model is loosely based on the model of Fernández et al. (2002), who showed that nongravitational accelerations may allow objects to decouple from Jupiter and migrate inwards. We assume that the nongravitational acceleration takes the form

𝑭=αA0(r0r)2𝒆^v.𝑭𝛼subscript𝐴0superscriptsubscript𝑟0𝑟2subscriptbold-^𝒆𝑣{\boldsymbol{F}=\alpha A_{0}\left(\frac{r_{0}}{r}\right)^{2}\boldsymbol{\hat{e% }}_{v}}\,.bold_italic_F = italic_α italic_A start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT ( divide start_ARG italic_r start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_r end_ARG ) start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT overbold_^ start_ARG bold_italic_e end_ARG start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT . (21)

In Eq. (21), 𝒆^vsubscriptbold-^𝒆𝑣\boldsymbol{\hat{e}}_{v}overbold_^ start_ARG bold_italic_e end_ARG start_POSTSUBSCRIPT italic_v end_POSTSUBSCRIPT is the unit velocity vector, r0subscript𝑟0r_{0}italic_r start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT is a scaling distance (typically 1 au), r𝑟ritalic_r is the heliocentric distance, and α𝛼\alphaitalic_α is a constant between -1 and 1. The parameter α𝛼\alphaitalic_α is introduced in order to account for stochastic variation in the direction and magnitude of the nongravitational acceleration. Here, this parameter is randomly chosen from a sinusoidal distribution and resampled every of 50 yr. We have implemented this nongravitational acceleration in REBOUNDx (Tamayo et al., 2020).

Appendix B Outgassing Mass Budget

In this section, we calculate a mass budget for the known dark comets, which we convert into an allowable nongravitational acceleration magnitude over a timescale.

The mass-loss rate by outgassing of a species X𝑋Xitalic_X is given by (Seligman et al., 2023)

dMdt=M|Ai|vgasζ.d𝑀d𝑡𝑀subscript𝐴𝑖subscript𝑣gas𝜁{\frac{{\rm d}M}{{\rm d}t}=\frac{M|A_{i}|}{v_{\rm gas}\zeta}}\,.divide start_ARG roman_d italic_M end_ARG start_ARG roman_d italic_t end_ARG = divide start_ARG italic_M | italic_A start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | end_ARG start_ARG italic_v start_POSTSUBSCRIPT roman_gas end_POSTSUBSCRIPT italic_ζ end_ARG . (22)

The variable ζ𝜁\zetaitalic_ζ indicates the collimation of the outgassing, and vgassubscript𝑣gasv_{\rm gas}italic_v start_POSTSUBSCRIPT roman_gas end_POSTSUBSCRIPT is the gas velocity, given by

vgas=(8kBTgasπmX)1/2.subscript𝑣gassuperscript8subscript𝑘𝐵subscript𝑇gas𝜋subscript𝑚𝑋12{v_{\rm gas}=\left(\frac{8k_{B}T_{\rm gas}}{\pi m_{X}}\right)^{1/2}}\,.italic_v start_POSTSUBSCRIPT roman_gas end_POSTSUBSCRIPT = ( divide start_ARG 8 italic_k start_POSTSUBSCRIPT italic_B end_POSTSUBSCRIPT italic_T start_POSTSUBSCRIPT roman_gas end_POSTSUBSCRIPT end_ARG start_ARG italic_π italic_m start_POSTSUBSCRIPT italic_X end_POSTSUBSCRIPT end_ARG ) start_POSTSUPERSCRIPT 1 / 2 end_POSTSUPERSCRIPT . (23)

Assuming that Aisubscript𝐴𝑖A_{i}italic_A start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT and vgassubscript𝑣gasv_{\rm gas}italic_v start_POSTSUBSCRIPT roman_gas end_POSTSUBSCRIPT are constants in time, Eq. (22) can be solved to find that

M(t)=M0exp(|Ai|tvgasζ).𝑀𝑡subscript𝑀0subscript𝐴𝑖𝑡subscript𝑣gas𝜁{M(t)=M_{0}\exp\left(-\frac{|A_{i}|t}{v_{\rm gas}\zeta}\right)}\,.italic_M ( italic_t ) = italic_M start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT roman_exp ( - divide start_ARG | italic_A start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | italic_t end_ARG start_ARG italic_v start_POSTSUBSCRIPT roman_gas end_POSTSUBSCRIPT italic_ζ end_ARG ) . (24)

Technically, there is no true time limit for a nongravitational acceleration with a fixed magnitude — as the object’s mass decreases, the acceleration requires less mass outflow, and so M0𝑀0M\rightarrow 0italic_M → 0 only when t𝑡t\rightarrow\inftyitalic_t → ∞. However, we are here restricted by our (assumed) knowledge of the initial and final states of the dark comets — objects of Rnuc200similar-tosubscript𝑅nuc200R_{\rm nuc}\sim 200italic_R start_POSTSUBSCRIPT roman_nuc end_POSTSUBSCRIPT ∼ 200 m become objects of Rnuc10similar-tosubscript𝑅nuc10R_{\rm nuc}\sim 10italic_R start_POSTSUBSCRIPT roman_nuc end_POSTSUBSCRIPT ∼ 10 m. Assuming that the dark comets are spheres, the acceleration and the timescale are therefore constrained by

|Ai|τ3vgasζln(R0R(τ)).subscript𝐴𝑖𝜏3subscript𝑣gas𝜁subscript𝑅0𝑅𝜏|A_{i}|\tau\leq 3v_{\rm gas}\zeta\ln\left(\frac{R_{0}}{R(\tau)}\right)\,.| italic_A start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | italic_τ ≤ 3 italic_v start_POSTSUBSCRIPT roman_gas end_POSTSUBSCRIPT italic_ζ roman_ln ( divide start_ARG italic_R start_POSTSUBSCRIPT 0 end_POSTSUBSCRIPT end_ARG start_ARG italic_R ( italic_τ ) end_ARG ) . (25)

This restriction implies that we can solve for a relationship between the nongravitational acceleration magnitude |Ai|subscript𝐴𝑖|A_{i}|| italic_A start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | and the time τ𝜏\tauitalic_τ. With mX=18subscript𝑚𝑋18m_{X}=18italic_m start_POSTSUBSCRIPT italic_X end_POSTSUBSCRIPT = 18 amu (for H2O outgassing) and Tgas100similar-to-or-equalssubscript𝑇gas100T_{\rm gas}\simeq 100italic_T start_POSTSUBSCRIPT roman_gas end_POSTSUBSCRIPT ≃ 100 K, we find that vgas350similar-to-or-equalssubscript𝑣gas350v_{\rm gas}\simeq 350italic_v start_POSTSUBSCRIPT roman_gas end_POSTSUBSCRIPT ≃ 350 m s-1. Assuming that ζ=1𝜁1\zeta=1italic_ζ = 1, we can require that

|Ai|τ1 au/yr.subscript𝐴𝑖𝜏1 au/yr{|A_{i}|\,\tau\leq 1\text{ au/yr}}\,.| italic_A start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | italic_τ ≤ 1 au/yr . (26)

The typical (nonradiative) nongravitational acceleration of a dark comet is |Ai|105similar-to-or-equalssubscript𝐴𝑖superscript105|A_{i}|\simeq 10^{-5}| italic_A start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | ≃ 10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT au yr-2, which can therefore be maintained for τ105similar-to-or-equals𝜏superscript105\tau\simeq 10^{5}italic_τ ≃ 10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT yr.

It is worth noting that the stochastic outgassing does not perturb an orbit as significantly as a constant nongravitational acceleration in a single direction. The stochastic model is more representative of an ensemble of objects, allowing us to effectively explore the diffusion of these objects in orbital parameter space. The constant acceleration model would instead cause the population to follow well-defined paths in parameter space.

Appendix C Simulation Details

In this section, we discuss the details of the numerical simulations that we use to investigate the dark comets’ dynamical origins.

As discussed in Sec. 6, each simulation is constructed from sources in NEOMOD. We use 11 of the 12 sources available in NEOMOD — for convenience, the high-inclination Hungarias and Phocaeas are combined into a single high-inclination source. The other 10 sources — the JFCs, the forest of weak resonances in the inner main belt, and the ν6subscript𝜈6\nu_{6}italic_ν start_POSTSUBSCRIPT 6 end_POSTSUBSCRIPT, 3:1, 5:2, 7:3, 8:3, 9:4, 11:5, and 2:1 resonances are included directly.

From each source distribution, 104superscript10410^{4}10 start_POSTSUPERSCRIPT 4 end_POSTSUPERSCRIPT objects are generated as an initial population. While the NEOMOD sample determines the semimajor axis, eccentricity, and inclination, the other three orbital parameters (longitude of perihelion, longitude of the ascending node, and true anomaly) are chosen uniformly from [0,2π)02𝜋[0,2\pi)[ 0 , 2 italic_π ). This initial sample is then loaded into REBOUND (Rein & Liu, 2012), along with the Sun and all 8 solar system planets. Using REBOUNDx (Tamayo et al., 2020), each object is given a nongravitational acceleration of the form described in Appendix A and the typical (statistically significant, nonradiative) magnitude of the dark comets (105superscript10510^{-5}10 start_POSTSUPERSCRIPT - 5 end_POSTSUPERSCRIPT au yr-2; see Table 2).

While we do not include radiation pressure or the Yarkovsky effect in these simulations, these forces will not significantly affect our results. These forces are one to three orders of magnitude smaller than the dark comet nongravitational accelerations (see Sec. 2). Second, the larger radiation pressure force only operates in the radial direction and so affects the mean anomaly at epoch of an orbit. For sufficiently small nongravitational accelerations, the radiation pressure will not modify the semimajor axis, eccentricity, or inclination.

Each source population is then integrated for 105superscript10510^{5}10 start_POSTSUPERSCRIPT 5 end_POSTSUPERSCRIPT yr, which is consistent with the mass budget described in Appendix B. Our simulations use the MERCURIUS integrator (Rein et al., 2019), which uses the WHFast integrator for most circumstances and the IAS15 integrator for close approaches. The sample bodies are set to have no mass and do not interact with each other, although they will interact with the Sun and the planets. If an object collides with a planet, passes within 102superscript10210^{-2}10 start_POSTSUPERSCRIPT - 2 end_POSTSUPERSCRIPT au of the Sun, or passes >70absent70>70> 70 au from the Sun, it is removed from the simulation. The simulation timestep is set to a small value of 48484848 h. Every 100100100100 yr, the semimajor axis, eccentricity, and inclination of each object is recorded, providing a significant sample population of N107similar-to𝑁superscript107N\sim 10^{7}italic_N ∼ 10 start_POSTSUPERSCRIPT 7 end_POSTSUPERSCRIPT points. We show the initial and final samples for all of these sources in Figs. 818. In Fig. 19, we show the samples for the main belt as a whole, which was obtained by combining and normalizing all of the samples from the main belt resonances.

Table 4: NEO Source Probabilities. The probability that a randomly chosen NEO originates from a certain source population, according to the NEOMOD model. This is equivalent to the fraction of all NEOs that originate in that population.
Pop. Prob.
JFCs 0.0128
High i𝑖iitalic_i 0.0137
ν6subscript𝜈6\nu_{6}italic_ν start_POSTSUBSCRIPT 6 end_POSTSUBSCRIPT 0.6003
3:1 0.2979
5:2 0.3570
7:3 0.0012
8:3 0.0082
9:4 0.0008
11:5 0.0072
2:1 0.0032
Forest 0.0187

Appendix D Source Probability Calculations

In this section, we describe our methodology to calculate the source probabilities provided in Table 3.

We determine the probability that an object in a parameter-space bin Bisubscript𝐵𝑖B_{i}italic_B start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT belongs to an arbitrary source population Sjsubscript𝑆𝑗S_{j}italic_S start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT. In other words, we calculate the dependent probability P(Sj|Bi)𝑃conditionalsubscript𝑆𝑗subscript𝐵𝑖P(S_{j}|B_{i})italic_P ( italic_S start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT | italic_B start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ). Our simulations provide parameter-space locations for each individual source population (see Appendix C). We then calculate P(Bi|Sj)𝑃conditionalsubscript𝐵𝑖subscript𝑆𝑗P(B_{i}|S_{j})italic_P ( italic_B start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | italic_S start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ) by binning and normalizing these output samples. Using Bayes’ theorem, we write that

P(Sj|Bi)=P(Bi|Sj)P(Sj)P(Bi).𝑃conditionalsubscript𝑆𝑗subscript𝐵𝑖𝑃conditionalsubscript𝐵𝑖subscript𝑆𝑗𝑃subscript𝑆𝑗𝑃subscript𝐵𝑖{P(S_{j}|B_{i})=\frac{P(B_{i}|S_{j})\,P(S_{j})}{P(B_{i})}\,.}italic_P ( italic_S start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT | italic_B start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) = divide start_ARG italic_P ( italic_B start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | italic_S start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ) italic_P ( italic_S start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ) end_ARG start_ARG italic_P ( italic_B start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) end_ARG . (27)

The probability of finding a sample in a given bin P(Bi)𝑃subscript𝐵𝑖P(B_{i})italic_P ( italic_B start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) is the sum of the probabilities of being in that bin (given a source population) times the probability of being in that source population. That is,

P(Bi)=jP(Bi|Sj)P(Sj).𝑃subscript𝐵𝑖subscript𝑗𝑃conditionalsubscript𝐵𝑖subscript𝑆𝑗𝑃subscript𝑆𝑗{P(B_{i})=\sum_{j}P(B_{i}|S_{j})\,P(S_{j})\,.}italic_P ( italic_B start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) = ∑ start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT italic_P ( italic_B start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | italic_S start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ) italic_P ( italic_S start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ) . (28)

Note that P(Bi|Sj)𝑃conditionalsubscript𝐵𝑖subscript𝑆𝑗P(B_{i}|S_{j})italic_P ( italic_B start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT | italic_S start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ) is an output of our simulations.

The NEOMOD data file provides PNMD(Sj|Bi)subscript𝑃NMDconditionalsubscript𝑆𝑗subscript𝐵𝑖P_{\rm NMD}(S_{j}|B_{i})italic_P start_POSTSUBSCRIPT roman_NMD end_POSTSUBSCRIPT ( italic_S start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT | italic_B start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) and PNMD(Bi)subscript𝑃NMDsubscript𝐵𝑖P_{\rm NMD}(B_{i})italic_P start_POSTSUBSCRIPT roman_NMD end_POSTSUBSCRIPT ( italic_B start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ). Analogously to Eq. (28), we therefore find that

P(Sj)=iPNMD(Sj|Bi)PNMD(Bi).𝑃subscript𝑆𝑗subscript𝑖subscript𝑃NMDconditionalsubscript𝑆𝑗subscript𝐵𝑖subscript𝑃NMDsubscript𝐵𝑖P(S_{j})=\sum_{i}P_{\rm NMD}(S_{j}|B_{i})\,P_{\rm NMD}(B_{i})\,.italic_P ( italic_S start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ) = ∑ start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT italic_P start_POSTSUBSCRIPT roman_NMD end_POSTSUBSCRIPT ( italic_S start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT | italic_B start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) italic_P start_POSTSUBSCRIPT roman_NMD end_POSTSUBSCRIPT ( italic_B start_POSTSUBSCRIPT italic_i end_POSTSUBSCRIPT ) . (29)

These values of P(Sj)𝑃subscript𝑆𝑗P(S_{j})italic_P ( italic_S start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT ) are given in Table 4. Note that the P(Sj)𝑃subscript𝑆𝑗P(S_{j})italic_P ( italic_S start_POSTSUBSCRIPT italic_j end_POSTSUBSCRIPT )’s calculated in Eq. (29) are based on the NEOMOD model, which has different population distributions from our data (see Figs. 819).

We use Eq. (27) to calculate the probability that a given dark comet is from a given source. These values are reported in Table 3. The limits and bins on our calculations are identical to that of NEOMOD — 42 bins in semimajor axis from 0 to 4.2 au, 25 bins in eccentricity from 0 to 1, and 22 bins in inclination from 0 to 90.

Note that the values shown in Table 3 are only order-of-magnitude estimates as a result of several limitations to this methodology. First, we use the normalization from NEOMOD, rather than computing our own normalization by comparing to observational data. Second, in our calculations we marginalize over the absolute magnitude, reducing a dimension in our sample. Third, our probability calculations do not account for any morphological properties of the dark comets, which may further constraint the dark comets’ source population.

Refer to caption
Figure 8: The initial and final population distribution for the JFC source, using NEOMOD and REBOUNDx. The colorbar shows the probability of being found in the given projected bin.
Refer to caption
Figure 9: The initial and final population distribution for the high-inclination Hungarias and Phocaeas, using NEOMOD and REBOUNDx. The colorbar shows the probability of being found in the given projected bin.
Refer to caption
Figure 10: The initial and final population distribution for the weak resonance forest source, using NEOMOD and REBOUNDx. The colorbar shows the probability of being found in the given projected bin.
Refer to caption
Figure 11: The initial and final population distribution for the ν6subscript𝜈6\nu_{6}italic_ν start_POSTSUBSCRIPT 6 end_POSTSUBSCRIPT resonance source, using NEOMOD and REBOUNDx. The colorbar shows the probability of being found in the given projected bin.
Refer to caption
Figure 12: The initial and final population distribution for the 3:1 resonance source, using NEOMOD and REBOUNDx. The colorbar shows the probability of being found in the given projected bin.
Refer to caption
Figure 13: The initial and final population distribution for the 5:2 resonance source, using NEOMOD and REBOUNDx. The colorbar shows the probability of being found in the given projected bin.
Refer to caption
Figure 14: The initial and final population distribution for the 7:3 resonance source, using NEOMOD and REBOUNDx. The colorbar shows the probability of being found in the given projected bin.
Refer to caption
Figure 15: The initial and final population distribution for the 8:3 resonance source, using NEOMOD and REBOUNDx. The colorbar shows the probability of being found in the given projected bin.
Refer to caption
Figure 16: The initial and final population distribution for the 9:4 resonance source, using NEOMOD and REBOUNDx. The colorbar shows the probability of being found in the given projected bin.
Refer to caption
Figure 17: The initial and final population distribution for the 11:5 resonance source, using NEOMOD and REBOUNDx. The colorbar shows the probability of being found in the given projected bin.
Refer to caption
Figure 18: The initial and final population distribution for the 2:1 resonance source, using NEOMOD and REBOUNDx. The colorbar shows the probability of being found in the given projected bin.
Refer to caption
Figure 19: The initial and final population distribution for a sample drawn from all main belt objects, using NEOMOD and REBOUNDx. The colorbar shows the probability of being found in the given projected bin.