# OSSOS. XXVI. On the Lack of Catastrophic Collisions in the Present Kuiper Belt

Abedin Abedin, J. Kavelaars, Jean-Marc Petit, Brett Gladman, Michele

Bannister, Mike Alexandersen, Ying-Tung Chen, Stephen Gwyn, Kathryn

Volk

► To cite this version:

Abedin Abedin, J. Kavelaars, Jean-Marc Petit, Brett Gladman, Michele Bannister, et al.. OSSOS. XXVI. On the Lack of Catastrophic Collisions in the Present Kuiper Belt. The Astronomical Journal, 2022, 164 (6), pp.261. 10.3847/1538-3881/ac9cdb. hal-04627165

HAL Id: hal-04627165

https://hal.science/hal-04627165

Submitted on 28 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

# OSSOS. XXVI. On the Lack of Catastrophic Collisions in the Present Kuiper Belt

Abedin Y. Abedin<fnref n="1" /> ID, J. J. Kavelaars<fnref n="1" /> ID, Jean-Marc Petit<fnref n="2" /> ID, Brett Gladman<fnref n="3" /> ID, Michele Bannister⁴ ID, Mike Alexandersen<fnref n="5" /> ID, Ying-Tung Chen<fnref n="6" />, Stephen Gwyn<fnref n="1" /> ID, and Kathryn Volk<fnref n="7" /> ID

*Received 2022 August 26; revised 2022 October 15; accepted 2022 October 20; published 2022 November 30*

## Abstract

We investigate different conditions, including the orbital and size-frequency distribution (SFD) of the early Kuiper Belt, that can trigger catastrophic planetesimal destruction. The goal of this study is to test if there is evidence for collisional grinding in the Kuiper Belt that has occurred since its formation. This analysis has important implications for whether the present-day SFD of the cold classical trans-Neptunian objects (TNOs) is a result of collisional equilibrium or if it reflects the primordial stage of planetesimal accretion. As an input to our modeling, we use the most up-to-date debiased OSSOS++ ensemble sample of the TNO population and orbital model based on the present-day architecture of the Kuiper Belt. We calculate the specific impact energies between impactor-target pairs from different TNO groups and compare our computed energies to catastrophic disruption results from smoothed particle hydrodynamics simulations. We explore different scenarios by considering different total primordial Kuiper Belt masses and power slopes of the SFD and allowing collisions to take place over different timescales. The collisional evolution of the Kuiper Belt is a strong function of the unknown initial mass in the trans-Neptunian region, where collisional grinding of planetesimals requires a total primordial Kuiper Belt mass of $M > 5 M_⊕$, collision speeds as high as $3 \text{ km s}^{-1}$, and collisions over at least 0.5 Gyr. We conclude that presently, most of the collisions in the trans-Neptunian region are in the cratering rather than disruption regime. Given the low collision rates among the cold classical Kuiper Belt objects, their SFD most likely represents the primordial planetesimal accretion.

*Unified Astronomy Thesaurus concepts: Trans-Neptunian objects (1705); Collisional processes (2286); Kuiper belt (893); Planetesimals (1259)*

## 1. Introduction

The size-frequency distribution (SFD) of planetesimals in the Kuiper Belt is a combined signature of the formation and collisional history of trans-Neptunian objects (TNOs), informing whether these objects are primordial or have undergone significant collisional erosion. Since the discovery of the first TNO (1992 QB₁; Jewitt & Luu 1993), the number of Kuiper Belt object (KBO) detections has increased markedly, owing to dedicated TNO surveys. These surveys, facilitated by TNO follow-up observations, allow orbital classification, which serves as a foundation for studying the dynamical and collisional history of that region of the solar system.

According to the classical Nice model (Gomes et al. 2005; Morbidelli 2005; Tsiganis et al. 2005) of the formation of the Kuiper Belt, the region of the solar system beyond 35 au was sparsely populated. Around 800–900 Myr after the giant planets formed, Jupiter and Saturn crossed the 1:2 mean-motion resonance (MMR), destabilizing the planetesimal disk and leading to major mass removal from the region between 15 and 34 au. A large fraction of planetesimals were scattering into the inner solar system, leading to the late heavy bombardment

(Gomes et al. 2005), while some were shoved outward into the Kuiper Belt. However, over the years, the “original” Nice model has undergone several modifications, especially concerning the starting giant planet orbital configuration and the timing of the dynamical instability; i.e., there are a number of pathways that may lead to the current configuration of the solar system. The main trigger of the dynamical instability of the giant planets has been attributed to planet–planet and planet–outer planetesimal disk interaction (Morbidelli 2005; Tsiganis et al. 2005; Levison et al. 2011; Nesvorný 2011). The outcome of this gravitational interaction is highly dependent on the initial heliocentric position of the planets and their separation and the distance between the outermost icy giant and the inner edge of the planetesimal disk. Recent studies find that the giant planet instability occurred on even shorter timescales, as originally proposed in the classical Nice model. Clement et al. (2018, 2019a, 2019b) argued that the giant planet instability may have occurred as early as the formation of the terrestrial planets, as opposed to ~100 Myr (Tsiganis et al. 2005), and the records of the late heavy bombardment on the Moon could be explained by accumulation of the remaining large planetesimals. De Sousa et al. (2020) found that the dynamical instability may actually have been triggered sometime between 4 and <60 Myr after the solar system’s formation, depending on whether the planets started in a stable or unstable orbital configuration.

<figure><img src="image_2.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ National Research Council of Canada, Herzberg Astronomy and Astrophysics, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada abedin.y.abedin@gmail.com</footnote>

<footnote n="2">² Institut UTINAM, CNRS-UMR 6213, Université Bourgogne Franche Comté BP 1615, F-25010 Besançon Cedex, France</footnote>

<footnote n="3">³ Department of Physics and Astronomy, 6224 Agricultural Road, University of British Columbia, Vancouver, BC, Canada
⁴ University of Canterbury, Christchurch, New Zealand</footnote>

<footnote n="5">⁵ Minor Planet Center, Smithsonian Astrophysical Observatory, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="6">⁶ Academia Sinica Institute of Astronomy and Astrophysics, Taiwan</footnote>

<footnote n="7">⁷ Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ, USA</footnote>
OPEN ACCESS

<figure><img src="image_2.png" /></figure>

CrossMark

# OSSOS. XXVI. On the Lack of Catastrophic Collisions in the Present Kuiper Belt

Abedin Y. Abedin<fnref n="1" /> ID, J. J. Kavelaars<fnref n="1" /> ID, Jean-Marc Petit<fnref n="2" /> ID, Brett Gladman<fnref n="3" /> ID, Michele Bannister⁴ ID, Mike Alexandersen<fnref n="5" /> ID, Ying-Tung Chen<fnref n="6" />, Stephen Gwyn<fnref n="1" /> ID, and Kathryn Volk<fnref n="7" /> ID

Received 2022 August 26; revised 2022 October 15; accepted 2022 October 20; published 2022 November 30

## Abstract

We investigate different conditions, including the orbital and size-frequency distribution (SFD) of the early Kuiper Belt, that can trigger catastrophic planetesimal destruction. The goal of this study is to test if there is evidence for collisional grinding in the Kuiper Belt that has occurred since its formation. This analysis has important implications for whether the present-day SFD of the cold classical trans-Neptunian objects (TNOs) is a result of collisional equilibrium or if it reflects the primordial stage of planetesimal accretion. As an input to our modeling, we use the most up-to-date debiased OSSOS++ ensemble sample of the TNO population and orbital model based on the present-day architecture of the Kuiper Belt. We calculate the specific impact energies between impactor-target pairs from different TNO groups and compare our computed energies to catastrophic disruption results from smoothed particle hydrodynamics simulations. We explore different scenarios by considering different total primordial Kuiper Belt masses and power slopes of the SFD and allowing collisions to take place over different timescales. The collisional evolution of the Kuiper Belt is a strong function of the unknown initial mass in the trans-Neptunian region, where collisional grinding of planetesimals requires a total primordial Kuiper Belt mass of $M > 5 M_⊕$, collision speeds as high as $3 \text{ km s}^{-1}$, and collisions over at least 0.5 Gyr. We conclude that presently, most of the collisions in the trans-Neptunian region are in the cratering rather than disruption regime. Given the low collision rates among the cold classical Kuiper Belt objects, their SFD most likely represents the primordial planetesimal accretion.

Unified Astronomy Thesaurus concepts: Trans-Neptunian objects (1705); Collisional processes (2286); Kuiper belt (893); Planetesimals (1259)

## 1. Introduction

The size-frequency distribution (SFD) of planetesimals in the Kuiper Belt is a combined signature of the formation and collisional history of trans-Neptunian objects (TNOs), informing whether these objects are primordial or have undergone significant collisional erosion. Since the discovery of the first TNO (1992 QB₁; Jewitt & Luu 1993), the number of Kuiper Belt object (KBO) detections has increased markedly, owing to dedicated TNO surveys. These surveys, facilitated by TNO follow-up observations, allow orbital classification, which serves as a foundation for studying the dynamical and collisional history of that region of the solar system.

According to the classical Nice model (Gomes et al. 2005; Morbidelli 2005; Tsiganis et al. 2005) of the formation of the Kuiper Belt, the region of the solar system beyond 35 au was sparsely populated. Around 800–900 Myr after the giant planets formed, Jupiter and Saturn crossed the 1:2 mean-motion resonance (MMR), destabilizing the planetesimal disk and leading to major mass removal from the region between 15 and 34 au. A large fraction of planetesimals were scattering into the inner solar system, leading to the late heavy bombardment

(Gomes et al. 2005), while some were shoved outward into the Kuiper Belt. However, over the years, the “original” Nice model has undergone several modifications, especially concerning the starting giant planet orbital configuration and the timing of the dynamical instability; i.e., there are a number of pathways that may lead to the current configuration of the solar system. The main trigger of the dynamical instability of the giant planets has been attributed to planet–planet and planet–outer planetesimal disk interaction (Morbidelli 2005; Tsiganis et al. 2005; Levison et al. 2011; Nesvorný 2011). The outcome of this gravitational interaction is highly dependent on the initial heliocentric position of the planets and their separation and the distance between the outermost icy giant and the inner edge of the planetesimal disk. Recent studies find that the giant planet instability occurred on even shorter timescales, as originally proposed in the classical Nice model. Clement et al. (2018, 2019a, 2019b) argued that the giant planet instability may have occurred as early as the formation of the terrestrial planets, as opposed to ~100 Myr (Tsiganis et al. 2005), and the records of the late heavy bombardment on the Moon could be explained by accumulation of the remaining large planetesimals. De Sousa et al. (2020) found that the dynamical instability may actually have been triggered sometime between 4 and <60 Myr after the solar system’s formation, depending on whether the planets started in a stable or unstable orbital configuration.

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote n="1">¹ National Research Council of Canada, Herzberg Astronomy and Astrophysics, 5071 West Saanich Road, Victoria, BC V9E 2E7, Canada abedin.y.abedin@gmail.com</footnote>

<footnote n="2">² Institut UTINAM, CNRS-UMR 6213, Université Bourgogne Franche Comté BP 1615, F-25010 Besançon Cedex, France</footnote>

<footnote n="3">³ Department of Physics and Astronomy, 6224 Agricultural Road, University of British Columbia, Vancouver, BC, Canada
⁴ University of Canterbury, Christchurch, New Zealand</footnote>

<footnote n="5">⁵ Minor Planet Center, Smithsonian Astrophysical Observatory, 60 Garden Street, Cambridge, MA 02138, USA</footnote>

<footnote n="6">⁶ Academia Sinica Institute of Astronomy and Astrophysics, Taiwan</footnote>

<footnote n="7">⁷ Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ, USA</footnote>
Attempts to determine the collisional history of the Kuiper Belt were made shortly after the detection of the first few dozen TNOs (e.g., Stern 1995, 1996; Farinella & Davis 1996; Davis & Farinella 1997). These initial works were based on calculating the collision rates within a disk of planetesimals beyond Neptune's orbit. The TNO SFD and total disk mass, however, were derived from discovery statistics of only ~20–30 objects in the size range 100 km $\lesssim r \lesssim$ 300 km. The abundance of smaller objects was obtained by direct extrapolation, significantly increasing the uncertainty in their number. These early investigations led to the conclusion that the objects with size $r \lesssim 10$ km were collisional fragments, whereas $r \gtrsim 100$ km class TNOs must be primordial. Unfortunately, most of these works were too plagued by low TNO number statistics to allow for more thorough Kuiper Belt orbital and SFD models.

Kenyon & Luu (1999) considered a more sophisticated model of planetesimal formation accounting for velocity evolution, as well as catastrophic and noncatastrophic interaction (merging, cratering, rebound, or disruption) between the colliding bodies. Their simulations indicated that, regardless of the initial conditions, the collisions led to a cumulative SFD described by a power law of the form $N \propto r^{-q}$ with two separate slope indices $q$. The authors found that the SFD of smaller TNOs ($r \lesssim 1$ km) is described by a power-law index of $q = 2.5$ and slightly steeper values ($q > 3$) for objects with $r \sim 1–1000$ km, concluding that fragmentation played a significant role in reshaping the early Kuiper Belt (see also Kenyon & Bromley 2004). Schlichting & Sari (2011) investigated protoplanet formation in the Kuiper Belt and concluded that the mass contained in large TNOs ($r > 50$ km) was primordial; i.e., these objects did not suffer catastrophic disruptions. Furthermore, their analysis indicated that the cumulative SFD for objects with $r \gtrsim 50$ km could be described with a single power law of $N(r) \propto r^{-4}$, though their simulations predicted much steeper SFD for objects with $r \lesssim 3$ km compared to New Horizons crater observations of the Pluto–Charon system (Singer et al. 2019) and the cold classical TNO Arrokoth (Spencer et al. 2020).

Other means of investigating the collisional evolution of solids focus on analytic and numerical treatment of the temporal evolution of the SFD of collisionally derived cascades (e.g., Kenyon & Bromley 2004, 2017, 2020; Pan & Sari 2005). Using a simple analytical approach but realistic collision scenarios (accounting for energy dissipation from the impact), Pan & Sari (2005) investigated the SFD of TNOs as a result of collisional fragmentation of ensembles of bodies. As with previous works, they fitted the TNO SFD with two different slopes. However, Pan & Sari (2005) found that the transition of the slope of the SFD is at $R \sim 40$ km, rather than $R = 70$ km (as derived from observations by Bernstein et al. 2004), where the slope of the fainter ($R < 40$ km) objects is described by a power-law index $q = 3$, which transitions to a steeper slope of $q = 5$ for larger TNOs ($R > 40$ km). That led the authors to conclude that this SFD and the location of the break are consistent with collisional evolution of the Kuiper Belt. However, deep sky surveys indicate multiple slope transitions and even different SFD for different TNO dynamical classes, which was not demonstrated by the analysis of Pan & Sari (2005). Recent analysis by Kenyon & Bromley (2020) of the temporal evolution of collisionally derived cascades indicates that collisions lead to a wavy equilibrium size distribution,

where the morphology of the waves is a function of the impact speed and the gravitational binding energy and bulk strength of solids. The authors found that an equilibrium wavy SFD matches well the New Horizons crater data on Charon (Singer et al. 2019) if TNO material ranges from weak to normal ice. Their simulations indicate that such an equilibrium SFD could be achieved within 100–300 Myr at a distance of $a = 45$ au for a protosolar nebula surface density profile $\Sigma(a) \approx 30 \text{ g cm}^{-2}(a/1 \text{ au})^{-3/2}$ and only 10–30 Myr at 25 au. Integrating this surface density profile function from $a = 40$ to 50 au yields a mass of the Kuiper Belt $M_{\text{KB}} \approx 0.1 M_{\oplus}$, almost an order of magnitude greater than the contemporary estimates of 0.01–0.02 $M_{\oplus}$ (e.g., Pitjeva & Pitjev 2018a, 2018b). One shortcoming of this model is the assumption of uniform collision speed, regardless of the orbits of the colliding bodies. A large fraction ($\gtrsim 30\%$) of the cold classical KBOs are loosely bound binaries (e.g., Noll et al. 2008; Grundy et al. 2011; Fraser et al. 2017, 2021), and their survival in a catastrophic collision environment is a question mark. Moreover, recent crater observations on the cold classical Kuiper Belt object (CCKBO) Arrokoth by the New Horizons mission (Spencer et al. 2020) provide direct evidence that collisions in that zone of the solar system were indeed sparse. Furthermore, the lack of large objects ($\gtrsim 400$ km) among the low-inclination component of the TNOs (Kavelaars et al. 2021) does not necessarily imply collisional grinding; rather, it appears to be consistent with a relatively new theory of planetesimal accretion via streaming instability (SI; e.g., Johansen et al. 2015; Simon et al. 2016; Schäfer et al. 2017; Rucska & Wadsley 2021). Independent studies show that the initial planetesimal mass function for objects smaller than ~100 km from SI simulations is well fit with a power law of the form $dN/dM \propto M^{-p}$, where $p \approx 1.6$, or $dN/dR \propto R^{-q}$, where $q \approx 2.5$. This results in a cumulative SFD $N(R) \propto R^{-q'}$ with ($q' \approx 3.5$), which is similar to the slope of 3.5 for collisionally derived cascades (Dohnanyi 1969). These indicate that the SFD of $R < 100$ km CCKBOs may not necessarily be a result of collisional grinding.

