[orcid=0000-0002-0140-4475]

[1]

[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] [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

[cor1]Corresponding author

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

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

Aster G. Taylor<EMAIL_ADDRESS>[  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 $\nu_{6}$ 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) Celestial mechanics (211) Comets (280) Near-Earth Objects
(1092) Orbital 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 $\rho\sim 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 $\nu_{6}$ 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 $q\leq 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 $\rho\sim 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 $a\geq 5$ au and perihelia of $q\leq 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 $\sim 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 $\sim$1/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 $10-100$ m), (ii)
rapidly rotating ($0.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 $A_{i}\,(\text{1 au}/r)^{2}$
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/R_{\rm N}$ dependence, where $R_{\rm N}$ is the radius of the nucleus and
is computed from the absolute magnitude $H$ by

$R_{\rm N}=\bigg{(}\frac{1329\text{ km}}{2\sqrt{p}}\bigg{)}10^{-0.2H}\,.$ (1)

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

$A_{i,{\rm est}}=10^{0.2H+b_{i}}\text{ au/yr${}^{2}$}\,.$ (2)

In Eq. (2), $b_{i}$ is a constant scaling factor, which is $-10.734\pm 0.075$
for radiation pressure ($i=1$) and $-12.315\pm 0.014$ for the Yarkovsky effect
($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$ and the empirically-
determined $b_{i}$ 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$)
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 |  $A_{1}$a | Signif.b |  $A_{2}$c | Signif.d |  $A_{3}$e | Signif.f | max($C_{i}$)g
---|---|---|---|---|---|---|---
| [$10^{-5}$ au yr-2] | [$\sigma$] | [$10^{-5}$ au yr-2] | [$\sigma$] | [$10^{-5}$ au yr-2] | [$\sigma$] |
1998 KY26 | 2.31$\pm$ | 1.21 | 2 | -0.00168$\pm$ | 0.00081 | 2 | 0.426$\pm$ | 0.153 | 3 | 66.86
2005 VL1 | -8.87$\pm$ | 10.7 |  $<$1 | -0.00947$\pm$ | 0.00789 | 1 | -0.320$\pm$ | 0.0546 | 6 | 33.90
2016 NJ33 | 12.4$\pm$ | 3.94 | 3 | -0.00754$\pm$ | 0.00257 | 3 | 1.13$\pm$ | 0.217 | 5 | 186.38
2010 VL65 | 8.75$\pm$ | 17.3 |  $<$1 | -0.00195$\pm$ | 0.00711 |  $<$1 | -1.22$\pm$ | 0.173 | 7 | 36.01
2010 RF12 | 0.650$\pm$ | 0.795 |  $<$1 | -0.00181$\pm$ | 0.00381 |  $<$1 | -0.224$\pm$ | 0.028 | 8 | 9.58
2006 RH120 | 1.84$\pm$ | 0.107 | 18 | -0.676$\pm$ | 0.0849 | 8 | -0.173$\pm$ | 0.0426 | 4 | 17.58
2003 RM | -1.39$\pm$ | 1.62 |  $<$1 | 0.0286$\pm$ | 0.0005 | 56 | 0.0200$\pm$ | 0.0723 |  $<$1 | 679.90

  * a

Radial nongravitational acceleration.

  * b

Significance of $A_{1}$.

  * c

Transverse nongravitational acceleration.

  * d

Significance of $A_{2}$.

  * e

Out-of-plane nongravitational acceleration.

  * f

Significance of $A_{3}$.

  * g

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

Figure 1: Cumulative fraction/number of objects with nongravitational
accelerations smaller than a given cutoff, which is a factor of $C$ larger
than the estimated nongravitational acceleration from radiation pressure
($i=1$) or the Yarkovsky effect ($i=2$, and $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 $\sim$10 – 680$\times$ the expected radiative values, which are
shown as a shaded region bounded by vertical dashed lines.

We define a cutoff value $CA_{i,{\rm est}}$, 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 $R\simeq 100$ m and an albedo of $p=0.1$ is expected to have a
radiation-pressure acceleration of $A_{1,{\rm est}}\sim 4\times 10^{-7}$ au
yr-2. For a given $C$, we count the number of objects which have
nongravitational acceleration magnitudes smaller than this cutoff by
$>3\sigma$. Therefore, the estimated fraction of objects with nongravitational
accelerations definitively lower than $CA_{i,{\rm est}}$ is given by the
fraction of objects that satisfy

$|A_{i}|+3\sigma_{i}\leq CA_{i,{\rm est}}\,.$ (3)

In Eq. (3), $\sigma_{i}$ is the uncertainty in the $i$th 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$. The dark comet accelerations are
$\sim$10 – 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$ – $680$, 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 $2s$, 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$:$2s$:$2s$ (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$, $2s$, and $2s$ and
density $\rho$ in contact with one another such that their centers are
separated by a distance $s$. The physical binding force $F_{\rm s}$ of the
nucleus is the product of the contact area of the two cubes and their tensile
strength $\sigma_{t}$. The gravitational force is $Gm^{2}/s^{2}$ ($G$ is the
gravitational constant), the tensile strength force is $\sigma_{t}(2s)^{2}$,
and the mass of each prism is $m=4\rho s^{3}$. The binding force is therefore
given by

$F_{\rm bind}=16G\rho^{2}s^{4}+4\sigma_{t}s^{2}\,.$ (4)

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

$F_{\rm r}=2\rho s^{4}\omega^{2}.$ (5)

In Eq. (5), $\omega$ is the angular velocity, which is related to the rotation
period $P$ by $P=2\pi/\omega$. 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

$\sigma_{t}=\rho s^{2}\left(\frac{1}{2}\omega^{2}-4G\rho\right)\,.$ (6)

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

Figure 2: A diagram of our model of a cometary nucleus. We represent a
simplified nucleus as a cube with a side length of $2s$. The cube is composed
of two identical rectangular prisms with sides of length $s$, $2s$, $2s$. 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$a |  $e$b |  $i$c | Hd |  $R$e |  $P$f |  $\sigma_{t,c}>$g |  $\sigma_{t,a}>$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.41$-$1.99$ |  $0.0$-$1.1$ |  $0.0$-$5.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 $\rho\simeq 500$ kg/m3.

  * h

Minimum material strength for an asteroidal density of $\rho\simeq 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$ 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 $\sigma_{t,c}>$ and $\sigma_{t,a}>$ respectively.

These computed minimum strengths are on the order of $\sigma_{t,c}\gtrsim
0.01$ – $18$ Pa or $\sigma_{t,a}\gtrsim 0.0$ – $91$ 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
$\sigma_{t}\sim 1$–$10$ Pa (Sekanina & Yeomans, 1985; Asphaug & Benz, 1996;
Steckloff et al., 2015; Attree et al., 2018), with other estimates reaching as
high as $\sigma_{t}\sim 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$ 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 $\sim 2.2$ hours for asteroids (Pravec et al.,
2002; Warner et al., 2009; Sánchez & Scheeres, 2014) or $\sim 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

$\omega_{\rm crit}^{\rm str}=\sqrt{\frac{2\sigma_{t}}{\rho s^{2}}}\,.$ (7)

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

$\omega_{\rm crit}^{\rm grav}=\sqrt{\frac{4\pi\rho G}{3}}\,.$ (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 — $R\sim 200$ m for asteroids and $R\sim
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

$s_{\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}\,.$ (9)

Similarly, the transition size for comets is

$s_{\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}\,.$ (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 ($s\sim 600$ m; A’Hearn et al., 2011) possibly
representing a transition object.

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 ($\sim$1–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. 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 $\sim$20% 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 ($s\sim 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 $a\simeq 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 $\sim 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
$\gtrsim 500$ m in radius, based on Safrit et al., 2021).999Bottke et al.
(2023) found that the strength-gravity transition was at $R\sim 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 $\sigma_{t}\simeq 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 $\tau_{\rm dev}$ 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

$\tau_{\rm dev}=\frac{R^{2}}{\alpha}\,.$ (11)

In Eq. (11), $R$ is the radius of the dark comet nucleus and $\alpha$ is the
thermal diffusivity of the material, typically on the order of $\alpha\sim
10^{-8}-10^{-7}$ 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, $\tau_{\rm dev}\sim 1-1000$ years.

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

$\tau_{\rm tor}=\frac{R\sqrt{32\rho\sigma_{t}}}{P_{S}Y_{S}}\,.$ (12)

In Eq. (12), $P_{S}$ is the material sublimation pressure (Steckloff et al.,
2015; Steckloff & Jacobson, 2016) and $Y_{S}$ 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 $P_{S}$ for water ice is on the order of $P_{S}\sim
0.01$ – $0.1$ Pa (Steckloff & Jacobson, 2016). SYORP parameters for objects of
this size are thought to be on the order of $Y_{S}\sim 10^{-4}$ – $10^{-3}$
(Steckloff & Jacobson, 2016), but have been measured in larger, strength-
dominated comet nuclei to be an order of magnitude lower, at $Y_{S}\sim
10^{-5}$ – $10^{-4}$ (Steckloff et al., 2021). This results in SYORP
timescales (time required to monotonically spin up from rest to disruption) of
$\sim 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

$\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}$ (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.

Figure 5: The normalized size frequency distribution (SFD) of the objects
evolving under this model, for a range of removal rate parameters $q_{0}$. We
show $q_{0}\in[10^{-14},10^{-6}]$ s-1. We expect $q_{0}=10^{-14}$ s-1 to be
the true removal rate, and $q_{0}=10^{-6}$ s-1 is chosen because the
distribution is flat at this point. The model SFD $f(R)$ is shown in terms of
the steady-state divisionless SFD, $p(R)/q_{0}$. We assume that the injection
SFD power-law slope is $\beta=1$, but this only affects the exact values
rather than the slope. The predicted small-R limit of $f(R)\propto R\,p(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$
with two objects of radius $R/(2^{1/3})$. We ignore collisional fragmentation
in these calculations because the collisional timescale (Bottke et al., 2005)
is

$\tau_{\rm col}=1\text{ Myr }\left(\frac{R}{10\text{ m}}\right)^{1/2}\,.$ (14)

Since $\tau_{\rm div}\ll\tau_{\rm col}$ for the relevant sizes, collisions are
negligible. We then define the differential object fraction per radius $R$ as

$f(R)\equiv\frac{{\rm d}n}{{\rm d}R}\,.$ (15)

This is normalized such that the total number of objects
$N=\int_{0}^{\infty}\\!\\!f(R)\,{\rm d}R$.

Next, we define a function $p(R)$, which is the number of objects of a size
$R$ that are injected into the population over a time step ${\rm d}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 $q_{0}$,
which is the fraction of objects that are removed from the population in a
given time ${\rm d}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)$, which is the fraction of
objects of size $R$ that divide in a time step ${\rm d}t$. Combining these
definitions, we find that the change in the population is given by

$\frac{{\rm d}f}{{\rm d}t}(R)=p(R)+2q(\epsilon R)f(\epsilon
R)-q(R)f(R)-q_{0}f(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 $\epsilon 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, $\epsilon=2^{1/3}$. In order to find the
steady-state SFD, we set ${\rm d}f/{\rm d}t=0$.

We assume that $p(R)$ is constant in time, so that the age distribution of
objects is uniform. Therefore, if objects live in the NEO environment for
$10^{6}$ yr (Gladman et al., 2000; Nesvorný et al., 2023), then $q_{0}\simeq
10^{-14}$ s-1.

There are two competing factors to determine $q(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){\rm d}t/\tau_{\rm div}(R)$ objects are removed, where $\tau_{\rm div}$
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, $\tau_{\rm div}=\tau_{\rm tor}(1-1/\epsilon)$,
where $\tau_{\rm tor}$ is given by Eq. (12), and the factor of $(1-2^{-1/3})$
is introduced by beginning at $\omega=\omega_{\rm crit}(\epsilon R)$ rather
than $\omega=0$. With the devolatilization timescale given by Eq. (11), $q(R)$
is given by

$q(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}$ (17)

Solving Eq. (16) for $f(R)$ is complicated by the fact that both $f(R)$ and
$f(\epsilon R)$ appear in the equation. Moreover, a single solution does not
exist because $q(R)$ is a piecewise function of radius $R$. We therefore
consider three domains of the function.

Case 1. First we consider the case where $R$ is small and $q(R)=q(\epsilon
R)=0$. In this instance, Eq. (16) can be easily solved to show that

$f(R)=\frac{p(R)}{q_{0}}\,.$ (18)

In this limit, $f(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$ is large and $q(R)\neq 0$, at the
extreme end from the previous case. We can rewrite Eq. (16), solving for
$f(r)$ in terms of $f(\epsilon R)$ to obtain

$f(R)=\frac{p(R)+2q(\epsilon R)f(\epsilon R)}{q_{0}+q(R)}\,.$ (19)

We continue this series, writing $f(\epsilon R)$ in terms of
$f(\epsilon^{2}R)$ up to infinity. From Eq. (17), $q(R)=A/R$ where
$A=(1-2^{-1/3})^{-1}P_{S}Y_{S}/\sqrt{32\rho\sigma_{t}}$ is a constant.
Therefore, $q(\epsilon^{k}R)=\epsilon^{-k}q(R)$. We also assume that $p(R)$ is
a power law such that $p(R)\propto R^{-\beta}$. Continuing out the series and
simplifying, we find that

$\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}$
(20)

In general, Eq. (20) can be used to compute $f(R)$ for the domain where
$q(R)\neq 0$, so long as the series converges.

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

We show $f(R)$ for a range of $q_{0}$ values in Fig. 5. We normalize the SFD
to the small-$R$ limit $p(R)/q_{0}$. We assume that $\beta=1$, $P_{S}=0.1$ Pa,
$Y_{S}=10^{-4}$, $\sigma_{t}=50$ Pa, $\rho=500$ kg/m3, $\alpha=10^{-7}$ m2
s-1, and $\epsilon=2^{1/3}$, identically to Eq. (13). For small $R$, the form
of the SFD is unchanged from the input function (Eq. 18). For large $R$, the
relationship between Eq. (20) and $p(R)$ depends on the removal parameter
$q_{0}$. If $q_{0}=0$, then $f(R)\propto Rp(R)$ and the power-law index is
increased by a factor of order unity. Similarly, if $q_{0}\neq 0$ then
$f(R)\propto p(R)$. As $q_{0}\rightarrow 0$, the small-$n$ terms will be
dominated by $Rp(R)$, before eventually converging to $p(R)$ at large $R$.
Therefore, the value of $q_{0}$ 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 $q_{0}$.
The spikes at $R\sim 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 | $\nu_{6}$ | 3:1 | 5:2 | 7:3 | 8:3 | 9:4 | 11:5 | 2:1 | Forest | High $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
Figure 6: The bin probabilities for JFCs and main belt objects (from the
$\nu_{6}$ 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 $0^{\circ}$ and $15^{\circ}$ 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 $\nu_{6}$ 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 $\nu_{6}$ 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 $10^{4}$
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 $|A_{0}|=10^{-5}$ au
yr-2, consistent with the nongravitational accelerations of the dark comets
(see Table 2). Each source population is then numerically integrated for
$10^{5}$ 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 $\nu_{6}$ resonance,
and the 5:2 resonance (see Appendix C for more complete source region
calculations). We restrict the PDF to an inclination of $i\in[0,15]^{\circ}$,
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 $\nu_{6}$ 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.

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 $R\sim 1$ – $10$ 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 ($R\sim 10$ m) and have
densities of $\rho\simeq(500\pm 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 $\sim 10^{3}$ yr,
significantly shorter than the residence time in the near-Earth environment of
$\sim 10^{6}$ 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 $\sigma_{t}\geq 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, , 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, , 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, , 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, , 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, , 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, , 2, 14, doi: 10.3847/PSJ/abc9c8
  * Sánchez & Scheeres (2014) Sánchez, P., & Scheeres, D. J. 2014, , 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, , 2, 234, doi: 10.3847/PSJ/ac2dee
  * Seligman et al. (2023) Seligman, D. Z., Farnocchia, D., Micheli, M., et al. 2023, , 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

##  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

${\boldsymbol{F}=\alpha
A_{0}\left(\frac{r_{0}}{r}\right)^{2}\boldsymbol{\hat{e}}_{v}}\,.$ (21)

In Eq. (21), $\boldsymbol{\hat{e}}_{v}$ is the unit velocity vector, $r_{0}$
is a scaling distance (typically 1 au), $r$ is the heliocentric distance, and
$\alpha$ is a constant between -1 and 1. The parameter $\alpha$ 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$ is given by (Seligman et
al., 2023)

${\frac{{\rm d}M}{{\rm d}t}=\frac{M|A_{i}|}{v_{\rm gas}\zeta}}\,.$ (22)

The variable $\zeta$ indicates the collimation of the outgassing, and $v_{\rm
gas}$ is the gas velocity, given by

${v_{\rm gas}=\left(\frac{8k_{B}T_{\rm gas}}{\pi m_{X}}\right)^{1/2}}\,.$ (23)

Assuming that $A_{i}$ and $v_{\rm gas}$ are constants in time, Eq. (22) can be
solved to find that

${M(t)=M_{0}\exp\left(-\frac{|A_{i}|t}{v_{\rm gas}\zeta}\right)}\,.$ (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 $M\rightarrow 0$ only when
$t\rightarrow\infty$. However, we are here restricted by our (assumed)
knowledge of the initial and final states of the dark comets — objects of
$R_{\rm nuc}\sim 200$ m become objects of $R_{\rm nuc}\sim 10$ m. Assuming
that the dark comets are spheres, the acceleration and the timescale are
therefore constrained by

$|A_{i}|\tau\leq 3v_{\rm gas}\zeta\ln\left(\frac{R_{0}}{R(\tau)}\right)\,.$
(25)

This restriction implies that we can solve for a relationship between the
nongravitational acceleration magnitude $|A_{i}|$ and the time $\tau$. With
$m_{X}=18$ amu (for H2O outgassing) and $T_{\rm gas}\simeq 100$ K, we find
that $v_{\rm gas}\simeq 350$ m s-1. Assuming that $\zeta=1$, we can require
that

${|A_{i}|\,\tau\leq 1\text{ au/yr}}\,.$ (26)

The typical (nonradiative) nongravitational acceleration of a dark comet is
$|A_{i}|\simeq 10^{-5}$ au yr-2, which can therefore be maintained for
$\tau\simeq 10^{5}$ 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 $\nu_{6}$, 3:1, 5:2, 7:3, 8:3, 9:4, 11:5, and 2:1
resonances are included directly.

From each source distribution, $10^{4}$ 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\pi)$. 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
($10^{-5}$ 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 $10^{5}$ 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 $10^{-2}$ au of the Sun, or passes $>70$ au from the
Sun, it is removed from the simulation. The simulation timestep is set to a
small value of $48$ h. Every $100$ yr, the semimajor axis, eccentricity, and
inclination of each object is recorded, providing a significant sample
population of $N\sim 10^{7}$ points. We show the initial and final samples for
all of these sources in Figs. 8–18. 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$ | 0.0137
$\nu_{6}$ | 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 $B_{i}$
belongs to an arbitrary source population $S_{j}$. In other words, we
calculate the dependent probability $P(S_{j}|B_{i})$. Our simulations provide
parameter-space locations for each individual source population (see Appendix
C). We then calculate $P(B_{i}|S_{j})$ by binning and normalizing these output
samples. Using Bayes’ theorem, we write that

${P(S_{j}|B_{i})=\frac{P(B_{i}|S_{j})\,P(S_{j})}{P(B_{i})}\,.}$ (27)

The probability of finding a sample in a given bin $P(B_{i})$ 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(B_{i})=\sum_{j}P(B_{i}|S_{j})\,P(S_{j})\,.}$ (28)

Note that $P(B_{i}|S_{j})$ is an output of our simulations.

The NEOMOD data file provides $P_{\rm NMD}(S_{j}|B_{i})$ and $P_{\rm
NMD}(B_{i})$. Analogously to Eq. (28), we therefore find that

$P(S_{j})=\sum_{i}P_{\rm NMD}(S_{j}|B_{i})\,P_{\rm NMD}(B_{i})\,.$ (29)

These values of $P(S_{j})$ are given in Table 4. Note that the $P(S_{j})$’s
calculated in Eq. (29) are based on the NEOMOD model, which has different
population distributions from our data (see Figs. 8–19).

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.

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.  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.  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.  Figure 11:
The initial and final population distribution for the $\nu_{6}$ resonance
source, using NEOMOD and REBOUNDx. The colorbar shows the probability of being
found in the given projected bin.  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.  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.  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.  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.
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.  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.  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.
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.