In this work, we aim to test if collisional destruction has been an important mechanism in shaping the SFD of TNOs or reflects the primordial planetesimal accretional stage. Furthermore, we explore which dynamical TNO population is more likely to experience catastrophic collisional destruction and initiate cascading fragmentation. To investigate this, we use the intrinsic collision probabilities (ICPs; defined in Wetherill 1967) calculated in Abedin et al. (2021; see also Greenstreet et al. 2015, 2016, 2019) and the most up-to-date debiased sample of dynamically different TNO populations from OSSOS++ (Kavelaars et al. 2009, 2020; Petit et al. 2011, 2021; Alexandersen et al. 2016; Bannister et al. 2016, 2018; Lawler et al. 2018b). We use that information to calculate the collision frequency and the expected radius $r$ of the largest impactor that has unit probability of colliding with a target of radius $R$ over some given timescale. These sets of calculations are performed for all different TNO populations. We next calculate the specific energies $Q_s$ involved in these collisions and compare them to the impact energy disruption thresholds $Q_D^*$ derived from smoothed particle hydrodynamic (SPH) simulations for icy bodies (e.g., Benz & Asphaug 1999; Leinhardt & Stewart 2009, 2012). That, altogether, sheds light on whether the impact energies in the present Kuiper Belt are high enough to initiate cascading fragmentation.
Attempts to determine the collisional history of the Kuiper Belt were made shortly after the detection of the first few dozen TNOs (e.g., Stern 1995, 1996; Farinella & Davis 1996; Davis & Farinella 1997). These initial works were based on calculating the collision rates within a disk of planetesimals beyond Neptune's orbit. The TNO SFD and total disk mass, however, were derived from discovery statistics of only ~20–30 objects in the size range 100 km $\lesssim r \lesssim$ 300 km. The abundance of smaller objects was obtained by direct extrapolation, significantly increasing the uncertainty in their number. These early investigations led to the conclusion that the objects with size $r \lesssim 10$ km were collisional fragments, whereas $r \gtrsim 100$ km class TNOs must be primordial. Unfortunately, most of these works were too plagued by low TNO number statistics to allow for more thorough Kuiper Belt orbital and SFD models.

Kenyon & Luu (1999) considered a more sophisticated model of planetesimal formation accounting for velocity evolution, as well as catastrophic and noncatastrophic interaction (merging, cratering, rebound, or disruption) between the colliding bodies. Their simulations indicated that, regardless of the initial conditions, the collisions led to a cumulative SFD described by a power law of the form $N \propto r^{-q}$ with two separate slope indices $q$. The authors found that the SFD of smaller TNOs ($r \lesssim 1$ km) is described by a power-law index of $q = 2.5$ and slightly steeper values ($q > 3$) for objects with $r \sim 1–1000$ km, concluding that fragmentation played a significant role in reshaping the early Kuiper Belt (see also Kenyon & Bromley 2004). Schlichting & Sari (2011) investigated protoplanet formation in the Kuiper Belt and concluded that the mass contained in large TNOs ($r > 50$ km) was primordial; i.e., these objects did not suffer catastrophic disruptions. Furthermore, their analysis indicated that the cumulative SFD for objects with $r \gtrsim 50$ km could be described with a single power law of $N(r) \propto r^{-4}$, though their simulations predicted much steeper SFD for objects with $r \lesssim 3$ km compared to New Horizons crater observations of the Pluto–Charon system (Singer et al. 2019) and the cold classical TNO Arrokoth (Spencer et al. 2020).

Other means of investigating the collisional evolution of solids focus on analytic and numerical treatment of the temporal evolution of the SFD of collisionally derived cascades (e.g., Kenyon & Bromley 2004, 2017, 2020; Pan & Sari 2005). Using a simple analytical approach but realistic collision scenarios (accounting for energy dissipation from the impact), Pan & Sari (2005) investigated the SFD of TNOs as a result of collisional fragmentation of ensembles of bodies. As with previous works, they fitted the TNO SFD with two different slopes. However, Pan & Sari (2005) found that the transition of the slope of the SFD is at $R \sim 40$ km, rather than $R = 70$ km (as derived from observations by Bernstein et al. 2004), where the slope of the fainter ($R < 40$ km) objects is described by a power-law index $q = 3$, which transitions to a steeper slope of $q = 5$ for larger TNOs ($R > 40$ km). That led the authors to conclude that this SFD and the location of the break are consistent with collisional evolution of the Kuiper Belt. However, deep sky surveys indicate multiple slope transitions and even different SFD for different TNO dynamical classes, which was not demonstrated by the analysis of Pan & Sari (2005). Recent analysis by Kenyon & Bromley (2020) of the temporal evolution of collisionally derived cascades indicates that collisions lead to a wavy equilibrium size distribution,

where the morphology of the waves is a function of the impact speed and the gravitational binding energy and bulk strength of solids. The authors found that an equilibrium wavy SFD matches well the New Horizons crater data on Charon (Singer et al. 2019) if TNO material ranges from weak to normal ice. Their simulations indicate that such an equilibrium SFD could be achieved within 100–300 Myr at a distance of $a = 45$ au for a protosolar nebula surface density profile $\Sigma(a) \approx 30 \text{ g cm}^{-2}(a/1\text{ au})^{-3/2}$ and only 10–30 Myr at 25 au. Integrating this surface density profile function from $a = 40$ to 50 au yields a mass of the Kuiper Belt $M_{\text{KB}} \approx 0.1 M_{\oplus}$, almost an order of magnitude greater than the contemporary estimates of 0.01–0.02 $M_{\oplus}$ (e.g., Pitjeva & Pitjev 2018a, 2018b). One shortcoming of this model is the assumption of uniform collision speed, regardless of the orbits of the colliding bodies. A large fraction ($\gtrsim 30\%$) of the cold classical KBOs are loosely bound binaries (e.g., Noll et al. 2008; Grundy et al. 2011; Fraser et al. 2017, 2021), and their survival in a catastrophic collision environment is a question mark. Moreover, recent crater observations on the cold classical Kuiper Belt object (CCKBO) Arrokoth by the New Horizons mission (Spencer et al. 2020) provide direct evidence that collisions in that zone of the solar system were indeed sparse. Furthermore, the lack of large objects ($\gtrsim 400$ km) among the low-inclination component of the TNOs (Kavelaars et al. 2021) does not necessarily imply collisional grinding; rather, it appears to be consistent with a relatively new theory of planetesimal accretion via streaming instability (SI; e.g., Johansen et al. 2015; Simon et al. 2016; Schäfer et al. 2017; Rucska & Wadsley 2021). Independent studies show that the initial planetesimal mass function for objects smaller than ~100 km from SI simulations is well fit with a power law of the form $dN/dM \propto M^{-p}$, where $p \approx 1.6$, or $dN/dR \propto R^{-q}$, where $q \approx 2.5$. This results in a cumulative SFD $N(R) \propto R^{-q'}$ with ($q' \approx 3.5$), which is similar to the slope of 3.5 for collisionally derived cascades (Dohnanyi 1969). These indicate that the SFD of $R < 100$ km CCKBOs may not necessarily be a result of collisional grinding.

In this work, we aim to test if collisional destruction has been an important mechanism in shaping the SFD of TNOs or reflects the primordial planetesimal accretional stage. Furthermore, we explore which dynamical TNO population is more likely to experience catastrophic collisional destruction and initiate cascading fragmentation. To investigate this, we use the intrinsic collision probabilities (ICPs; defined in Wetherill 1967) calculated in Abedin et al. (2021; see also Greenstreet et al. 2015, 2016, 2019) and the most up-to-date debiased sample of dynamically different TNO populations from OSSOS++ (Kavelaars et al. 2009, 2020; Petit et al. 2011, 2021; Alexandersen et al. 2016; Bannister et al. 2016, 2018; Lawler et al. 2018b). We use that information to calculate the collision frequency and the expected radius $r$ of the largest impactor that has unit probability of colliding with a target of radius $R$ over some given timescale. These sets of calculations are performed for all different TNO populations. We next calculate the specific energies $Q_s$ involved in these collisions and compare them to the impact energy disruption thresholds $Q_D^*$ derived from smoothed particle hydrodynamic (SPH) simulations for icy bodies (e.g., Benz & Asphaug 1999; Leinhardt & Stewart 2009, 2012). That, altogether, sheds light on whether the impact energies in the present Kuiper Belt are high enough to initiate cascading fragmentation.
In Section 2, we briefly revisit the orbital and SFD of TNOs (for full details, see, e.g., Bernstein et al. 2004; Kavelaars et al. 2009, 2021; Petit et al. 2011; Lawler et al. 2018a). In Section 3, we present our approach of calculating the collision frequencies in the Kuiper Belt and the associated impact energies, while in Section 4, we present our results and compare them to impact simulations. Finally, we conclude with a discussion of our key findings and summarize them in Section 5.

## 2. SFD and Orbital Model of TNOs

Presently, there are over 3500 reported TNOs with semimajor axis $a \gtrsim 30$ au in the Minor Planet Center database https://www.minorplanetcenter.net/iau/mpc.html, but only half of them have been tracked over multiple oppositions to allow orbital classification. Tracking of TNOs usually requires a significant amount of telescope time, as the mean motion of these objects (~0".1–0".3 day⁻¹) is so small that the sky motion is mostly due to Earth’s orbital motion. Therefore, the orbital distribution of TNOs is best characterized by dedicated sky surveys, which are usually allocated enough time to characterize their orbits. Orbital models combined with photometric observations allow for the calculation of the collision rates between dynamically different TNO populations.

In this work, we utilize the SFD of TNOs based on the debiased OSSOS++ ensemble sample (Kavelaars et al. 2009; Petit et al. 2011, 2017; Gladman et al. 2012; Alexandersen et al. 2016; Bannister et al. 2018; Lawler et al. 2018b). According to the nomenclature described in Gladman et al. (2008), the TNOs fall into several distinct dynamical groups.

1. Resonant objects—TNOs that are currently in MMR with Neptune.

2. Scattering objects—TNOs that experience close encounters with Neptune and consequently have an excursion in semimajor axis $|\Delta a| > 1.5$ au over 10 Myr of numerical orbit integration.

3. Classical or detached TNOs—Everything that does not fall into the previous two groups. The classical Kuiper Belt is further subdivided into the following groups, based on the orbital phase space the objects occupy.

(a) Inner classical belt—objects with semimajor axis interior to the 3:2 MMR with Neptune.

(b) Main classical belt—objects with semimajor axis between the 3:2 and 2:1 MMR with Neptune.

(c) Outer classical belt—TNOs with semimajor axis beyond the 2:1 MMR and orbital eccentricity $e < 0.24$.

(d) Detached objects—TNOs with semimajor axis beyond the 2:1 MMR with Neptune and orbital eccentricity $e > 0.24$.

For a thorough analysis of the orbital distribution of each TNO population, the reader is referred to, e.g., Petit et al. (2011). Throughout the text, we sometimes refer to dynamically cold and excited TNO populations. The former are objects with low inclination, $i < 5°$, whereas objects with $i > 5°$ are referred to as “excited” TNOs; i.e., their orbits have been affected significantly by Neptune. In essence, the resonant, hot classical, outer belt, detached, and scattering TNOs are considered as excited TNOs in this work.

Generally, solar system small-body observations yield the apparent magnitude *m* in a specific photometric band. Then, the cumulative apparent magnitude distribution of TNOs is given

by

$$N(<m) = \alpha \ln(10) 10^{\alpha(m-m_o)}, \quad (1)$$

where $N(<m)$ is the cumulative number of TNOs with apparent magnitude less than *m*, $m_o$ is a normalization coefficient, and $\alpha$ is referred to as the logarithmic slope of the distribution (see, e.g., Trujillo et al. 2001; Bernstein et al. 2004; Fraser & Kavelaars 2009; Fraser et al. 2014). Equation (1) can also be presented in terms of absolute magnitude *H* if the heliocentric distances of the objects are known,

$$N(<H) = \alpha \ln(10) 10^{\alpha(H-H_o)}, \quad (2)$$

where similarly, $N(<H)$ is the cumulative number of objects brighter than *H*, and $H_o$ is a normalization. Note that Equation (2) holds over some interval from *H* and $H+\Delta H$ for a particular slope $\alpha$ and normalization $H_o$; i.e., different absolute magnitude intervals may have different slopes and normalization coefficients. In fact, deep sky surveys indicate that the cumulative number of objects $\log(N)$ of TNOs as a function of the absolute magnitude *H* is well modeled with a series of power laws with different logarithmic slopes $\alpha$ (Bernstein et al. 2004; Kavelaars et al. 2009, 2021; Petit et al. 2021). Thus, the number of objects $N(<H)$ brighter than absolute magnitude *H* can be calculated as follows:

$$N(<H) = \sum_i \alpha_i \ln(10) 10^{\alpha_i(H-H_{oi})}, \quad (3)$$

where the summation is over each slope break, and $\alpha_i$ and $H_{oi}$ are the associated power slopes and normalization coefficients. Note that $H = H_{bi}$ for each slope break *i*; i.e., if the desired *H* is beyond the slope break, the equation is first solved for the portion until the break and then from the break to the desired *H*, hence the summation sign in Equation (3). Finally, the two portions are summed to obtain the cumulative $N(<H)$. If the albedo of a TNO in a given photometric band is known, the absolute magnitude can be converted to size via

$$R = C_x \frac{10^{-H_x/5}}{2\sqrt{p_x}} \text{ km,} \\ C_x \approx 2.99 \cdot 10^8 \times 10^{0.2X_{\odot}} \text{ km}, \quad (4)$$

where $R$, $H_x$, and $p_x$ are the radius, the absolute magnitude in a given photometric band *x*, and the albedo in the same photometric band (e.g., Masiero et al. 2021). The constant $C_x$ is photometric band-dependent and a function of the apparent stellar magnitude of the Sun in the *x* band, for instance, the constant $C_r \approx 1168$ for the *r* band of the Legacy Survey of the Canada-France-Hawaii Telescope, which gives the Sun’s apparent *r*-band magnitude as $r_{\odot} = -27.07$ (Gwyn 2012). To convert from absolute magnitude *H* to TNO radius, we assume a geometric albedo of 14% as deduced from the Herschel Space Observatory and Spitzer Space Telescope for the cold classical (Vilenius et al. 2014) and ~8% for the excited TNO populations (Vilenius et al. 2014; Kovalenko et al. 2017). For example, a cold classical TNO with $H_r = 3$ has a radius of roughly $R \approx 50$ km for an albedo of 14%, whereas a hot classical object of the same radius and albedo of 8% corresponds to $H_r \approx 2.5$. Thus, following Equation (3), the cumulative SFD of TNOs greater than some radius *R* is given
In Section 2, we briefly revisit the orbital and SFD of TNOs (for full details, see, e.g., Bernstein et al. 2004; Kavelaars et al. 2009, 2021; Petit et al. 2011; Lawler et al. 2018a). In Section 3, we present our approach of calculating the collision frequencies in the Kuiper Belt and the associated impact energies, while in Section 4, we present our results and compare them to impact simulations. Finally, we conclude with a discussion of our key findings and summarize them in Section 5.

## 2. SFD and Orbital Model of TNOs

Presently, there are over 3500 reported TNOs with semimajor axis $a \gtrsim 30$ au in the Minor Planet Center database https://www.minorplanetcenter.net/iau/mpc.html, but only half of them have been tracked over multiple oppositions to allow orbital classification. Tracking of TNOs usually requires a significant amount of telescope time, as the mean motion of these objects (~0".1–0".3 day⁻¹) is so small that the sky motion is mostly due to Earth’s orbital motion. Therefore, the orbital distribution of TNOs is best characterized by dedicated sky surveys, which are usually allocated enough time to characterize their orbits. Orbital models combined with photometric observations allow for the calculation of the collision rates between dynamically different TNO populations.

In this work, we utilize the SFD of TNOs based on the debiased OSSOS++ ensemble sample (Kavelaars et al. 2009; Petit et al. 2011, 2017; Gladman et al. 2012; Alexandersen et al. 2016; Bannister et al. 2018; Lawler et al. 2018b). According to the nomenclature described in Gladman et al. (2008), the TNOs fall into several distinct dynamical groups.

1. Resonant objects—TNOs that are currently in MMR with Neptune.

2. Scattering objects—TNOs that experience close encounters with Neptune and consequently have an excursion in semimajor axis $|\Delta a| > 1.5$ au over 10 Myr of numerical orbit integration.

3. Classical or detached TNOs—Everything that does not fall into the previous two groups. The classical Kuiper Belt is further subdivided into the following groups, based on the orbital phase space the objects occupy.

(a) Inner classical belt—objects with semimajor axis interior to the 3:2 MMR with Neptune.

(b) Main classical belt—objects with semimajor axis between the 3:2 and 2:1 MMR with Neptune.

(c) Outer classical belt—TNOs with semimajor axis beyond the 2:1 MMR and orbital eccentricity $e < 0.24$.

(d) Detached objects—TNOs with semimajor axis beyond the 2:1 MMR with Neptune and orbital eccentricity $e > 0.24$.

For a thorough analysis of the orbital distribution of each TNO population, the reader is referred to, e.g., Petit et al. (2011). Throughout the text, we sometimes refer to dynamically cold and excited TNO populations. The former are objects with low inclination, $i < 5°$, whereas objects with $i > 5°$ are referred to as “excited” TNOs; i.e., their orbits have been affected significantly by Neptune. In essence, the resonant, hot classical, outer belt, detached, and scattering TNOs are considered as excited TNOs in this work.

Generally, solar system small-body observations yield the apparent magnitude *m* in a specific photometric band. Then, the cumulative apparent magnitude distribution of TNOs is given

by

$$N(<m) = \alpha \ln(10) 10^{\alpha(m-m_o)}, \quad (1)$$

where $N(<m)$ is the cumulative number of TNOs with apparent magnitude less than *m*, $m_o$ is a normalization coefficient, and $\alpha$ is referred to as the logarithmic slope of the distribution (see, e.g., Trujillo et al. 2001; Bernstein et al. 2004; Fraser & Kavelaars 2009; Fraser et al. 2014). Equation (1) can also be presented in terms of absolute magnitude *H* if the heliocentric distances of the objects are known,

$$N(<H) = \alpha \ln(10) 10^{\alpha(H-H_o)}, \quad (2)$$

where similarly, $N(<H)$ is the cumulative number of objects brighter than *H*, and $H_o$ is a normalization. Note that Equation (2) holds over some interval from *H* and $H + \Delta H$ for a particular slope $\alpha$ and normalization $H_o$; i.e., different absolute magnitude intervals may have different slopes and normalization coefficients. In fact, deep sky surveys indicate that the cumulative number of objects $\log(N)$ of TNOs as a function of the absolute magnitude *H* is well modeled with a series of power laws with different logarithmic slopes $\alpha$ (Bernstein et al. 2004; Kavelaars et al. 2009, 2021; Petit et al. 2021). Thus, the number of objects $N(<H)$ brighter than absolute magnitude *H* can be calculated as follows:

$$N(<H) = \sum_i \alpha_i \ln(10) 10^{\alpha_i(H-H_{oi})}, \quad (3)$$

where the summation is over each slope break, and $\alpha_i$ and $H_{oi}$ are the associated power slopes and normalization coefficients. Note that $H = H_{bi}$ for each slope break $i$; i.e., if the desired $H$ is beyond the slope break, the equation is first solved for the portion until the break and then from the break to the desired $H$, hence the summation sign in Equation (3). Finally, the two portions are summed to obtain the cumulative $N(<H)$. If the albedo of a TNO in a given photometric band is known, the absolute magnitude can be converted to size via

$$R = C_x \frac{10^{-H_x/5}}{2\sqrt{p_x}} \text{ km,} \\ C_x \approx 2.99 \cdot 10^8 \times 10^{0.2X_{\odot}} \text{ km,} \quad (4)$$

where $R$, $H_x$, and $p_x$ are the radius, the absolute magnitude in a given photometric band $x$, and the albedo in the same photometric band (e.g., Masiero et al. 2021). The constant $C_x$ is photometric band-dependent and a function of the apparent stellar magnitude of the Sun in the $x$ band, for instance, the constant $C_r \approx 1168$ for the $r$ band of the Legacy Survey of the Canada-France-Hawaii Telescope, which gives the Sun’s apparent $r$-band magnitude as $r_{\odot} = -27.07$ (Gwyn 2012). To convert from absolute magnitude $H$ to TNO radius, we assume a geometric albedo of 14% as deduced from the Herschel Space Observatory and Spitzer Space Telescope for the cold classical (Vilenius et al. 2014) and ~8% for the excited TNO populations (Vilenius et al. 2014; Kovalenko et al. 2017). For example, a cold classical TNO with $H_r = 3$ has a radius of roughly $R \approx 50$ km for an albedo of 14%, whereas a hot classical object of the same radius and albedo of 8% corresponds to $H_r \approx 2.5$. Thus, following Equation (3), the cumulative SFD of TNOs greater than some radius *R* is given
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Schematic of the cumulative $H_r$ (r-band absolute magnitude distribution) of the excited (left) and cold classical (right) TNOs (Kavelaars et al. 2021; Petit et al. 2021). The colored labels denote the logarithmic slopes $\alpha$. The horizontal lines show the range of TNOs detected by surveys and inferred from crater mapping of the Pluto-Charon system and Arrokoth by the New Horizons mission science team (Singer et al. 2019; Spencer et al. 2020).</figcaption></figure>

by

$$N(>R) = \sum_{i} C_{i} R^{-q_{i}} = \sum_{i} \left( \frac{R}{R_{oi}} \right)^{-q_{i}}, \qquad (5)$$

where $q_i = 5\alpha_i + 1$ is referred to as the “power-law index” of the distribution, and $R_{oi}$ is a normalization associated with a particular slope. Thus, the cumulative number of TNOs $N(<H_r)$ or $N(>R)$ can be calculated for any $H_r$ using Equations (3) or (5), given that the cumulative population estimate is known at some absolute magnitude or radius ($H_r$ or $R$). The population estimates derived from the OSSOS+ + ensemble sample that we use to calculate the ICPs in Abedin et al. (2021) are scaled to absolute magnitude $H_r < 8.5$ and presented in Table 1.

Figure 1 depicts the cumulative absolute magnitude distribution of TNOs $N(<H_r)$ in the r band based on the debiased OSSOS++ ensemble sample for the cold and excited TNO populations (Kavelaars et al. 2021; Petit et al. 2021). Note that there are no cold classical objects brighter than $H_r = 3$, where the logarithmic slope dips to $\alpha \sim 10$ for $H_r \lesssim 5$. That implies that there are no dynamically cold TNOs with $R \gtrsim 200$ km.

The number of objects with $H_r \lesssim 10-11$ is derived from telescopic observations, while the SFD of fainter TNOs ($H_r > 10$) has been inferred from crater counts in the Pluto-Charon system (Singer et al. 2019) and the cold classical KBO Arrokoth (Spencer et al. 2020).

# 3. Modeling Collisions between TNOs

## 3.1. ICPs and Collision Rates

The present orbital structure of the Kuiper Belt (Kavelaars et al. 2009; Petit et al. 2011, 2017; Gladman et al. 2012; Alexandersen et al. 2016; Bannister et al. 2018; Lawler et al. 2018b) is such that collisions are rather infrequent. In this work, we assume that this orbital configuration has persisted for ~4 Gyr. Using the debiased TNO orbital model from OSSOS++, Abedin et al. (2021) calculated the ICPs within and between seven different dynamical TNO groups. Briefly, these groups are as follows.

<table><caption>**Table 1**
TNO Subcomponent Population Estimates</caption><thead><tr><th>Orbit Class</th><th>Fraction<br>% of Total</th><th>Population<br>N(H<sub>r</sub> &lt; 8.5)</th></tr></thead><tbody><tr><td>Cold main belt</td><td>8</td><td>16,400 ± 6000</td></tr><tr><td>Hot main belt</td><td>14</td><td>28,000 ± 9000</td></tr><tr><td>Detached/outer</td><td>25</td><td>50,000<sup>+35,000</sup><sub>-25,000</sub></td></tr><tr><td>Scattering</td><td>39</td><td>78,000 ± 17,000</td></tr><tr><td>Inner resonances</td><td></td><td></td></tr><tr><td>4:3</td><td>0.2</td><td>400<sup>+800</sup><sub>-120</sub></td></tr><tr><td>3:2</td><td>3.0</td><td>5700<sup>+2800</sup><sub>-3400</sub></td></tr><tr><td>Main resonances</td><td></td><td></td></tr><tr><td>5:3</td><td>2.6</td><td>5000<sup>+5200</sup><sub>-3000</sub></td></tr><tr><td>7:4</td><td>1.5</td><td>3000<sup>+4000</sup><sub>-2000</sub></td></tr><tr><td>Outer resonances</td><td></td><td></td></tr><tr><td>2:1</td><td>2.3</td><td>4400<sup>+900</sup><sub>-1200</sub></td></tr><tr><td>7:3</td><td>1.4</td><td>2500<sup>+4200</sup><sub>-2000</sub></td></tr><tr><td>5:2</td><td>2.8</td><td>5500<sup>+3500</sup><sub>-2500</sub></td></tr></tbody></table>

**References.** Gladman et al. (2012), Volk et al. (2016), Lawler et al. (2018b), Chen et al. (2019).

(i) Cold classical TNOs—objects with semimajor axis $40 \text{ au} < a < 47 \text{ au}$ and orbital inclination $i < 5^\circ$ (Petit et al. 2011).

(ii) Hot classical TNOs—objects with perihelion distance $35 \text{ au} < q < 40 \text{ au}$ and orbital inclination $i > 5^\circ$ (e.g., Petit et al. 2011).

We choose to subdivide the resonant TNO population into three distinct groups, which are in order of increasing semimajor axis as follows.

(iii) Inner resonant population—TNOs in the 4:3 and 3:2 MMR with Neptune.

(iv) Main resonant population—TNOs in the 5:3 and 7:4 MMR.

(v) Outer resonant population—TNOs in the 2:1, 7:3, and 5:2 MMR.

Splitting the resonant TNO population into smaller groups enables us to calculate the ICPs with a better resolution, i.e.,
<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Schematic of the cumulative $H_r$ (r-band absolute magnitude distribution) of the excited (left) and cold classical (right) TNOs (Kavelaars et al. 2021; Petit et al. 2021). The colored labels denote the logarithmic slopes $\alpha$. The horizontal lines show the range of TNOs detected by surveys and inferred from crater mapping of the Pluto–Charon system and Arrokoth by the New Horizons mission science team (Singer et al. 2019; Spencer et al. 2020).</figcaption></figure>

by

$$ N(>R) = \sum_{i} C_{i} R^{-q_{i}} \\ = \sum_{i} \left( \frac{R}{R_{oi}} \right)^{-q_{i}}, \qquad (5) $$

where $q_i = 5\alpha_i + 1$ is referred to as the “power-law index” of the distribution, and $R_{oi}$ is a normalization associated with a particular slope. Thus, the cumulative number of TNOs $N(<H_r)$ or $N(>R)$ can be calculated for any $H_r$ using Equations (3) or (5), given that the cumulative population estimate is known at some absolute magnitude or radius ($H_r$ or $R$). The population estimates derived from the OSSOS+ + ensemble sample that we use to calculate the ICPs in Abedin et al. (2021) are scaled to absolute magnitude $H_r < 8.5$ and presented in Table 1.

Figure 1 depicts the cumulative absolute magnitude distribution of TNOs $N(<H_r)$ in the r band based on the debiased OSSOS++ ensemble sample for the cold and excited TNO populations (Kavelaars et al. 2021; Petit et al. 2021). Note that there are no cold classical objects brighter than $H_r = 3$, where the logarithmic slope dips to $\alpha \sim 10$ for $H_r \lesssim 5$. That implies that there are no dynamically cold TNOs with $R \gtrsim 200$ km.

The number of objects with $H_r \lesssim 10–11$ is derived from telescopic observations, while the SFD of fainter TNOs ($H_r > 10$) has been inferred from crater counts in the Pluto–Charon system (Singer et al. 2019) and the cold classical KBO Arrokoth (Spencer et al. 2020).

# 3. Modeling Collisions between TNOs

## 3.1. ICPs and Collision Rates

The present orbital structure of the Kuiper Belt (Kavelaars et al. 2009; Petit et al. 2011, 2017; Gladman et al. 2012; Alexandersen et al. 2016; Bannister et al. 2018; Lawler et al. 2018b) is such that collisions are rather infrequent. In this work, we assume that this orbital configuration has persisted for ~4 Gyr. Using the debiased TNO orbital model from OSSOS++, Abedin et al. (2021) calculated the ICPs within and between seven different dynamical TNO groups. Briefly, these groups are as follows.

<table><caption>**Table 1**
TNO Subcomponent Population Estimates</caption><thead><tr><th>Orbit Class</th><th>Fraction<br>% of Total</th><th>Population<br>N(H<sub>r</sub> &lt; 8.5)</th></tr></thead><tbody><tr><td>Cold main belt</td><td>8</td><td>16,400 ± 6000</td></tr><tr><td>Hot main belt</td><td>14</td><td>28,000 ± 9000</td></tr><tr><td>Detached/outer</td><td>25</td><td>50,000<sup>+35,000</sup><sub>-25,000</sub></td></tr><tr><td>Scattering</td><td>39</td><td>78,000 ± 17,000</td></tr><tr><td>Inner resonances</td><td></td><td></td></tr><tr><td>4:3</td><td>0.2</td><td>400<sup>+800</sup><sub>-120</sub></td></tr><tr><td>3:2</td><td>3.0</td><td>5700<sup>+2800</sup><sub>-3400</sub></td></tr><tr><td>Main resonances</td><td></td><td></td></tr><tr><td>5:3</td><td>2.6</td><td>5000<sup>+5200</sup><sub>-3000</sub></td></tr><tr><td>7:4</td><td>1.5</td><td>3000<sup>+4000</sup><sub>-2000</sub></td></tr><tr><td>Outer resonances</td><td></td><td></td></tr><tr><td>2:1</td><td>2.3</td><td>4400<sup>+900</sup><sub>-1200</sub></td></tr><tr><td>7:3</td><td>1.4</td><td>2500<sup>+4200</sup><sub>-2000</sub></td></tr><tr><td>5:2</td><td>2.8</td><td>5500<sup>+3500</sup><sub>-2500</sub></td></tr></tbody></table>

**References.** Gladman et al. (2012), Volk et al. (2016), Lawler et al. (2018b), Chen et al. (2019).

(i) Cold classical TNOs—objects with semimajor axis $40 \text{ au} < a < 47 \text{ au}$ and orbital inclination $i < 5^\circ$ (Petit et al. 2011).

(ii) Hot classical TNOs—objects with perihelion distance $35 \text{ au} < q < 40 \text{ au}$ and orbital inclination $i > 5^\circ$ (e.g., Petit et al. 2011).

We choose to subdivide the resonant TNO population into three distinct groups, which are in order of increasing semimajor axis as follows.

(iii) Inner resonant population—TNOs in the 4:3 and 3:2 MMR with Neptune.

(iv) Main resonant population—TNOs in the 5:3 and 7:4 MMR.

(v) Outer resonant population—TNOs in the 2:1, 7:3, and 5:2 MMR.

Splitting the resonant TNO population into smaller groups enables us to calculate the ICPs with a better resolution, i.e.,
allowing us to discern which resonant subgroups are more likely to experience collisions.

Finally, we choose to combine the detached and outer belts into a single group, which we refer to as “detached,” and we separately consider the “scattering” population. We thus have

(vi) detached (detached + outer belts) and

(vii) scattering.

This splitting of the TNOs into somewhat arbitrary groups also allows us to pinpoint the heliocentric distance range where collisional erosion is more likely to occur. The ICPs are only a function of the orbital elements and do not involve the size of the colliding objects (e.g., Wetherill 1967). The calculated ICPs have units of collisions of km⁻² yr⁻¹. Thus, the collision rate Ṅ and number of collisions N between two objects are found from

$$ 
\begin{aligned}
\dot{N} &= P_i \cdot (R+r)^2 \text{ yr}^{-1} \\
\text{or} \\
N &= P_i \cdot \Delta T \cdot (R+r)^2,
\end{aligned}
\quad (6) 
$$

where Pᵢ is the ICP between a target of radius R and an impactor of radius r, and ΔT is the time interval. The term (R+r)² accounts for the cross section of the impactor and the target. We note that a π term in the cross section has been added to Pᵢ (see Abedin et al. 2021).

The collision rate Ṅ of a target of radius R with an ensemble of impactors with a size distribution dN/dr ∝ r⁻q is

$$ 
\begin{aligned}
\dot{N} &= P_i \int_{r_1}^{r_2} (R+r)^2 \left(\frac{dN}{dr}\right) dr \text{ yr}^{-1}, \\
\dot{N} &= P_i \int_{r_1}^{r_2} \int_{R_1}^{R_2} (R+r)^2 \left(\frac{dN}{dr}\right) dr \, dR \text{ yr}^{-1}.
\end{aligned}
\quad (7) 
$$

Keeping in mind the discussion above, Equation (7) reflects Ṅ over a given [r, r + Δr] or [R, R + ΔR] range for a particular power-law index q. To obtain the collision rate of a single target of radius R with an ensemble of impactors with a size spanning a range of different q, Equation (7) needs to be summed over different slope transitions of the SFD, i.e.,

$$ 
\begin{aligned}
\dot{N} &= P_i \sum_i \int_{r_{1i}}^{r_{2i}} (R+r)^2 \left(\frac{dN}{dr}\right)_i dr \text{ yr}^{-1} \\
\text{or} \\
\dot{N} &= P_i \sum_i C_i \int_{r_{1i}}^{r_{2i}} (R+r)^2 r^{-q_i} dr \text{ yr}^{-1},
\end{aligned}
\quad (8) 
$$

where each segment (slope) in the TNO SFD is associated with a given constant Cᵢ, a power-law index qᵢ, and integration limits corresponding to the break points r₁ᵢ and r₂ᵢ. Finally, when Equation (8) is multiplied by time ΔT, we obtain the number of collisions a target with radius R would experience with an ensemble of impactors with a given SFD over that time interval.

We use Equation (8) to calculate the largest impactor expected to hit a target of radius R over a given time frame. The

equation actually yields the number of impacts N(r) onto a given target of radius R, and we find the single largest impactor radius through interpolation to where N = 1. Thus, for each R, there is a unique value of r, and we obtain a function r(R), where r denotes the radius of the largest single impactor. However, there will be a different function r(R) for a given target depending on the dynamical TNO population it belongs to, the impactors it can collide with, and the time interval over which collisions are considered.

## 3.2. Specific Impact Energy, Scaling Laws, and Experiments

An important parameter that describes a collisional outcome is the specific impact energy Qs (the ratio of the projectile’s kinetic energy to the mass of the target),

$$ Q_s = \frac{E}{M} = \frac{1}{2} \frac{mv^2}{M} \text{ J kg}^{-1}, \quad (9) $$

where M and m are the target and impactor mass, and v is the impactor speed relative to the target.

The collision outcome can be classified into a few general categories depending on the specific impact energy (Benz & Asphaug 1999): merging, cratering, disruption, and dispersal of collisional fragments. In the merging and cratering regimes, the specific impact energy Qs must be smaller than the catastrophic disruption energy threshold Qs*. The latter is defined as the energy required to shatter the target and disperse half of the combined mass of the projectile and the target. The general form of this relationship as a function of impact speed and material properties has been derived from laboratory experiments (Housen & Holsapple 1990),

$$ Q_D^* \propto \left(\frac{\rho}{\delta}\right)^{(3\nu-1)} \nu^{(2-3\mu)} \times \left[ \left(\frac{S}{\rho}\right)^{\frac{-3\mu}{\nu-2}} R^{\frac{3\mu(\lambda+\tau)}{\nu-2}} + (\rho G)^{\frac{3\mu}{2}} R^{3\mu} \right], \quad (10) $$

where the first term in the brackets is the contribution of the strength regime, and the second term represents the gravity regime. The parameters ρ and δ are the target and impactor bulk densities, respectively; ν is the impact speed relative to the target; S is a measure of the material strength of the target; R is the radius of the target; and G is the gravitational constant. The dimensionless exponents λ and τ describe the size and strain-rate dependence of the material, whereas μ and ν are coupling parameters. For pure energy scaling of the impact event, μ = 2/3 and ν = 1/3, whereas for pure momentum scaling, μ = 1/3 and ν = 1/3. Typically, for most impact events, 1/3 < μ < 2/3, depending on the material porosity.

Modeling of planetary-scale collisions relies on SPH simulations, in which the target and impactor are represented by a large number of individual particles (Love & Ahrens 1996; Benz & Asphaug 1999; Jutzi et al. 2009; Leinhardt & Stewart 2009, 2012) held together by gravity. The dynamics of the postcollision remnants is then tracked by an N-body simulation package. Thus, laboratory impact experiments and SPH codes enable us to bridge the gap for collisions between centimeter- and meter-sized bodies to planetary scales (Holsapple & Schmidt 1987; Housen & Holsapple 1990; Housen et al. 1991, 2018; Holsapple 1993).
allowing us to discern which resonant subgroups are more likely to experience collisions.

Finally, we choose to combine the detached and outer belts into a single group, which we refer to as “detached,” and we separately consider the “scattering” population. We thus have

(vi) detached (detached + outer belts) and

(vii) scattering.

This splitting of the TNOs into somewhat arbitrary groups also allows us to pinpoint the heliocentric distance range where collisional erosion is more likely to occur. The ICPs are only a function of the orbital elements and do not involve the size of the colliding objects (e.g., Wetherill 1967). The calculated ICPs have units of collisions of km⁻² yr⁻¹. Thus, the collision rate Ṅ and number of collisions N between two objects are found from

$$ \begin{aligned} \dot{N} &= P_i \cdot (R+r)^2 \text{ yr}^{-1} \\ \text{or} \\ N &= P_i \cdot \Delta T \cdot (R+r)^2, \end{aligned} \qquad (6) $$

where Pᵢ is the ICP between a target of radius R and an impactor of radius r, and ΔT is the time interval. The term (R + r)² accounts for the cross section of the impactor and the target. We note that a π term in the cross section has been added to Pᵢ (see Abedin et al. 2021).

The collision rate Ṅ of a target of radius R with an ensemble of impactors with a size distribution dN/dr ∝ r⁻q is

$$ \begin{aligned} \dot{N} &= P_i \int_{r_1}^{r_2} (R+r)^2 \left(\frac{dN}{dr}\right) dr \text{ yr}^{-1}, \\ \dot{N} &= P_i \int_{r_1}^{r_2} \int_{R_1}^{R_2} (R+r)^2 \left(\frac{dN}{dr}\right) dr \, dR \text{ yr}^{-1}. \end{aligned} \qquad (7) $$

Keeping in mind the discussion above, Equation (7) reflects Ṅ over a given [r, r + Δr] or [R, R + ΔR] range for a particular power-law index q. To obtain the collision rate of a single target of radius R with an ensemble of impactors with a size spanning a range of different q, Equation (7) needs to be summed over different slope transitions of the SFD, i.e.,

$$ \begin{aligned} \dot{N} &= P_i \sum_i \int_{r_{1i}}^{r_{2i}} (R+r)^2 \left(\frac{dN}{dr}\right)_i dr \text{ yr}^{-1} \\ \text{or} \\ \dot{N} &= P_i \sum_i C_i \int_{r_{1i}}^{r_{2i}} (R+r)^2 r^{-q_i} dr \text{ yr}^{-1}, \end{aligned} \qquad (8) $$

where each segment (slope) in the TNO SFD is associated with a given constant Cᵢ, a power-law index qᵢ, and integration limits corresponding to the break points r₁ᵢ and r₂ᵢ. Finally, when Equation (8) is multiplied by time ΔT, we obtain the number of collisions a target with radius R would experience with an ensemble of impactors with a given SFD over that time interval.

We use Equation (8) to calculate the largest impactor expected to hit a target of radius R over a given time frame. The

equation actually yields the number of impacts N(r) onto a given target of radius R, and we find the single largest impactor radius through interpolation to where N = 1. Thus, for each R, there is a unique value of r, and we obtain a function r(R), where r denotes the radius of the largest single impactor. However, there will be a different function r(R) for a given target depending on the dynamical TNO population it belongs to, the impactors it can collide with, and the time interval over which collisions are considered.

## 3.2. Specific Impact Energy, Scaling Laws, and Experiments

An important parameter that describes a collisional outcome is the specific impact energy Qs (the ratio of the projectile’s kinetic energy to the mass of the target),

$$ Q_s = \frac{E}{M} = \frac{1}{2} \frac{mv^2}{M} \text{ J kg}^{-1}, \qquad (9) $$

where M and m are the target and impactor mass, and v is the impactor speed relative to the target.

The collision outcome can be classified into a few general categories depending on the specific impact energy (Benz & Asphaug 1999): merging, cratering, disruption, and dispersal of collisional fragments. In the merging and cratering regimes, the specific impact energy Qs must be smaller than the catastrophic disruption energy threshold Qs*. The latter is defined as the energy required to shatter the target and disperse half of the combined mass of the projectile and the target. The general form of this relationship as a function of impact speed and material properties has been derived from laboratory experiments (Housen & Holsapple 1990),

$$ Q_D^* \propto \left(\frac{\rho}{\delta}\right)^{(3\nu-1)} \nu^{(2-3\mu)} \times \left[ \left(\frac{S}{\rho}\right)^{\frac{-3\mu}{\nu-2}} R^{\frac{3\mu(\lambda+\tau)}{\nu-2}} + (\rho G)^{\frac{3\mu}{2}} R^{3\mu} \right], \qquad (10) $$

where the first term in the brackets is the contribution of the strength regime, and the second term represents the gravity regime. The parameters ρ and δ are the target and impactor bulk densities, respectively; ν is the impact speed relative to the target; S is a measure of the material strength of the target; R is the radius of the target; and G is the gravitational constant. The dimensionless exponents λ and τ describe the size and strain-rate dependence of the material, whereas μ and ν are coupling parameters. For pure energy scaling of the impact event, μ = 2/3 and ν = 1/3, whereas for pure momentum scaling, μ = 1/3 and ν = 1/3. Typically, for most impact events, 1/3 < μ < 2/3, depending on the material porosity.

Modeling of planetary-scale collisions relies on SPH simulations, in which the target and impactor are represented by a large number of individual particles (Love & Ahrens 1996; Benz & Asphaug 1999; Jutzi et al. 2009; Leinhardt & Stewart 2009, 2012) held together by gravity. The dynamics of the postcollision remnants is then tracked by an N-body simulation package. Thus, laboratory impact experiments and SPH codes enable us to bridge the gap for collisions between centimeter- and meter-sized bodies to planetary scales (Holsapple & Schmidt 1987; Housen & Holsapple 1990; Housen et al. 1991, 2018; Holsapple 1993).
<table><caption>**Table 2**
Values of ($Q_o$, $B$, $a$, $b$) Parameters for Different Impact Speeds for an Icy Target for Head-on Collisions</caption><thead><tr><th>Material</th><th><i>v</i><sub>impact</sub><br>(km s<sup>-1</sup>)</th><th><i>Q</i><sub>o</sub><br>(J kg<sup>-1</sup>)</th><th><i>B</i><br>(J m<sup>3</sup> kg<sup>-2</sup>)</th><th><i>a</i></th><th><i>b</i></th></tr></thead><tbody><tr><td>Ice</td><td>0.5</td><td>7.0 &times; 10<sup>3</sup></td><td>2.1 &times; 10<sup>-7</sup></td><td>-0.45</td><td>1.19</td></tr><tr><td>Ice</td><td>3.0</td><td>1.6 &times; 10<sup>3</sup></td><td>1.2 &times; 10<sup>-7</sup></td><td>-0.39</td><td>1.26</td></tr></tbody></table>

Note. The values are taken from Table 3 of Benz & Asphaug (1999).

Benz & Asphaug (1999) used an SPH method to calculate the disruption energy threshold for head-on collisions over 8 orders of magnitudes in target size (centimeters to hundreds of kilometers). Their model included dynamical fracture of brittle material, proving to match laboratory experiments exceptionally well. Benz & Asphaug (1999) found that their model disruption thresholds are well matched by the form

$$Q_D^* = Q_o \left(\frac{R}{1 \text{ m}}\right)^a + B \rho \left(\frac{R}{1 \text{ m}}\right)^b, \quad (11)$$

where the first term on the right describes the strength regime, and the second one is the gravity regime. The quantities ($Q_o$, $B$, $a$, $b$) are fit parameters to the modeling. For convenience, we reproduce the best-fit values in Table 2 for different impact speeds.

Figures 6 and 7 in Benz & Asphaug (1999) illustrate the dominance of the material strength for objects with $R \lesssim 150$ m. For $R \gtrsim 150$ m, the self-gravity of the target becomes important and dictates the collisional outcome. However, most of the impacts between TNOs are unlikely to be head-on; rather, an oblique angle is more probable. Leinhardt & Stewart (2012) expanded on the work by Benz & Asphaug (1999) by deriving equations for the catastrophic disruption threshold $Q_D^*$ for gravity-dominated bodies in oblique impacts. Considering oblique impacts compared to a head-on collision increases the disruption energy threshold $Q_D^*$. The reason is that only a fraction of the impactor mass will interact with the target, and hence less kinetic energy is deposited into the target (Leinhardt & Stewart 2012).

Using the radius $r$ of the largest impactor, derived in Section 3.1, impacting a target of radius $R$ from a given dynamical TNO class, we calculate the specific impact energy using Equation (9). Throughout this work, we assume an equal impactor and target density of $\rho = 500 \text{ kg m}^{-3}$ and a uniform impact speed of $v = 3 \text{ km s}^{-1}$. The choice of $v = 3 \text{ km s}^{-1}$ is motivated by the fact that most of the collisions in the trans-Neptunian region are likely to occur at lower impact speeds (Greenstreet et al. 2015, 2016; Abedin et al. 2021). Hence, the impact speed of $v = 3 \text{ km s}^{-1}$ will serve as an upper limit of $Q_s$. Thus, for each target from a given TNO population, we calculate a collection of functions $Q_s = Q_s(R)$ for each impactor population and compare them to SPH simulations for icy bodies by Benz & Asphaug (1999) and Leinhardt & Stewart (2012).

### 3.3. Exploring Different Cases of the Kuiper Belt Structure

Current estimates place the mass of the Kuiper Belt in the range $M_{KB} \approx 0.01-0.02 M_\oplus$ (Pitjeva & Pitjev 2018a, 2018b), though during the major planet migration, the mass of the primordial belt was likely higher. The Nice model (Gomes et al. 2005; Morbidelli 2005; Tsiganis et al. 2005) suggests that

the total mass of planetesimals between 15 and 35 au (before the dynamical instability of the giant planets' orbits) was $\sim 35 M_\oplus$. More recent simulations indicate a total planetesimal mass between 20 and 30 au (beyond Neptune's initial orbit) of $15-20 M_\oplus$ (Nesvorný & Vokrouhlický 2016; Nesvorný 2018). Most of that mass was scattered into the inner solar system, while a small fraction of the remaining mass formed the hot component of the Kuiper Belt. The exact amount of initial mass scattered into the trans-Neptunian region is uncertain and is the key to understanding if collisions have played a major role in reshaping the Kuiper Belt.

The cases we consider are presented in Table 3. Note that all cases presented below are based on the ICPs, calculated assuming the current orbital architecture of the Kuiper Belt. The different cases can be summarized as follows.

1. Nominal case—Under this case, the orbital distribution and SFD of TNOs is the present one, as derived from OSSOS++, and collisions are considered over a 4 Gyr timescale. This scenario explores if TNOs have suffered significant collisional fragmentation since the giant planets attained their present orbits approximately 4 Gyr ago. The population estimates for different TNO dynamical groups are those presented in Abedin et al. (2021).

2. Case_x200_4GYR—This case assumes a planetesimal disk mass 200 times the current mass of the Kuiper Belt and collisions over 4 Gyr. This case explores the required planetesimal mass that leads to collisional grinding over the history of the solar system. Of course, there is no observational indication of such high present-day Kuiper Belt mass, though this case is meant to demonstrate some of the conditions needed for catastrophic disruption.

3. Case_x500_300MYR—The giant planet instability was likely followed by removal of a large fraction of the planetesimal mass from the region of the icy giants. However, the rate of mass depletion is uncertain, motivating us to explore an initial planetesimal mass equivalent to $\sim 500$ times (or $\sim 5 M_\oplus$) the present Kuiper Belt mass and a depletion timescale of 300 Myr. This timescale is longer than might be expected for the instability period, but this longer timescale partially compensates for the lower collision rates we compute based on the current structure of the Kuiper Belt as compared to the more compact form in the Nice model. This scenario could shed light on the efficiency of mass removal versus the collisional grinding of TNOs. That is, if the mass depletion rate was higher than the catastrophic collision lifetime of planetesimals, then most remaining TNOs must be collisionally pristine.

4. Case_x2000_300MYR—This is the same as in case_x500_300MYR, except the initial planetesimal mass immediately following the giant planet instability was $\sim 10-20 M_\oplus$.

5. Case_H85_07_300MYR—This case explores the collisional evolution (catastrophic versus cratering) of planetesimals if the slope of the SFD above the knee ($H_r > 8.5$) was steeper ($\alpha = 0.7$) compared to the present-day value of $\alpha = 0.45$. Note that the SFD of planetesimals brighter than the knee is assumed to be the present one. The steeper power-law slope above the knee effectively increases the number of objects with radii $r \lesssim 50$ km, which, on the other hand, increases the total collisional
<table><caption>**Table 2**
Values of ($Q_o$, $B$, $a$, $b$) Parameters for Different Impact Speeds for an Icy Target for Head-on Collisions</caption><thead><tr><th>Material</th><th><i>v</i><sub>impact</sub><br>(km s<sup>-1</sup>)</th><th><i>Q</i><sub>o</sub><br>(J kg<sup>-1</sup>)</th><th><i>B</i><br>(J m<sup>3</sup> kg<sup>-2</sup>)</th><th><i>a</i></th><th><i>b</i></th></tr></thead><tbody><tr><td>Ice</td><td>0.5</td><td>7.0 &times; 10<sup>3</sup></td><td>2.1 &times; 10<sup>-7</sup></td><td>-0.45</td><td>1.19</td></tr><tr><td>Ice</td><td>3.0</td><td>1.6 &times; 10<sup>3</sup></td><td>1.2 &times; 10<sup>-7</sup></td><td>-0.39</td><td>1.26</td></tr></tbody></table>

Note. The values are taken from Table 3 of Benz & Asphaug (1999).

Benz & Asphaug (1999) used an SPH method to calculate the disruption energy threshold for head-on collisions over 8 orders of magnitudes in target size (centimeters to hundreds of kilometers). Their model included dynamical fracture of brittle material, proving to match laboratory experiments exceptionally well. Benz & Asphaug (1999) found that their model disruption thresholds are well matched by the form

$$Q_D^* = Q_o \left(\frac{R}{1 \text{ m}}\right)^a + B \rho \left(\frac{R}{1 \text{ m}}\right)^b, \quad (11)$$

where the first term on the right describes the strength regime, and the second one is the gravity regime. The quantities ($Q_o$, $B$, $a$, $b$) are fit parameters to the modeling. For convenience, we reproduce the best-fit values in Table 2 for different impact speeds.

Figures 6 and 7 in Benz & Asphaug (1999) illustrate the dominance of the material strength for objects with $R \lesssim 150$ m. For $R \gtrsim 150$ m, the self-gravity of the target becomes important and dictates the collisional outcome. However, most of the impacts between TNOs are unlikely to be head-on; rather, an oblique angle is more probable. Leinhardt & Stewart (2012) expanded on the work by Benz & Asphaug (1999) by deriving equations for the catastrophic disruption threshold $Q_D^*$ for gravity-dominated bodies in oblique impacts. Considering oblique impacts compared to a head-on collision increases the disruption energy threshold $Q_D^*$. The reason is that only a fraction of the impactor mass will interact with the target, and hence less kinetic energy is deposited into the target (Leinhardt & Stewart 2012).

Using the radius $r$ of the largest impactor, derived in Section 3.1, impacting a target of radius $R$ from a given dynamical TNO class, we calculate the specific impact energy using Equation (9). Throughout this work, we assume an equal impactor and target density of $\rho = 500 \text{ kg m}^{-3}$ and a uniform impact speed of $v = 3 \text{ km s}^{-1}$. The choice of $v = 3 \text{ km s}^{-1}$ is motivated by the fact that most of the collisions in the trans-Neptunian region are likely to occur at lower impact speeds (Greenstreet et al. 2015, 2016; Abedin et al. 2021). Hence, the impact speed of $v = 3 \text{ km s}^{-1}$ will serve as an upper limit of $Q_s$. Thus, for each target from a given TNO population, we calculate a collection of functions $Q_s = Q_s(R)$ for each impactor population and compare them to SPH simulations for icy bodies by Benz & Asphaug (1999) and Leinhardt & Stewart (2012).

### 3.3. Exploring Different Cases of the Kuiper Belt Structure

Current estimates place the mass of the Kuiper Belt in the range $M_{KB} \approx 0.01-0.02 M_\oplus$ (Pitjeva & Pitjev 2018a, 2018b), though during the major planet migration, the mass of the primordial belt was likely higher. The Nice model (Gomes et al. 2005; Morbidelli 2005; Tsiganis et al. 2005) suggests that

the total mass of planetesimals between 15 and 35 au (before the dynamical instability of the giant planets' orbits) was $\sim 35 M_\oplus$. More recent simulations indicate a total planetesimal mass between 20 and 30 au (beyond Neptune's initial orbit) of $15-20 M_\oplus$ (Nesvorný & Vokrouhlický 2016; Nesvorný 2018). Most of that mass was scattered into the inner solar system, while a small fraction of the remaining mass formed the hot component of the Kuiper Belt. The exact amount of initial mass scattered into the trans-Neptunian region is uncertain and is the key to understanding if collisions have played a major role in reshaping the Kuiper Belt.

The cases we consider are presented in Table 3. Note that all cases presented below are based on the ICPs, calculated assuming the current orbital architecture of the Kuiper Belt. The different cases can be summarized as follows.

1. Nominal case—Under this case, the orbital distribution and SFD of TNOs is the present one, as derived from OSSOS++, and collisions are considered over a 4 Gyr timescale. This scenario explores if TNOs have suffered significant collisional fragmentation since the giant planets attained their present orbits approximately 4 Gyr ago. The population estimates for different TNO dynamical groups are those presented in Abedin et al. (2021).

2. Case_x200_4GYR—This case assumes a planetesimal disk mass 200 times the current mass of the Kuiper Belt and collisions over 4 Gyr. This case explores the required planetesimal mass that leads to collisional grinding over the history of the solar system. Of course, there is no observational indication of such high present-day Kuiper Belt mass, though this case is meant to demonstrate some of the conditions needed for catastrophic disruption.

3. Case_x500_300MYR—The giant planet instability was likely followed by removal of a large fraction of the planetesimal mass from the region of the icy giants. However, the rate of mass depletion is uncertain, motivating us to explore an initial planetesimal mass equivalent to $\sim 500$ times (or $\sim 5 M_\oplus$) the present Kuiper Belt mass and a depletion timescale of 300 Myr. This timescale is longer than might be expected for the instability period, but this longer timescale partially compensates for the lower collision rates we compute based on the current structure of the Kuiper Belt as compared to the more compact form in the Nice model. This scenario could shed light on the efficiency of mass removal versus the collisional grinding of TNOs. That is, if the mass depletion rate was higher than the catastrophic collision lifetime of planetesimals, then most remaining TNOs must be collisionally pristine.

4. Case_x2000_300MYR—This is the same as in case_x500_300MYR, except the initial planetesimal mass immediately following the giant planet instability was $\sim 10-20 M_\oplus$.

5. Case_H85_07_300MYR—This case explores the collisional evolution (catastrophic versus cratering) of planetesimals if the slope of the SFD above the knee ($H_r > 8.5$) was steeper ($\alpha = 0.7$) compared to the present-day value of $\alpha = 0.45$. Note that the SFD of planetesimals brighter than the knee is assumed to be the present one. The steeper power-law slope above the knee effectively increases the number of objects with radii $r \lesssim 50$ km, which, on the other hand, increases the total collisional
<table><caption>**Table 3**

Different Cases for the Collision Rates within the TNO Populations</caption><thead><tr><th></th><th>P.E.<sup>a</sup></th><th>α<sub>cold</sub><sup>b</sup></th><th>α<sub>hot</sub><sup>b</sup></th><th>T<br>(Gyr)</th><th>ρ<br>(kg m<sup>-3</sup>)</th><th>ν<br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Nominal case</td><td>Nominal</td><td>Nominal</td><td>Nominal</td><td>4.0</td><td>500</td><td>3.0</td></tr><tr><td>Case_x200_4GYR</td><td>x200</td><td>Nominal</td><td>Nominal</td><td>4.0</td><td>500</td><td>3.0</td></tr><tr><td>Case_x500_300MYR</td><td>x500</td><td>Nominal</td><td>Nominal</td><td>0.3</td><td>500</td><td>3.0</td></tr><tr><td>Case_x2000_300MYR</td><td>x2000</td><td>Nominal</td><td>Nominal</td><td>0.3</td><td>500</td><td>3.0</td></tr><tr><td>Case_H85_06_4GYR</td><td>Nominal</td><td>{10, 10, 10, 1.1, 0.7, 0.6, 0.15}</td><td>{10, 0.2, 0.6, 1.0, 0.7, 0.6, 0.15}</td><td>4.0</td><td>500</td><td>3.0</td></tr><tr><td>Case_H85_07_300MYR</td><td>Nominal</td><td>{10, 10, 10, 1.1, 0.7, 0.7, 0.15}</td><td>{10, 0.2, 0.6, 1.0, 0.7, 0.7, 0.15}</td><td>4.0</td><td>500</td><td>3.0</td></tr></tbody></table>

**Notes.** T—integration time of the collision rates in Gyr. ρ—assumed bulk density for impactors and targets in kg m⁻³. ν—assumed relative collision speeds in km s⁻¹.

a P. E.—Population estimate N(Hr < 8.5) of different TNO classes based on the OSSOS++ ensemble sample (see Table 1 in Section 2). P. E. = x200 and x2000 correspond to an upscale of the nominal N(Hr < 8.5) by a factor of 200 and 2000, respectively.

b Nominal values of the power-law slopes of the SFD for the cold (α<sub>cold</sub> = {10, 10, 10, 1.1, 0.7, 0.45, 0.15}) and excited (α<sub>hot</sub> = {10, 0.2, 0.6, 1.0, 0.7, 0.45, 0.15}) TNO populations (see also Figure 1).

cross section. This orbital and SFD configuration is modeled as lasting 300 Myr. This model explores whether the increase of the number of smaller TNOs (smaller than r = 50 km) would lead to catastrophic disruptions over moderate timescales. This scenario assumes that the bulk of the planetesimal mass was scattered away immediately following the giant planet instability, leaving behind only the present-day mass of the Kuiper Belt but with a slightly different SFD.

6. Case_H85_06_4GYR—This is the same as in case_-H85_07_300MYR, except the power-law slope above the knee was assumed to be α = 0.6, but collisions lasted over 4 Gyr. Although TNO surveys indicate a shallower power-law slope (α = 0.45) above the knee, this case informs whether collisions over a timescale equal to the age of the solar system and the slightly greater number of TNOs with r < 50 km are a sufficient condition to trigger collisional grinding of planetesimals.

For all collision scenario cases, we first calculate the radius r of the largest impactor expected to hit a target of radius R and the specific impact energies Qs over some timescale. For each case and each target TNO belonging to a specific TNO population, we obtain an ensemble of functions r = r(R) and Qs = Qs(R) for the specific impact energies. When compared to SPH simulations, this informs us about the expected outcome from these collisions, e.g., cratering versus complete destruction of the target.

## 4. Results

Figures 2 and 3 show the distribution of the radius r of the largest impactor from a given TNO population that can hit a target of radius R from another TNO population for our nominal case (see Table 3). Superimposed are the theoretical curves for the radius r of the largest impactor required to catastrophically disrupt a target of radius R according to numerical simulations by Benz & Asphaug (1999) and Leinhardt & Stewart (2012). We recall that, in the derivation of the theoretical curves, a bulk density of ρ = 500 kg m⁻³ for both impactors and targets was assumed, as well as a uniform collision speed of v = 3 km s⁻¹. It is evident that given the current architecture of the Kuiper Belt, the impactors hitting a given target of radius R are too small (by almost 2 orders of magnitude) to result in catastrophic disruption, even for collisions over 4 Gyr. This strongly suggests that collisions in

the present-day configuration of the Kuiper Belt are likely to be in the cratering regime rather than the disruptive regime. Moreover, if the initial Kuiper Belt mass was the present one, the SFD of TNOs under these assumptions must reflect the primordial planetesimal accretion SFD. We note that our analysis ignores the collisional evolution during the dynamical instability phase itself. This assumption is likely valid for the “cold” Kuiper Belt, which appears to have formed in situ. For the “hot” objects, the inner disk mass was likely substantially higher and is explored in our other cases. The period of dense orbital configuration for the hot objects likely lasted for T ≲ 200 Myr (e.g., Gomes et al. 2005; Nesvorný et al. 2018), with more recent studies reporting even shorter timescales of ~10–30 Myr (e.g., Nesvorný & Morbidelli 2012; Nesvorný et al. 2018; de Sousa et al. 2020).

Figures 4 and 5 are similar to Figures 2 and 3 but show, for each of our model cases in Table 3, the weighted average radius r of the largest impactor from all TNO populations expected to collide with a target of radius R. The weighting has been performed on the basis of the population estimate of different dynamical TNO groups. As in Figures 2 and 3, collisions do not appear to be in a destruction regime for our nominal case. Some of our exploratory cases, however, do reach the collisional regime.

From Figures 2 and 3, we see that, for collisions to initiate cascading fragmentation, the mass of the Kuiper Belt needs to be 200–2000 times more massive (2–20 M⊕) than current. In the x200_4GYR and x2000_300MYR scenarios, TNOs with sizes R ≲ 40 km can be collisionally evolved for a uniform collision speed of v = 3 km s⁻¹, whereas collisions with v = 0.5 km s⁻¹ shift the transition radius R down to ~12 km. This implies that TNOs with R ≳ 40 km must be collisionally pristine, even for these more massive belt scenarios. The critical radius at which catastrophic disruptions are possible for the detached and scattering populations corresponds to much lower sizes of R ≲ 20 and 10 km, respectively, implying that TNOs with R ≳ 10–20 km do not experience collisional destruction while on these orbits, even if the initial masses in these populations were 200–2000 times larger than seen today.

In scenario case_x500_300MYR, the mass of the early Kuiper Belt is assumed to be ~5–10 M⊕, and collisions lasted for only 300 Myr. In this scenario, the targets of all TNO populations barely reach the disruption threshold of Benz & Asphaug (1999) and Leinhardt & Stewart (2012). Similar conclusions can be drawn for case_H85_06_4GYR and
<table><caption>**Table 3**

Different Cases for the Collision Rates within the TNO Populations</caption><thead><tr><th></th><th>P.E.<sup>a</sup></th><th>α<sub>cold</sub><sup>b</sup></th><th>α<sub>hot</sub><sup>b</sup></th><th>T<br>(Gyr)</th><th>ρ<br>(kg m<sup>-3</sup>)</th><th>ν<br>(km s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Nominal case</td><td>Nominal</td><td>Nominal</td><td>Nominal</td><td>4.0</td><td>500</td><td>3.0</td></tr><tr><td>Case_x200_4GYR</td><td>x200</td><td>Nominal</td><td>Nominal</td><td>4.0</td><td>500</td><td>3.0</td></tr><tr><td>Case_x500_300MYR</td><td>x500</td><td>Nominal</td><td>Nominal</td><td>0.3</td><td>500</td><td>3.0</td></tr><tr><td>Case_x2000_300MYR</td><td>x2000</td><td>Nominal</td><td>Nominal</td><td>0.3</td><td>500</td><td>3.0</td></tr><tr><td>Case_H85_06_4GYR</td><td>Nominal</td><td>{10, 10, 10, 1.1, 0.7, 0.6, 0.15}</td><td>{10, 0.2, 0.6, 1.0, 0.7, 0.6, 0.15}</td><td>4.0</td><td>500</td><td>3.0</td></tr><tr><td>Case_H85_07_300MYR</td><td>Nominal</td><td>{10, 10, 10, 1.1, 0.7, 0.7, 0.15}</td><td>{10, 0.2, 0.6, 1.0, 0.7, 0.7, 0.15}</td><td>4.0</td><td>500</td><td>3.0</td></tr></tbody></table>

**Notes.** T—integration time of the collision rates in Gyr. ρ—assumed bulk density for impactors and targets in kg m⁻³. ν—assumed relative collision speeds in km s⁻¹.

a P. E.—Population estimate N(Hr < 8.5) of different TNO classes based on the OSSOS++ ensemble sample (see Table 1 in Section 2). P. E. = x200 and x2000 correspond to an upscale of the nominal N(Hr < 8.5) by a factor of 200 and 2000, respectively.

b Nominal values of the power-law slopes of the SFD for the cold (α<sub>cold</sub> = {10, 10, 10, 1.1, 0.7, 0.45, 0.15}) and excited (α<sub>hot</sub> = {10, 0.2, 0.6, 1.0, 0.7, 0.45, 0.15}) TNO populations (see also Figure 1).

cross section. This orbital and SFD configuration is modeled as lasting 300 Myr. This model explores whether the increase of the number of smaller TNOs (smaller than r = 50 km) would lead to catastrophic disruptions over moderate timescales. This scenario assumes that the bulk of the planetesimal mass was scattered away immediately following the giant planet instability, leaving behind only the present-day mass of the Kuiper Belt but with a slightly different SFD.

6. Case_H85_06_4GYR—This is the same as in case_-H85_07_300MYR, except the power-law slope above the knee was assumed to be α = 0.6, but collisions lasted over 4 Gyr. Although TNO surveys indicate a shallower power-law slope (α = 0.45) above the knee, this case informs whether collisions over a timescale equal to the age of the solar system and the slightly greater number of TNOs with r < 50 km are a sufficient condition to trigger collisional grinding of planetesimals.

For all collision scenario cases, we first calculate the radius r of the largest impactor expected to hit a target of radius R and the specific impact energies Qs over some timescale. For each case and each target TNO belonging to a specific TNO population, we obtain an ensemble of functions r = r(R) and Qs = Qs(R) for the specific impact energies. When compared to SPH simulations, this informs us about the expected outcome from these collisions, e.g., cratering versus complete destruction of the target.

## 4. Results

Figures 2 and 3 show the distribution of the radius r of the largest impactor from a given TNO population that can hit a target of radius R from another TNO population for our nominal case (see Table 3). Superimposed are the theoretical curves for the radius r of the largest impactor required to catastrophically disrupt a target of radius R according to numerical simulations by Benz & Asphaug (1999) and Leinhardt & Stewart (2012). We recall that, in the derivation of the theoretical curves, a bulk density of ρ = 500 kg m⁻³ for both impactors and targets was assumed, as well as a uniform collision speed of v = 3 km s⁻¹. It is evident that given the current architecture of the Kuiper Belt, the impactors hitting a given target of radius R are too small (by almost 2 orders of magnitude) to result in catastrophic disruption, even for collisions over 4 Gyr. This strongly suggests that collisions in

the present-day configuration of the Kuiper Belt are likely to be in the cratering regime rather than the disruptive regime. Moreover, if the initial Kuiper Belt mass was the present one, the SFD of TNOs under these assumptions must reflect the primordial planetesimal accretion SFD. We note that our analysis ignores the collisional evolution during the dynamical instability phase itself. This assumption is likely valid for the “cold” Kuiper Belt, which appears to have formed in situ. For the “hot” objects, the inner disk mass was likely substantially higher and is explored in our other cases. The period of dense orbital configuration for the hot objects likely lasted for T ≲ 200 Myr (e.g., Gomes et al. 2005; Nesvorný et al. 2018), with more recent studies reporting even shorter timescales of ~10–30 Myr (e.g., Nesvorný & Morbidelli 2012; Nesvorný et al. 2018; de Sousa et al. 2020).

Figures 4 and 5 are similar to Figures 2 and 3 but show, for each of our model cases in Table 3, the weighted average radius r of the largest impactor from all TNO populations expected to collide with a target of radius R. The weighting has been performed on the basis of the population estimate of different dynamical TNO groups. As in Figures 2 and 3, collisions do not appear to be in a destruction regime for our nominal case. Some of our exploratory cases, however, do reach the collisional regime.

From Figures 2 and 3, we see that, for collisions to initiate cascading fragmentation, the mass of the Kuiper Belt needs to be 200–2000 times more massive (2–20 M⊕) than current. In the x200_4GYR and x2000_300MYR scenarios, TNOs with sizes R ≲ 40 km can be collisionally evolved for a uniform collision speed of v = 3 km s⁻¹, whereas collisions with v = 0.5 km s⁻¹ shift the transition radius R down to ~12 km. This implies that TNOs with R ≳ 40 km must be collisionally pristine, even for these more massive belt scenarios. The critical radius at which catastrophic disruptions are possible for the detached and scattering populations corresponds to much lower sizes of R ≲ 20 and 10 km, respectively, implying that TNOs with R ≳ 10–20 km do not experience collisional destruction while on these orbits, even if the initial masses in these populations were 200–2000 times larger than seen today.

In scenario case_x500_300MYR, the mass of the early Kuiper Belt is assumed to be ~5–10 M⊕, and collisions lasted for only 300 Myr. In this scenario, the targets of all TNO populations barely reach the disruption threshold of Benz & Asphaug (1999) and Leinhardt & Stewart (2012). Similar conclusions can be drawn for case_H85_06_4GYR and
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Nominal case (current mass of the Kuiper Belt) for the radius *r* of the largest impactor belonging to a particular TNO population expected to impact a particular target of radius *R* (different subplots) over 4 Gyr. All of our modeled curves (colored solid lines) are estimated for a uniform impact speed of *v* = 3.0 km s⁻¹. The gray and cyan dashed lines indicate the radius *r* of the impactor required to collisionally disrupt a target of radius *R*, derived from Benz & Asphaug (1999), for collisions with impact speeds of *v* = 0.5 and 3.0 km s⁻¹, respectively. Similarly, the yellow and blue dashed lines correspond to the radius *r* of an impactor required to fragment a target of radius *R* according to Leinhardt & Stewart (2012) for two different impact speeds. All experimental curves are derived for icy bodies with an assumed bulk density of *ρ* = 500 kg m⁻³.</figcaption></figure>

case_H85_07_300MYR, in which we change the nominal power slope of the SFD for $H_r > 8.5$ (Figure 1) from $\alpha = 0.45$ to $\alpha = 0.6$ and $0.7$, effectively accounting for a greater number of TNOs with $5 \text{ km} < R < 50 \text{ km}$ (see Section 2). Regardless of the integration length (4 Gyr versus 300 Myr), the specific impact energies remain low for both cases.

It is worth noting that if the theoretical curves are $v = 0.5 \text{ km s}^{-1}$ for both Benz & Asphaug (1999) and Leinhardt & Stewart (2012), the impactor radius, required to disrupt a target of radius *R* beyond 100 km, should be interpreted with caution. Especially for *R* ~ 1000 km, the expected impactor radius that can break apart the target appears to be greater than *R*. The reason for that is that the Benz & Asphaug (1999) simulations span ranges of *R* < 100 km, and extrapolations to larger sizes may be inappropriate.

The scenarios considering a large primordial Kuiper Belt mass are difficult to reconcile with the catastrophic destruction survival of wide binary systems that make up a significant fraction of the cold classicals. The cold classical objects are likely to have formed near 40 au from the Sun, and the presence of more impactors in that region would have increased the collision rates and significantly depleted the wide binary population (Parker & Kavelaars 2010).

Our results demand explanation of the following puzzling questions.

i. If the mass of the planetesimal belt beyond Pluto was ~2 $M_⊕$ after the dynamical instability, that requires that the collisions must have persisted for 4 Gyr (case_-x200_4GYR). However, then TNOs fainter than the “knee” in the SFD would have acquired a collisional equilibrium slope and cratering records on the surface of Arroloth, though Spencer et al. (2020) found no evidence of a such collisional steady state. In fact, New Horizons finds a deficit of small impactors on both Pluto and Charon and the small CCKBO Arroloth. Unless there was a major resurfacing process that has erased old craters and the “footprint” of early TNO SFD, it is reasonable to accept that the small TNOs experienced very few catastrophic collisions.

ii. If the mass of the primordial Kuiper Belt was ~20 $M_⊕$ and collisions lasted for only 300 Myr (case_-x2000_300MYR), then the SFD of TNOs with $R < 40$ km may have reached collisional equilibrium. However, once again, crater observations by the New Horizons mission do not support this outcome.

As an alternative way of displaying our results, we calculate the specific impact energy distribution $Q_s$ for the largest impactors that can hit a given target of radius *R* for the cases described in Table 3 and compare them to the curves derived
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Nominal case (current mass of the Kuiper Belt) for the radius *r* of the largest impactor belonging to a particular TNO population expected to impact a particular target of radius *R* (different subplots) over 4 Gyr. All of our modeled curves (colored solid lines) are estimated for a uniform impact speed of *v* = 3.0 km s⁻¹. The gray and cyan dashed lines indicate the radius *r* of the impactor required to collisionally disrupt a target of radius *R*, derived from Benz & Asphaug (1999), for collisions with impact speeds of *v* = 0.5 and 3.0 km s⁻¹, respectively. Similarly, the yellow and blue dashed lines correspond to the radius *r* of an impactor required to fragment a target of radius *R* according to Leinhardt & Stewart (2012) for two different impact speeds. All experimental curves are derived for icy bodies with an assumed bulk density of *ρ* = 500 kg m⁻³.</figcaption></figure>

case_H85_07_300MYR, in which we change the nominal power slope of the SFD for $H_r > 8.5$ (Figure 1) from $\alpha = 0.45$ to $\alpha = 0.6$ and $0.7$, effectively accounting for a greater number of TNOs with $5 \text{ km} < R < 50 \text{ km}$ (see Section 2). Regardless of the integration length (4 Gyr versus 300 Myr), the specific impact energies remain low for both cases.

It is worth noting that if the theoretical curves are $v = 0.5 \text{ km s}^{-1}$ for both Benz & Asphaug (1999) and Leinhardt & Stewart (2012), the impactor radius, required to disrupt a target of radius *R* beyond 100 km, should be interpreted with caution. Especially for *R* ~ 1000 km, the expected impactor radius that can break apart the target appears to be greater than *R*. The reason for that is that the Benz & Asphaug (1999) simulations span ranges of *R* < 100 km, and extrapolations to larger sizes may be inappropriate.

The scenarios considering a large primordial Kuiper Belt mass are difficult to reconcile with the catastrophic destruction survival of wide binary systems that make up a significant fraction of the cold classicals. The cold classical objects are likely to have formed near 40 au from the Sun, and the presence of more impactors in that region would have increased the collision rates and significantly depleted the wide binary population (Parker & Kavelaars 2010).

Our results demand explanation of the following puzzling questions.

i. If the mass of the planetesimal belt beyond Pluto was ~2 $M_⊕$ after the dynamical instability, that requires that the collisions must have persisted for 4 Gyr (case_-x200_4GYR). However, then TNOs fainter than the “knee” in the SFD would have acquired a collisional equilibrium slope and cratering records on the surface of Arroloth, though Spencer et al. (2020) found no evidence of a such collisional steady state. In fact, New Horizons finds a deficit of small impactors on both Pluto and Charon and the small CCKBO Arroloth. Unless there was a major resurfacing process that has erased old craters and the “footprint” of early TNO SFD, it is reasonable to accept that the small TNOs experienced very few catastrophic collisions.

ii. If the mass of the primordial Kuiper Belt was ~20 $M_⊕$ and collisions lasted for only 300 Myr (case_-x2000_300MYR), then the SFD of TNOs with $R < 40$ km may have reached collisional equilibrium. However, once again, crater observations by the New Horizons mission do not support this outcome.

As an alternative way of displaying our results, we calculate the specific impact energy distribution $Q_s$ for the largest impactors that can hit a given target of radius *R* for the cases described in Table 3 and compare them to the curves derived
<figure><img src="image_5.png" /><figcaption>Figure 3. Same as Figure 2.</figcaption></figure>

from numerical simulations by Benz & Asphaug (1999) and Leinhardt & Stewart (2012). These results are presented in Figures 6 and 7 as the weighted average $Q_s$ for all impactor populations, where the weighting has been performed based on the ICP between different TNO populations and the relative population estimate for each impactor population. It is evident that the specific impact energies for the nominal case are well below the catastrophic disruption threshold $Q_D^*$ curves. Once again, that argues against a collisional reshaping of the Kuiper Belt for our nominal case. In order for collisional cascade to commence, the mass of the early Kuiper disk must have been ~20 $M_⊕$ if collisions lasted for 300 Myr and ~2 $M_⊕$ if collisions persisted for 4 Gyr. It turns out that an early Kuiper Belt mass of ~5–10 $M_⊕$, case_x500_300MYR, is required to barely trigger cascading fragmentation for objects with $R \lesssim 12$ km, even if collisions persisted for only 300 Myr, leaving larger objects collisionally unevolved. Recall that these timescales are long compared to the typical modeling of Nice model–like instabilities, but in those configurations, the more tightly packed disk is likely to have evolved more quickly; thus, our long timescales provide an indication of evolution that is likely to have occurred early in the instability period.

The location of the knee in the SFD of TNO populations is well constrained by surveys and occurs at $R \sim 50$ km. Changing the slope of the SFD (case_H85_06_4GYR and case_H85_07_300MYR) increases the probable maximum specific impact energies $Q_s$ owing to a greater number of impactors available to collide with a target of radius $R$, though

the impact energy remains below $Q_D^*$, and collisional disruption is insignificant for these scenarios.

Figures 8 and 9 present the expected weighted average number, $N_c$, of catastrophic impacts that a target of radius $R$ from a particular TNO population (different panels) experiences under the different cases explored in this study. To calculate these functions, we utilized the catastrophic disruption energy model from Benz & Asphaug (1999), which yields similar results to Leinhardt & Stewart (2012; see, e.g., Figures 4 and 6). It is evident that in the nominal case, only about 2.5% of the cold classical, hot, inner, and main resonant TNO populations experience catastrophic collisions over the age of the solar system. On the other hand, case_x2000_300MYR and case_x200_4GYR suggest that all CCKBOs with $R \lesssim 40–50$ km must be fragments from previous collisions. Similar conclusions can be drawn for the hot, inner, and main resonant TNOs. The radius $R$ slowly transitions to smaller values for the outer resonant and detached populations, reaching $R \sim 10$ km for the scattering TNOs. As before, these cases require an early Kuiper Belt of mass $2M_⊕ > M_{KB} < 20M_⊕$, which is difficult to reconcile given the present low mass of planetesimals beyond Neptune’s orbit and the low crater frequency on the Pluto–Charon system (Singer et al. 2019) and the CCKBO Arrokoth (Spencer et al. 2020).

The increase in the number of impactors with $5 \text{ km} < R < 50 \text{ km}$ (case_H85_06_4GYR and case_H85_07_300MYR) sharply increases the number of catastrophic collisions, $N_c$, implying that ~20%–50% of the cold classical, hot, inner, and main resonant targets with radii $R \lesssim 10$ km would have
<figure><img src="image_7.png" /><figcaption>Figure 3. Same as Figure 2.</figcaption></figure>

from numerical simulations by Benz & Asphaug (1999) and Leinhardt & Stewart (2012). These results are presented in Figures 6 and 7 as the weighted average $Q_s$ for all impactor populations, where the weighting has been performed based on the ICP between different TNO populations and the relative population estimate for each impactor population. It is evident that the specific impact energies for the nominal case are well below the catastrophic disruption threshold $Q_D^*$ curves. Once again, that argues against a collisional reshaping of the Kuiper Belt for our nominal case. In order for collisional cascade to commence, the mass of the early Kuiper disk must have been ~20 $M_⊕$ if collisions lasted for 300 Myr and ~2 $M_⊕$ if collisions persisted for 4 Gyr. It turns out that an early Kuiper Belt mass of ~5–10 $M_⊕$, case_x500_300MYR, is required to barely trigger cascading fragmentation for objects with $R \lesssim 12$ km, even if collisions persisted for only 300 Myr, leaving larger objects collisionally unevolved. Recall that these timescales are long compared to the typical modeling of Nice model–like instabilities, but in those configurations, the more tightly packed disk is likely to have evolved more quickly; thus, our long timescales provide an indication of evolution that is likely to have occurred early in the instability period.

The location of the knee in the SFD of TNO populations is well constrained by surveys and occurs at $R \sim 50$ km. Changing the slope of the SFD (case_H85_06_4GYR and case_H85_07_300MYR) increases the probable maximum specific impact energies $Q_s$ owing to a greater number of impactors available to collide with a target of radius $R$, though

the impact energy remains below $Q_D^*$, and collisional disruption is insignificant for these scenarios.

Figures 8 and 9 present the expected weighted average number, $N_c$, of catastrophic impacts that a target of radius $R$ from a particular TNO population (different panels) experiences under the different cases explored in this study. To calculate these functions, we utilized the catastrophic disruption energy model from Benz & Asphaug (1999), which yields similar results to Leinhardt & Stewart (2012; see, e.g., Figures 4 and 6). It is evident that in the nominal case, only about 2.5% of the cold classical, hot, inner, and main resonant TNO populations experience catastrophic collisions over the age of the solar system. On the other hand, case_x2000_300MYR and case_x200_4GYR suggest that all CCKBOs with $R \lesssim 40–50$ km must be fragments from previous collisions. Similar conclusions can be drawn for the hot, inner, and main resonant TNOs. The radius $R$ slowly transitions to smaller values for the outer resonant and detached populations, reaching $R \sim 10$ km for the scattering TNOs. As before, these cases require an early Kuiper Belt of mass $2M_⊕ > M_{KB} < 20M_⊕$, which is difficult to reconcile given the present low mass of planetesimals beyond Neptune’s orbit and the low crater frequency on the Pluto–Charon system (Singer et al. 2019) and the CCKBO Arrokoth (Spencer et al. 2020).

The increase in the number of impactors with $5 \text{ km} < R < 50 \text{ km}$ (case_H85_06_4GYR and case_H85_07_300MYR) sharply increases the number of catastrophic collisions, $N_c$, implying that ~20%–50% of the cold classical, hot, inner, and main resonant targets with radii $R \lesssim 10$ km would have
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Average radius *r* of the largest impactor that can collide with a target of radius *R* over a range of collision cases. All curves (colored solid lines) are weighted by the relative population estimate of each TNO impactor population. The radius of the largest impactor, by definition, is *r* ≤ *R*. The different panels denote targets belonging to different TNO populations. The dashed lines indicate the minimum radius *r* of an impactor required to catastrophically disrupt a target body of radius *R* based on the analysis of Benz & Asphaug (1999) and Leinhardt & Stewart (2012) for two different impact speeds of *v* = 0.5 and 3 km s⁻¹. All four theoretical curves have been calculated for icy bodies with an assumed bulk density of *ρ* = 500 kg m⁻³.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Same as Figure 4.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 4.** Average radius *r* of the largest impactor that can collide with a target of radius *R* over a range of collision cases. All curves (colored solid lines) are weighted by the relative population estimate of each TNO impactor population. The radius of the largest impactor, by definition, is *r* ≤ *R*. The different panels denote targets belonging to different TNO populations. The dashed lines indicate the minimum radius *r* of an impactor required to catastrophically disrupt a target body of radius *R* based on the analysis of Benz & Asphaug (1999) and Leinhardt & Stewart (2012) for two different impact speeds of *v* = 0.5 and 3 km s⁻¹. All four theoretical curves have been calculated for icy bodies with an assumed bulk density of *ρ* = 500 kg m⁻³.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 5.** Same as Figure 4.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Specific impact energy $Q_s$ of collision between a target of radius $R$ and the weighted average largest impactor of radius $r \le R$ it can collide with for a range of different cases in Table 3. All cases are calculated for a uniform impact speed of $v = 3$ km s$^{-1}$. The different panels denote targets belonging to different TNO populations, with different colors indicating different cases. The dashed lines denote the disruption energy threshold $Q_D^p$ for a target of radius $R$ based on the analysis of Benz & Asphaug (1999) and Leinhardt & Stewart (2012) for two different collision speeds, $v = 0.5$ and $3$ km s$^{-1}$. All four theoretical curves are derived for ice, with an assumed bulk density of $\rho = 500$ kg m$^{-3}$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Same as Figure 6.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Specific impact energy $Q_s$ of collision between a target of radius $R$ and the weighted average largest impactor of radius $r \le R$ it can collide with for a range of different cases in Table 3. All cases are calculated for a uniform impact speed of $v = 3$ km s$^{-1}$. The different panels denote targets belonging to different TNO populations, with different colors indicating different cases. The dashed lines denote the disruption energy threshold $Q_D^p$ for a target of radius $R$ based on the analysis of Benz & Asphaug (1999) and Leinhardt & Stewart (2012) for two different collision speeds, $v = 0.5$ and $3$ km s$^{-1}$. All four theoretical curves are derived for ice, with an assumed bulk density of $\rho = 500$ kg m$^{-3}$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Same as Figure 6.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Weighted average number of catastrophic impacts $N_c$ (colored lines) a single target of radius $R$ experiences with an impactor of radius $r \le R$ for a range of different collision cases in Table 3. Each panel refers to a target body belonging to a particular TNO population. The dashed horizontal line corresponds to $N_c = 1$. (Author note—In calculating the curves for all cases, the catastrophic disruption threshold of Benz & Asphaug (1999) was used. Using the disruption threshold formalism of Leinhardt & Stewart (2012) does not change the results, as they produce similar disruption threshold curves for a given material and impact speed; see previous figures.) For all cases, a uniform impact speed of $v = 3 \text{ km s}^{-1}$ and target bulk density of $\rho = 500 \text{ kg m}^{-3}$ were assumed.</figcaption></figure>

experienced at least one catastrophic collision for these scenarios, while $R \gtrsim 10$ km planetesimals must be collisionally pristine. However, these two cases fail to match the observed location of the knee in the TNO SFD and thus could be ruled out as plausible initial Kuiper Belt structure.

**5. Discussion and Conclusions**

In this work, we investigate the rate of collisions in the trans-Neptunian region with implications for past collisional erosion of the KBOs. As an input to our modeling, we use the SFD of TNOs from the OSSOS++ ensemble sample and the ICPs reported in Abedin et al. (2021). We then calculate the radius $r$ of the largest impactor that has a unit probability of collision with a given target of radius $R$ over some time interval $\Delta T$. From this, and assuming a uniform impact speed of $3 \text{ km s}^{-1}$, we determine the largest specific impact energy $Q_s$ between each impactor–target pair considered here, and we then compare that to the catastrophic disruption energy threshold $Q_D^*$ determined from SPH simulations by Benz & Asphaug (1999) and Leinhardt & Stewart (2012). Furthermore, due to a lack of precise knowledge of the primordial (during and after giant planet instability and migration) Kuiper Belt mass, size, and orbital distribution, we investigate different cases with different TNO population estimates and SFD, with all cases based on the present OSSOS++ orbital model of the Kuiper Belt. That is, we only model the ICPs after the excited/hot TNO populations have been emplaced into their present orbits.

Based on our analysis, we find no evidence for catastrophic disruption going on in the present Kuiper Belt. In fact, the

collision frequencies in the current populations are so low that the size of most impactors colliding with a target of radius $R$ are well below the impact disruption threshold (Figures 4–7). Likely most of these collisions are in the crater-forming mode, rather than in a disruptive regime. That could also be inferred from Figures 8 and 9, which imply that under the present orbital and SFD of the KBOs, only $\lesssim 2.5\%$ of the initial cold classicals with radii $R \sim 10$ km experienced catastrophic disruption. This conclusion also holds for all dynamical subpopulations considered here. Catastrophic collisions are absent from the present-day Kuiper Belt. There may have been, however, catastrophic disruption early on in the Kuiper Belt’s history. Collisional destruction of TNOs, according to our analysis, is possible if the early Kuiper Belt was as massive as $2 M_\oplus$ and collisions persisted for over 4 Gyr, but these collisions would not be sufficient to remove 2 orders of magnitude of mass from the region. Collisional fragmentation for objects in the $5 \text{ km} < r < 30 \text{ km}$ size range, depending on the impactor speed, would occur if the primordial mass was at least $20 M_\oplus$ and collisions between planetesimals lasted only $\sim 300$ Myr (case_x200_4GYR and case_x2000_300MYR), a situation similar to the conditions in the Nice model (Gomes et al. 2005; Morbidelli et al. 2005; Tsiganis et al. 2005) type scenarios. Another case leading to catastrophic destruction is case_x500_300MYR, with a mass of the primordial Kuiper Belt equivalent to $5–10 M_\oplus$. This case indicates that the critical size of TNOs below which objects must be collisional fragments occurs around $R \sim 10$ km. In this scenario, for example, we determine that at least $\sim 10\%$ of the cold classical,
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Weighted average number of catastrophic impacts $N_c$ (colored lines) a single target of radius $R$ experiences with an impactor of radius $r \le R$ for a range of different collision cases in Table 3. Each panel refers to a target body belonging to a particular TNO population. The dashed horizontal line corresponds to $N_c = 1$. (Author note—In calculating the curves for all cases, the catastrophic disruption threshold of Benz & Asphaug (1999) was used. Using the disruption threshold formalism of Leinhardt & Stewart (2012) does not change the results, as they produce similar disruption threshold curves for a given material and impact speed; see previous figures.) For all cases, a uniform impact speed of $v = 3 \text{ km s}^{-1}$ and target bulk density of $\rho = 500 \text{ kg m}^{-3}$ were assumed.</figcaption></figure>

experienced at least one catastrophic collision for these scenarios, while $R \gtrsim 10$ km planetesimals must be collisionally pristine. However, these two cases fail to match the observed location of the knee in the TNO SFD and thus could be ruled out as plausible initial Kuiper Belt structure.

**5. Discussion and Conclusions**

In this work, we investigate the rate of collisions in the trans-Neptunian region with implications for past collisional erosion of the KBOs. As an input to our modeling, we use the SFD of TNOs from the OSSOS++ ensemble sample and the ICPs reported in Abedin et al. (2021). We then calculate the radius $r$ of the largest impactor that has a unit probability of collision with a given target of radius $R$ over some time interval $\Delta T$. From this, and assuming a uniform impact speed of $3 \text{ km s}^{-1}$, we determine the largest specific impact energy $Q_s$ between each impactor–target pair considered here, and we then compare that to the catastrophic disruption energy threshold $Q_D^*$ determined from SPH simulations by Benz & Asphaug (1999) and Leinhardt & Stewart (2012). Furthermore, due to a lack of precise knowledge of the primordial (during and after giant planet instability and migration) Kuiper Belt mass, size, and orbital distribution, we investigate different cases with different TNO population estimates and SFD, with all cases based on the present OSSOS++ orbital model of the Kuiper Belt. That is, we only model the ICPs after the excited/hot TNO populations have been emplaced into their present orbits.

Based on our analysis, we find no evidence for catastrophic disruption going on in the present Kuiper Belt. In fact, the

collision frequencies in the current populations are so low that the size of most impactors colliding with a target of radius $R$ are well below the impact disruption threshold (Figures 4–7). Likely most of these collisions are in the crater-forming mode, rather than in a disruptive regime. That could also be inferred from Figures 8 and 9, which imply that under the present orbital and SFD of the KBOs, only $\lesssim 2.5\%$ of the initial cold classicals with radii $R \sim 10$ km experienced catastrophic disruption. This conclusion also holds for all dynamical subpopulations considered here. Catastrophic collisions are absent from the present-day Kuiper Belt. There may have been, however, catastrophic disruption early on in the Kuiper Belt’s history. Collisional destruction of TNOs, according to our analysis, is possible if the early Kuiper Belt was as massive as $2 M_\oplus$ and collisions persisted for over 4 Gyr, but these collisions would not be sufficient to remove 2 orders of magnitude of mass from the region. Collisional fragmentation for objects in the $5 \text{ km} < r < 30 \text{ km}$ size range, depending on the impactor speed, would occur if the primordial mass was at least $20 M_\oplus$ and collisions between planetesimals lasted only $\sim 300$ Myr (case_x200_4GYR and case_x2000_300MYR), a situation similar to the conditions in the Nice model (Gomes et al. 2005; Morbidelli et al. 2005; Tsiganis et al. 2005) type scenarios. Another case leading to catastrophic destruction is case_x500_300MYR, with a mass of the primordial Kuiper Belt equivalent to $5–10 M_\oplus$. This case indicates that the critical size of TNOs below which objects must be collisional fragments occurs around $R \sim 10$ km. In this scenario, for example, we determine that at least $\sim 10\%$ of the cold classical,
<figure><img src="image_11.png" /><figcaption>Figure 9. Same as Figure 8.</figcaption></figure>

hot, inner resonant, and main resonant TNO populations with radius $R \gtrsim 10$ km are likely to be primordial (not fragments; see Figures 4–9), while 90% are likely collisional fragments. Under this case, only ~10% of the detached and scattering TNO populations with $R \lesssim 10$ km must be collisional fragments. We note, however, that these differences in the collisional disruption sizes are a result of the different orbital distributions. Increasing the slope of the TNO SFD for objects faintward of $H_r = 8.5$ from $\alpha = 0.45$ to 0.6 and even 0.7 (case_-H85_06_4GYR and case_H85_07_300MYR) does not increase the collision rates to a degree where catastrophic destruction may become important, even over a timescale equal to the age of the solar system. Increasing the SFD slope fainter than $H_r = 8.5$ is equivalent to increasing the number density of impactors and targets with radii $\lesssim 50$ km. The reason that does not initiate a collisional cascade is simply because of the low ICPs, given the present orbital structure of the Kuiper Belt, and the transition of the impactor SFD slope to $\alpha \sim 0.15-0.2$

(Singer et al. 2019; Spencer et al. 2020) for objects smaller than ~5 km. In other words, there are not enough objects in the range $5 \text{ km} < R < 50 \text{ km}$ in case_H85_06_4GYR and case_-H85_07_300MYR to initiate a collisional cascade. The shallow slope seen today cannot be the result of collisional evolution within the in situ population and may be a relic of the planetesimal accretion stage.

In Abedin et al. (2021; see also Greenstreet et al. 2015, 2016), we demonstrate that collisions between TNOs are happening with speeds well below $5 \text{ km s}^{-1}$. For instance, collisions between cold classicals are most likely to occur with $v < 0.5 \text{ km s}^{-1}$, whereas the excited TNOs are most likely to collide with $v \lesssim 2 \text{ km s}^{-1}$. Throughout our analysis, we have adopted a uniform collision speed of $v = 3 \text{ km s}^{-1}$, exceeding even the most energetic impacts expected, and thus $v = 3 \text{ km s}^{-1}$ serving as an upper limit of the impact energy for an impactor–target pair. Most certainly, increasing the impact speed to $v = 5 \text{ km s}^{-1}$ would increase the specific impact energies and perhaps bring them closer to the catastrophic disruption threshold, though only $\lesssim 5\%$ of TNOs collide with such high velocities in the present structure of the Kuiper Belt. Furthermore, we have adopted a constant mass density of $\rho = 500 \text{ kg m}^{-3}$ (appropriate for cometary nuclei) in our simulation for both impactors and targets. New Horizons measurements of the density of Charon and Pluto suggest values exceeding $1000 \text{ kg m}^{-3}$, though there would be fewer objects with such density assuming density scales with the mass of the object. However, since the catastrophic disruption threshold as derived by Benz & Asphaug (1999) and Leinhardt & Stewart (2012) depends linearly on the mass density, considering denser TNOs would result in even higher $Q_D^*$.

As a result of our study, we conclude that there are no catastrophic collisions happening in present Kuiper Belt. The specific impact energies between a target from any dynamical TNO class and the largest impactor it can collide with are well below the catastrophic disruption threshold $Q_D^*$. Moreover, as the cold classical TNOs are likely to have formed in situ, their current SFD must resemble the primordial planetesimal accretional stage. We also find that for a Nice model–like scenario, some collisional evolution would have been possible up to sizes as large as 20 km, depending on the amount of time before the instability. If the early Kuiper Belt was as massive as 20 $M_\oplus$, then 100% of the TNO populations with $R \lesssim 40$ km must have experienced at least one catastrophic collision over the assumed time intervals. Comparison of the size distributions of the excited and cold classical populations in the 5–50 km size range may reveal the signature of the timing of the instability.

A successful massive primordial (after giant planet migration) Kuiper Belt initial condition scenario for collisional evolution must be able to explain the presently low mass of the Kuiper Belt, the observed location of the knee in the TNO SFD, and the deficit of small impactors within the cold and excited classical KBOs. However, if the primordial Kuiper Belt was indeed massive and catastrophic collisions reached a steady state, it remains a mystery why this has not been imprinted on the present TNO SFD and the cratering record on the Pluto–Charon system and Arrokoth and how two to three orders of mass (the current Kuiper Belt mass does not lead to collisional destruction of TNOs, and if collisions were important, it must have been $10^2–10^3$ times more massive) was removed from the Kuiper Belt.
<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>Figure 9. Same as Figure 8.</figcaption></figure>

hot, inner resonant, and main resonant TNO populations with radius $R \gtrsim 10$ km are likely to be primordial (not fragments; see Figures 4–9), while 90% are likely collisional fragments. Under this case, only ~10% of the detached and scattering TNO populations with $R \lesssim 10$ km must be collisional fragments. We note, however, that these differences in the collisional disruption sizes are a result of the different orbital distributions. Increasing the slope of the TNO SFD for objects faintward of $H_r = 8.5$ from $\alpha = 0.45$ to 0.6 and even 0.7 (case_-H85_06_4GYR and case_H85_07_300MYR) does not increase the collision rates to a degree where catastrophic destruction may become important, even over a timescale equal to the age of the solar system. Increasing the SFD slope fainter than $H_r = 8.5$ is equivalent to increasing the number density of impactors and targets with radii $\lesssim 50$ km. The reason that does not initiate a collisional cascade is simply because of the low ICPs, given the present orbital structure of the Kuiper Belt, and the transition of the impactor SFD slope to $\alpha \sim 0.15-0.2$

(Singer et al. 2019; Spencer et al. 2020) for objects smaller than ~5 km. In other words, there are not enough objects in the range $5 \text{ km} < R < 50 \text{ km}$ in case_H85_06_4GYR and case_-H85_07_300MYR to initiate a collisional cascade. The shallow slope seen today cannot be the result of collisional evolution within the in situ population and may be a relic of the planetesimal accretion stage.

In Abedin et al. (2021; see also Greenstreet et al. 2015, 2016), we demonstrate that collisions between TNOs are happening with speeds well below $5 \text{ km s}^{-1}$. For instance, collisions between cold classicals are most likely to occur with $v < 0.5 \text{ km s}^{-1}$, whereas the excited TNOs are most likely to collide with $v \lesssim 2 \text{ km s}^{-1}$. Throughout our analysis, we have adopted a uniform collision speed of $v = 3 \text{ km s}^{-1}$, exceeding even the most energetic impacts expected, and thus $v = 3 \text{ km s}^{-1}$ serving as an upper limit of the impact energy for an impactor–target pair. Most certainly, increasing the impact speed to $v = 5 \text{ km s}^{-1}$ would increase the specific impact energies and perhaps bring them closer to the catastrophic disruption threshold, though only $\lesssim 5\%$ of TNOs collide with such high velocities in the present structure of the Kuiper Belt. Furthermore, we have adopted a constant mass density of $\rho = 500 \text{ kg m}^{-3}$ (appropriate for cometary nuclei) in our simulation for both impactors and targets. New Horizons measurements of the density of Charon and Pluto suggest values exceeding $1000 \text{ kg m}^{-3}$, though there would be fewer objects with such density assuming density scales with the mass of the object. However, since the catastrophic disruption threshold as derived by Benz & Asphaug (1999) and Leinhardt & Stewart (2012) depends linearly on the mass density, considering denser TNOs would result in even higher $Q_D^*$.

As a result of our study, we conclude that there are no catastrophic collisions happening in present Kuiper Belt. The specific impact energies between a target from any dynamical TNO class and the largest impactor it can collide with are well below the catastrophic disruption threshold $Q_D^*$. Moreover, as the cold classical TNOs are likely to have formed in situ, their current SFD must resemble the primordial planetesimal accretional stage. We also find that for a Nice model–like scenario, some collisional evolution would have been possible up to sizes as large as 20 km, depending on the amount of time before the instability. If the early Kuiper Belt was as massive as 20 $M_\oplus$, then 100% of the TNO populations with $R \lesssim 40$ km must have experienced at least one catastrophic collision over the assumed time intervals. Comparison of the size distributions of the excited and cold classical populations in the 5–50 km size range may reveal the signature of the timing of the instability.

A successful massive primordial (after giant planet migration) Kuiper Belt initial condition scenario for collisional evolution must be able to explain the presently low mass of the Kuiper Belt, the observed location of the knee in the TNO SFD, and the deficit of small impactors within the cold and excited classical KBOs. However, if the primordial Kuiper Belt was indeed massive and catastrophic collisions reached a steady state, it remains a mystery why this has not been imprinted on the present TNO SFD and the cratering record on the Pluto–Charon system and Arrokoth and how two to three orders of mass (the current Kuiper Belt mass does not lead to collisional destruction of TNOs, and if collisions were important, it must have been $10^2–10^3$ times more massive) was removed from the Kuiper Belt.
K.V. acknowledges funding from the NSF (grant AST-1824869) and NASA (grants 80NSSC19K0785 and 80NSSC21K0376).

## ORCID iDs

Abedin Y. Abedin https://orcid.org/0000-0002-7025-0975

J. J. Kavelaars https://orcid.org/0000-0001-7032-5255

Jean-Marc Petit https://orcid.org/0000-0003-0407-2266

Brett Gladman https://orcid.org/0000-0002-0283-2260

Michele Bannister https://orcid.org/0000-0003-3257-4490

Mike Alexandersen https://orcid.org/0000-0003-4143-8589

Stephen Gwyn https://orcid.org/0000-0001-8221-8406

Kathryn Volk https://orcid.org/0000-0001-8736-236X

## References

Abedin, A. Y., Kavelaars, J. J., Greenstreet, S., et al. 2021, AJ, 161, 195

Alexandersen, M., Gladman, B., Kavelaars, J. J., et al. 2016, AJ, 152, 111

Bannister, M. T., Gladman, B. J., Kavelaars, J. J., et al. 2018, ApJS, 236, 18

Bannister, M. T., Kavelaars, J. J., Petit, J.-M., et al. 2016, AJ, 152, 70

Benz, W., & Asphaug, E. 1999, Icar, 142, 5

Bernstein, G. M., Trilling, D. E., Allen, R. L., et al. 2004, AJ, 128, 1364

Chen, Y.-T., Gladman, B., Volk, K., et al. 2019, AJ, 158, 214

Clement, M. S., Kaib, N. A., Raymond, S. N., Chambers, J. E., & Walsh, K. J. 2019a, Icar, 321, 778

Clement, M. S., Kaib, N. A., Raymond, S. N., & Walsh, K. J. 2018, Icar, 311, 340

Clement, M. S., Raymond, S. N., & Kaib, N. A. 2019b, AJ, 157, 38

DAVIS, D. R., & Farinella, P. 1997, Icar, 125, 50

de Sousa, R. R., Morbidelli, A., Raymond, S. N., et al. 2020, Icar, 339, 113605

Dohnanyi, J. S. 1969, JGR, 74, 2531

Farinella, P., & Davis, D. R. 1996, Sci, 273, 938

Fraser, W. C., Bannister, M. T., Pike, R. E., et al. 2017, NatAs, 1, 0088

Fraser, W. C., Benecchi, S. D., Kavelaars, J. J., et al. 2021, PSJ, 2, 90

Fraser, W. C., Brown, M. E., Morbidelli, A., Parker, A., & Batygin, K. 2014, ApJ, 782, 100

Fraser, W. C., & Kavelaars, J. J. 2009, AJ, 137, 72

Gladman, B., Lawler, S. M., Petit, J.-M., et al. 2012, AJ, 144, 23

Gladman, B., Marsden, B. G., & Vanlaerhoven, C. 2008, Nomenclature in the Outer Solar System (Tucson, AZ: Univ. Arizona Press), 43

Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005, Natur, 435, 466

Greenstreet, S., Gladman, B., & McKinnon, W. B. 2015, Icar, 258, 267

Greenstreet, S., Gladman, B., & McKinnon, W. B. 2016, Icar, 274, 366

Greenstreet, S., Gladman, B., McKinnon, W. B., Kavelaars, J. J., & Singer, K. N. 2019, ApJL, 872, L5

Grundy, W. M., Noll, K. S., Nimmo, F., et al. 2011, Icar, 213, 678

Gwyn, S. D. J. 2012, AJ, 143, 38

Holsapple, K. A. 1993, AREPS, 21, 333

Holsapple, K. A., & Schmidt, R. M. 1987, JGR, 92, 6350

Housen, K. R., & Holsapple, K. A. 1990, Icar, 84, 226

Housen, K. R., Schmidt, R. M., & Holsapple, K. A. 1991, Icar, 94, 180

Housen, K. R., Sweet, W. J., & Holsapple, K. A. 2018, Icar, 300, 72

Jewitt, D., & Luu, J. 1993, Natur, 362, 730

Johansen, A., Jacquet, E., Cuzzi, J. N., Morbidelli, A., & Gounelle, M. 2015, New Paradigms for Asteroid Formation (Tucson, AZ: Univ. Arizona Press), 471

Jutzi, M., Michel, P., Benz, W., & Richardson, D. C. 2009, AAS/DPS Meeting Abstracts, 41, 50.06

Kavelaars, J. J., Jones, R. L., Gladman, B. J., et al. 2009, AJ, 137, 4917

Kavelaars, J. J., Lawler, S. M., Bannister, M. T., & Shankman, C. 2020, Perspectives on the Distribution of Orbits of Distant Trans-Neptunian Objects (Amsterdam: Elsevier), 61

Kavelaars, J. J., Petit, J.-M., Gladman, B., et al. 2021, ApJL, 920, L28

Kenyon, S. J., & Bromley, B. C. 2004, AJ, 128, 1916

Kenyon, S. J., & Bromley, B. C. 2017, ApJ, 839, 38

Kenyon, S. J., & Bromley, B. C. 2020, PSJ, 1, 40

Kenyon, S. J., & Luu, J. X. 1999, AJ, 118, 1101

Kovalenko, I. D., Doressoundiram, A., Lellouch, E., et al. 2017, A&A, 608, A19

Lawler, S. M., Kavelaars, J. J., Alexandersen, M., et al. 2018a, FrASS, 5, 14

Lawler, S. M., Shankman, C., Kavelaars, J. J., et al. 2018b, AJ, 155, 197

Leinhardt, Z. M., & Stewart, S. T. 2009, Icar, 199, 542

Leinhardt, Z. M., & Stewart, S. T. 2012, ApJ, 745, 79

Levison, H. F., Morbidelli, A., Tsiganis, K., Nesvorný, D., & Gomes, R. 2011, AJ, 142, 152

Love, S. G., & Ahrens, T. J. 1996, Icar, 124, 141

Masiero, J. R., Wright, E. L., & Mainzer, A. K. 2021, PSJ, 2, 32

Morbidelli, A. 2005, arXiv:astro-ph/0512256

Morbidelli, A., Levison, H. F., Tsiganis, K., & Gomes, R. 2005, Natur, 435, 462

Nesvorný, D. 2011, ApJL, 742, L22

Nesvorný, D. 2018, ARA&A, 56, 137

Nesvorný, D., & Morbidelli, A. 2012, AJ, 144, 117

Nesvorný, D., & Vokrouhlický, D. 2016, ApJ, 825, 94

Nesvorný, D., Vokrouhlický, D., Bottke, W. F., & Levison, H. F. 2018, NatAs, 2, 878

Noll, K. S., Grundy, W. M., Chiang, E. I., Margot, J. L., & Kern, S. D. 2008, Binaries in the Kuiper Belt (Tucson, AZ: Univ. Arizona Press), 345

Pan, M., & Sari, R. 2005, Icar, 173, 342

Parker, A. H., & Kavelaars, J. J. 2010, ApJL, 722, L204

Petit, J.-M., Kavelaars, J. J., Bannister, M. T., et al. 2021, OSSOS XXX: The Size Distribution of the Kuiper Belt, in press

Petit, J.-M., Kavelaars, J. J., Gladman, B. J., et al. 2011, AJ, 142, 131

Petit, J.-M., Kavelaars, J. J., Gladman, B. J., et al. 2017, AJ, 153, 236

Pitjeva, E. V., & Pitjev, N. P. 2018a, AstL, 44, 554

Pitjeva, E. V., & Pitjev, N. P. 2018b, CeMDA, 130, 57

Rucska, J. J., & Wadsley, J. W. 2021, MNRAS, 500, 520

Schäfer, U., Yang, C.-C., & Johansen, A. 2017, A&A, 597, A69

Schlichting, H. E., & Sari, R. 2011, ApJ, 728, 68

Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, ApJ, 822, 55

Singer, K. N., McKinnon, W. B., Gladman, B., et al. 2019, Sci, 363, 955

Spencer, J. R., Stern, S. A., Moore, J. M., et al. 2020, Sci, 367, aay3999

Stern, S. A. 1995, AJ, 110, 856

Stern, S. A. 1996, AJ, 112, 1203

Trujillo, C. A., Jewitt, D. C., & Luu, J. X. 2001, AJ, 122, 457

Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Natur, 435, 459

Vilenius, E., Kiss, C., Müller, T., et al. 2014, A&A, 564, A35

Volk, K., Murray-Clay, R., Gladman, B., et al. 2016, AJ, 152, 23

Wetherill, G. W. 1967, JGR, 72, 2429
K.V. acknowledges funding from the NSF (grant AST-1824869) and NASA (grants 80NSSC19K0785 and 80NSSC21K0376).

## ORCID iDs

Abedin Y. Abedin https://orcid.org/0000-0002-7025-0975

J. J. Kavelaars https://orcid.org/0000-0001-7032-5255

Jean-Marc Petit https://orcid.org/0000-0003-0407-2266

Brett Gladman https://orcid.org/0000-0002-0283-2260

Michele Bannister https://orcid.org/0000-0003-3257-4490

Mike Alexandersen https://orcid.org/0000-0003-4143-8589

Stephen Gwyn https://orcid.org/0000-0001-8221-8406

Kathryn Volk https://orcid.org/0000-0001-8736-236X

## References

Abedin, A. Y., Kavelaars, J. J., Greenstreet, S., et al. 2021, AJ, 161, 195

Alexandersen, M., Gladman, B., Kavelaars, J. J., et al. 2016, AJ, 152, 111

Bannister, M. T., Gladman, B. J., Kavelaars, J. J., et al. 2018, ApJS, 236, 18

Bannister, M. T., Kavelaars, J. J., Petit, J.-M., et al. 2016, AJ, 152, 70

Benz, W., & Asphaug, E. 1999, Icar, 142, 5

Bernstein, G. M., Trilling, D. E., Allen, R. L., et al. 2004, AJ, 128, 1364

Chen, Y.-T., Gladman, B., Volk, K., et al. 2019, AJ, 158, 214

Clement, M. S., Kaib, N. A., Raymond, S. N., Chambers, J. E., & Walsh, K. J. 2019a, Icar, 321, 778

Clement, M. S., Kaib, N. A., Raymond, S. N., & Walsh, K. J. 2018, Icar, 311, 340

Clement, M. S., Raymond, S. N., & Kaib, N. A. 2019b, AJ, 157, 38

DAVIS, D. R., & Farinella, P. 1997, Icar, 125, 50

de Sousa, R. R., Morbidelli, A., Raymond, S. N., et al. 2020, Icar, 339, 113605

Dohnanyi, J. S. 1969, JGR, 74, 2531

Farinella, P., & Davis, D. R. 1996, Sci, 273, 938

Fraser, W. C., Bannister, M. T., Pike, R. E., et al. 2017, NatAs, 1, 0088

Fraser, W. C., Benecchi, S. D., Kavelaars, J. J., et al. 2021, PSJ, 2, 90

Fraser, W. C., Brown, M. E., Morbidelli, A., Parker, A., & Batygin, K. 2014, ApJ, 782, 100

Fraser, W. C., & Kavelaars, J. J. 2009, AJ, 137, 72

Gladman, B., Lawler, S. M., Petit, J.-M., et al. 2012, AJ, 144, 23

Gladman, B., Marsden, B. G., & Vanlaerhoven, C. 2008, Nomenclature in the Outer Solar System (Tucson, AZ: Univ. Arizona Press), 43

Gomes, R., Levison, H. F., Tsiganis, K., & Morbidelli, A. 2005, Natur, 435, 466

Greenstreet, S., Gladman, B., & McKinnon, W. B. 2015, Icar, 258, 267

Greenstreet, S., Gladman, B., & McKinnon, W. B. 2016, Icar, 274, 366

Greenstreet, S., Gladman, B., McKinnon, W. B., Kavelaars, J. J., & Singer, K. N. 2019, ApJL, 872, L5

Grundy, W. M., Noll, K. S., Nimmo, F., et al. 2011, Icar, 213, 678

Gwyn, S. D. J. 2012, AJ, 143, 38

Holsapple, K. A. 1993, AREPS, 21, 333

Holsapple, K. A., & Schmidt, R. M. 1987, JGR, 92, 6350

Housen, K. R., & Holsapple, K. A. 1990, Icar, 84, 226

Housen, K. R., Schmidt, R. M., & Holsapple, K. A. 1991, Icar, 94, 180

Housen, K. R., Sweet, W. J., & Holsapple, K. A. 2018, Icar, 300, 72

Jewitt, D., & Luu, J. 1993, Natur, 362, 730

Johansen, A., Jacquet, E., Cuzzi, J. N., Morbidelli, A., & Gounelle, M. 2015, New Paradigms for Asteroid Formation (Tucson, AZ: Univ. Arizona Press), 471

Jutzi, M., Michel, P., Benz, W., & Richardson, D. C. 2009, AAS/DPS Meeting Abstracts, 41, 50.06

Kavelaars, J. J., Jones, R. L., Gladman, B. J., et al. 2009, AJ, 137, 4917

Kavelaars, J. J., Lawler, S. M., Bannister, M. T., & Shankman, C. 2020, Perspectives on the Distribution of Orbits of Distant Trans-Neptunian Objects (Amsterdam: Elsevier), 61

Kavelaars, J. J., Petit, J.-M., Gladman, B., et al. 2021, ApJL, 920, L28

Kenyon, S. J., & Bromley, B. C. 2004, AJ, 128, 1916

Kenyon, S. J., & Bromley, B. C. 2017, ApJ, 839, 38

Kenyon, S. J., & Bromley, B. C. 2020, PSJ, 1, 40

Kenyon, S. J., & Luu, J. X. 1999, AJ, 118, 1101

Kovalenko, I. D., Doressoundiram, A., Lellouch, E., et al. 2017, A&A, 608, A19

Lawler, S. M., Kavelaars, J. J., Alexandersen, M., et al. 2018a, FrASS, 5, 14

Lawler, S. M., Shankman, C., Kavelaars, J. J., et al. 2018b, AJ, 155, 197

Leinhardt, Z. M., & Stewart, S. T. 2009, Icar, 199, 542

Leinhardt, Z. M., & Stewart, S. T. 2012, ApJ, 745, 79

Levison, H. F., Morbidelli, A., Tsiganis, K., Nesvorný, D., & Gomes, R. 2011, AJ, 142, 152

Love, S. G., & Ahrens, T. J. 1996, Icar, 124, 141

Masiero, J. R., Wright, E. L., & Mainzer, A. K. 2021, PSJ, 2, 32

Morbidelli, A. 2005, arXiv:astro-ph/0512256

Morbidelli, A., Levison, H. F., Tsiganis, K., & Gomes, R. 2005, Natur, 435, 462

Nesvorný, D. 2011, ApJL, 742, L22

Nesvorný, D. 2018, ARA&A, 56, 137

Nesvorný, D., & Morbidelli, A. 2012, AJ, 144, 117

Nesvorný, D., & Vokrouhlický, D. 2016, ApJ, 825, 94

Nesvorný, D., Vokrouhlický, D., Bottke, W. F., & Levison, H. F. 2018, NatAs, 2, 878

Noll, K. S., Grundy, W. M., Chiang, E. I., Margot, J. L., & Kern, S. D. 2008, Binaries in the Kuiper Belt (Tucson, AZ: Univ. Arizona Press), 345

Pan, M., & Sari, R. 2005, Icar, 173, 342

Parker, A. H., & Kavelaars, J. J. 2010, ApJL, 722, L204

Petit, J.-M., Kavelaars, J. J., Bannister, M. T., et al. 2021, OSSOS XXX: The Size Distribution of the Kuiper Belt, in press

Petit, J.-M., Kavelaars, J. J., Gladman, B. J., et al. 2011, AJ, 142, 131

Petit, J.-M., Kavelaars, J. J., Gladman, B. J., et al. 2017, AJ, 153, 236

Pitjeva, E. V., & Pitjev, N. P. 2018a, AstL, 44, 554

Pitjeva, E. V., & Pitjev, N. P. 2018b, CeMDA, 130, 57

Rucska, J. J., & Wadsley, J. W. 2021, MNRAS, 500, 520

Schäfer, U., Yang, C.-C., & Johansen, A. 2017, A&A, 597, A69

Schlichting, H. E., & Sari, R. 2011, ApJ, 728, 68

Simon, J. B., Armitage, P. J., Li, R., & Youdin, A. N. 2016, ApJ, 822, 55

Singer, K. N., McKinnon, W. B., Gladman, B., et al. 2019, Sci, 363, 955

Spencer, J. R., Stern, S. A., Moore, J. M., et al. 2020, Sci, 367, aay3999

Stern, S. A. 1995, AJ, 110, 856

Stern, S. A. 1996, AJ, 112, 1203

Trujillo, C. A., Jewitt, D. C., & Luu, J. X. 2001, AJ, 122, 457

Tsiganis, K., Gomes, R., Morbidelli, A., & Levison, H. F. 2005, Natur, 435, 459

Vilenius, E., Kiss, C., Müller, T., et al. 2014, A&A, 564, A35

Volk, K., Murray-Clay, R., Gladman, B., et al. 2016, AJ, 152, 23

Wetherill, G. W. 1967, JGR, 72, 2429