J J Kavelaars
email: jj.kavelaars@nrc-cnrc.gc.ca

Jean-Marc Petit

Brett Gladman

Michele T Bannister

Mike Alexandersen

Ying-Tung Chen

Stephen D J Gwyn

Typeset using

Keywords:

The cold main classical Kuiper Belt consists of those non-resonant small solar system bodies with low orbital inclinations and orbital semi-major axes between 42.4 and 47.

INTRODUCTION

Cold main classical Kuiper belt objects (KBOs) appear to be unevolved products of the initial planetesimal formation process in this region of the Solar System. The current number density of cold objects is such that collisions between these planetesimals are infrequent (e.g., Greenstreet et al. 2019;Abedin et al. 2021). The cold KBOs are known to contain a large number of loosely bound binary pairs [START_REF] Noll | Binaries in the Kuiper Belt[END_REF]; such pairs are very likely to be destroyed if collisions among KBOs are common (Petit & Mousis 2004) implying that the number density at the epoch of formation was similar to that we see today. The cold KBO pairs are so loosely bound that many would not have survived gravitational scattering into this zone of the solar system, implying they formed in-situ [START_REF] Parker | [END_REF]. Recently, observations of the cold classical KBO 486958 Arrokoth by the New Horizons mission have provided direct evidence of the low collision rate in this region (McKinnon et al. 2020). Arrokoth impactors are dominantly cold classical KBOs and the observed low crater density, well below the crater saturation threshold, is consistent with the ancient number density of material being within a factor of a few of that currently observed in this region (Greenstreet et al. 2019;Abedin et al. 2021). In addition, the photometric properties of the cold belt members appear distinct from the rest of the KBOs (e.g., Tegler et al. 2003;Pike et al. 2017;Schwamb et al. 2019).

The cold objects thus provide a window into the processes of planetesimal formation.

An examination today of the H r (absolute magnitude) distribution of cold objects larger than 20 km may provide a direct measurement of the distribution that emerged from the initial planetesimal formation processes. This population, unlike the collisionally evolved asteroid belt, never experienced the runaway growth to proto-or dwarf-planets, experienced minimal collisional erosion and has dynamical and surface properties that are distinct from the rest of the KBOs. The cold KBOs mass function today is the most likely to resemble the initial mass function of planetesimals.

The outcomes of planetesimal formation modeling are reaching a point where guidance from rigorous observational constraints are needed. The mass range in current model outputs overlaps with the well-observed range probed by the Outer Solar System Origins Survey (OSSOS). Here we present the high-fidelity measurement of the absolute magnitude distribution of the observed cold classical Kuiper belt as determined from an ensemble of survey samples [START_REF] Kavelaars | [END_REF]Petit et al. 2011;Alexandersen et al. 2016;Petit et al. 2017;Bannister et al. 2018) and associated detection characterizations (hereafter referred to as OSSOS++).

THE OSSOS++ COLD-BELT ABSOLUTE MAGNITUDE DISTRIBUTION

From the OSSOS++ sample (see Table 3 of Bannister et al. 2018, for full details) we select those 321 objects with free inclination i f ree < 4 • and semimajor axis in the range 42.4 au < a < 47.7 au to provide a relatively clean sample of the properties of cold main classical Kuiper belt objects. Van Laerhoven et al. (2019) found the orbital parameter that best separates the cold population from the background of the main classical Kuiper belt is the inclination with respect to the a-dependent Laplace plane, i.e., i f ree , and that i f ree < 4 • provides a reasonable split between the cold and excited populations. The 39.9 au a 42.4 au zone is usually included in the nominal definition of the main classical Kuiper belt. This zone, however, was likely destabilized by the passage of the ν 8 resonance during Neptune migration and KBOs on low-i orbits would have been removed. The few low-i f ree KBOs in this zone today are unlikely to have formed in situ and for that reason we exclude them when considering the cold classicals. Our sample criterion may exclude a small number of cold members from our analysis but ensures minimal contamination from other populations that may not have formed in situ and would distort the view of the unevolved H r distribution.

Sample characterization

Using the OSSOS++ characterization (Bannister et al. 2018) we debias the orbit and H r distributions of detected objects. The low eccentricity of the cold KBO orbits result in these objects exploring a limited range of solar distances within a constrained phase-space volume. As a result of this constraint, the OSSOS++ sample of cold objects provide a complete sampling of the orbit distribution (there are no hidden or unseen populations) and we can robustly debiased the detected sample. We consider the 4-dimensional phase space (a, q, sin (i f ree ), H r ) and slice it into cells small enough (0.2 au, 0.2 au, 0.001, 0.1) such that the distribution of elements within a cell are likely uniform (see Figure 5 of Bannister et al. 2018, to The black lines represent two exponentially tapered functions matched (see Section 3) to the debiased OSSOS++ data, with forced large-H r (small object) asymptotic slopes (dotted: α = 0.5; solid: α = 0.4). For H r < 9 the two model curves are nearly identical. The debiased OSSOS++ measurements are well matched by the exponential taper form.

The boxes represent literature derived estimates, see Table 1 andSection of the classical region. Although we computed bias factors for objects with larger H r values we do not use those in our analysis as the bias correction factors grow rapidly near the limit of detection.

For each of the observed cold objects in the sample, we determine the element cell, k, the detection belongs to and add 1/B k objects to our model, with the specific elements of those model objects drawn randomly across the cell's element distribution. Using this procedure provides an estimate of the number of objects in the cold population required to generate the OSSOS++ detections and the distribution of the population over each of the orbit and absolute magnitude (a, q, sin (i f ree ), H r ) parameters. Figure 1 presents the resulting OSSOS++ cumulative H r distribution.

2.2. Other observational constraints.

We compare with other cold KBO samples to verify our total population estimates at both the faint and bright ends of the H r distribution. Determination of population statistics requires carefully tracking each detected object to ensure that the orbit is accurately determined and object correctly classified. This tracking requires significant investment in telescope time if one is to avoid ephemeris bias entering the sample (Jones et al. 2006;[START_REF] Kavelaars | The Orbital and Spatial Distribution of the Kuiper Belt[END_REF]. This also places significant constraint on the faintness of objects that are allowed into a particular survey as the cost of tracking will rapidly increase. The OSSOS++ sample provides a high-quality sampling of the cold Kuiper belt due to the near 100% effectiveness in tracking detections to obtain high-quality orbits. The desire to achieve this high success rate in tracking, however, also limited the flux range accessible to the survey at both the bright and faint ends of the H r distribution. Other surveys which do not have precise orbit and distance estimates for all their objects, however, can be used to estimate the cold classical H-magnitude distribution and extend the flux range explored. We find that these additional samples agree well with the results measured via the OSSOS++ sample, providing a verification of the absolute calibration of our study. For each survey, we examined the published detection efficiency curves and determined the limiting apparent magnitude M r (all surveys used here reported limits in r) up to which the efficiency of detection, η, is roughly constant (see Table 1). For each survey we:

Deep studies

• count the number n of cold objects (estimated

1 inclination ≤ 4 • ) brighter than M r
• estimate the actual number of objects present in the field of view of the survey and brighter than M r as n/η

• determine the fraction of the full cold population that is in the field of view of a survey at any given time, F , using the CFEPS orbit model (Petit et al. 2011).

• compute an estimate of the implied full population brighter than M r :

N (m < M r ) = n/(F η).
The 95% confidence range for N is computed from the 95% Poisson confidence range for n and listed as N -for the lower end of the confidence range and N + for the upper in Table 1. To add these density estimates to our H r distribution we must convert the observed apparent magnitude limit B04 provides precise distance estimates and even rough orbital elements and the characterization of the detection efficiency. Thus, for B04, we also debiased the 3 detections of that survey following the same procedure as for OSSOS++. For B04, the full debiasing results in a population estimate at the largest H r of the detected B04 objects that is compatible with the estimate using the process outlined in the preceding paragraph (Figure 1).

We discuss implications for the observed H r distribution in Sections 3.

Inventory of the Brightest Cold Classicals

The now nearly complete, and sparse, inventory of the lowest H (largest) cold classical KBOs provides a further constraint on the size distribution. We select from the Minor Planet Center (MPC) database KBOs consistent with the a and i f ree cuts given above and with good orbits (2 oppositions or more, 6 observations or more, and MPC orbit uncertainty parameter less than or equal to 6). The MPC database provides the visual absolute magnitude H which we convert into H r using H -H r = 0.19, the mean value for the cold objects from OSSOS++ present in the MPC database. We then determine the number of known cold classical KBOs with H r < 5.1 (3) and H r < 5.5 (11) and include this as an estimate of the cumulative H r distribution; see Figure 1. These MPC-derived population estimates are nearly identical to the OSSOS++ based estimate of the total numbers at these H r values, confirming the global population estimates obtained by our debiasing.

That the OSSOS++ estimate of the total population of cold objects is the same size as the currently known sample, suggests that, as previously noted (e.g. Sheppard et al. 2011), the MPC database has reached (near) completeness for H around 5-6 mag. The expectation that the brightest members of the cold population would be in this range was also noted in Figure 7 of Bernstein et al. (2004) and is apparent in Figure 9 the prediction from the OSSOS++ H r -distribution, appears to be complete, and is small in number.

3. AN EXPONENTIALLY TAPERED H R DISTRIBUTION.

The shape of the cold population H r distribution presented in Figure 1 is inconsistent with twocomponent power-law2 fits. Using our de-biased model we imposed various two-component powerlaws (see Fraser et al. 2014) onto our model and compared the resulting orbit and H r distributions with our observed sample. No acceptable models (rejected by AD test statistic at more 99% confidence, see Lawler et al. 2018a, for details of our statistical procedure) were found. Moving to a multi-component power-law could provide a solution, but is difficult to physically motivate. The shape is also not consistent with the rolling power-law utilized in Bernstein et al. (2004). The steep and continuously evolving shape at the bright-end of the H r distribution is inconsistent with these representations.

A tapered power-law form appears to be emerging as a preferred functional form for planetesimal mass distributions. Disk instability mechanisms, and the streaming instability (SI) process in particular, have recently become highly favored solutions to overcoming various physical barriers in planetesimal formation (e.g., [START_REF] Safronov | Evolution of the protoplanetary cloud and formation of the earth and planets[END_REF]Johansen et al. 2007; Lyra & Umurhan 2019) and enable planetesimal formation to proceed more rapidly at lower surface densities, like those in the primordial cold Kuiper belt. Independent groups have investigated the initial mass function of planetesimals resulting from the SI (Johansen et al. 2015;Simon et al. 2016;[START_REF] Schäfer | [END_REF]Abod et al. 2019;Li et al. 2019;Liu et al. 2020;Rucska & Wadsley 2021). Due to limitation in the simulated mass resolution, the low-mass end of the mass function that emerges from SI models is not strongly constrained. The existing works, however, exhibit similar-shaped mass distributions that a power-law can roughly approximate at the small mass end but require a rather sharp, exponential, cut-off at the large mass end. Several functional forms to fit the mass distributions have been proposed. We select the exponentially tapered [START_REF] Schäfer | [END_REF] or variably tapered, Li et al. ( 2019)) power-law form.

According to Bayes criterion these forms provide better matches to the SI outcomes than simpler, single parameter functions [START_REF] Schäfer | [END_REF]Li et al. 2019).

The exponentially tapered power-law mass distribution [START_REF] Schäfer | [END_REF]) can be transformed into an H r distribution assuming a constant albedo and density, and spherical shape:

N (< H r ) = 10 3 5 α SI (Hr-Ho) exp 10 -3 5 β SI (Hr-H B ) (1)
where N will be the total population number and H o is a normalization, α SI is the asymptotic slope at large H r , β SI is the strength of the exponential tapering and H B is the H r value at which the exponential taper begins to dominate. The α in the single exponential H r distribution, N (< H r ) = 10 α(Hr-Ho) , is given by α = 3 5 α SI and can be seen as related to the faint/small object exponent when considering multi-component exponential distributions.

The OSSOS++ sample confirms the general shape of an exponential taper and can be used to determine the strength of that tapering (β SI ). The OSSOS++ sample is, however, limited to the large H r end where our debiasing factors are small; the value of α SI is not robustly constrained by the OSSOS++ sample. Crater counts on Pluto and Charon (Singer et al. 2019a), the observed sizes of Jupiter-family comets (Solontoi et al. 2012) and results from deep surveys shown in Figure 1 indicate that α 0.3-0.5 faint-ward of H r ∼ 9. In Figure 1 we present fits of Equation 1 with fixed values of α ∈ {0.4, 0.5} (α SI ∈ {0.66, 0.83}), with the value of β SI and the other free parameters determined using maximum likelihood MCMC parameter exploration via the emcee package (Foreman-Mackey et al. 2013) over the range H r ∈ {5.0, 8.3} The OSSOS++ estimates of β SI ∈ {0.42 +0.12 -0.16 , 0.59 +0.13 -0.27 }, H o ∈ {-2.6 +0.4 -0.9 , 0.0 +0.2 -0.4 } and H B ∈ {8.1 +1.7 -0.6 , 7.1 +0.9 -0.4 } provide a remarkably smooth match to our debiased observations and our estimates of the value of β SI are steeper than, but consistent with, the range of values found in SI modeling (0.28-0.37; e.g. [START_REF] Schäfer | [END_REF].

3.1. Comparison with published distributions.

Figure 2 presents our measured debiased H r distribution along with a number of results from the literature. The double or broken exponential forms presented in Figure 2 do not provide as compelling a match to the observations as the tapered exponential. There is good agreement that asymptotic small object exponential must be around α ∼ 0.4 but the forms do not provide a good match to the data at small H r which exhibits a continuously steepening slope. The OSSOS++ data is consistent with the bright end exponential slopes over some limited range of H r then drops away from those curves for smaller H r and larger H r . The bright/large object end of the H r distribution is not single sloped value and thus double exponential functions are not a good match.

DISCUSSION

The H r distribution of the OSSOS++ sample clearly demonstrates an exponentially tapered shape.

OSSOS was designed to be an absolutely calibrated survey that could be debiased to measure intrinsic absolute distributions inside each dynamical group, with the group membership based on highprecision orbits. These orbits permit the computation of free inclinations and the exclusion of resonant objects enabling the selection of relatively uncontaminated sample of cold belt members. We find the debiased H r distribution of this sample is inconsistent with a single power-law at the bright end (see Figures 1 and2). The functional form of the cold component H r distribution is well represented by an exponential taper of the type seen in numerical simulations of streaming instability driven planetesimal formation.

The OSSOS++ derived H r distribution matches well onto independent constraints at both ends of the distribution. At the faint (large H r ) end the debiased OSSOS++ sample connects smoothly to the faintest pencil-beam studies. Although each individual study is difficult to map to a precise N (< H r ) value, the ensemble of deep apparent-magnitude studies, collectively, match OSSOS++, to reflect difference in survey filters and for differences in selection function for cold classical KBOs. In particular we use (r -R) = 0.25 (Jordi et al. 2006) for (V -R) = 0.6 cold classical KBOs and we scale the apparent magnitude distribution given in B04 using a fixed distance of 42 au to transform from r to H r . The A14 total population is slightly low compared to OSSOS++ sample, this may be due to tracking losses reported in A14. The F14 fit has been scaled to match the OSSOS++ sample at H r = 8 as we were uncertain of the scaling from the surface density reported in F14 and the absolute total numbers reported here.

implying that the asymptotic form of a single exponential from H r > 9 down to at least H r ≈ 12 does not violate known constraints. The cratering records on Pluto and Charon (Robbins et al. 2017;Singer et al. 2019a) and Arrokoth (Spencer et al. 2020), however, indicate at least one further transition to an even shallower exponent beyond H r > 17. The OSSOS++ sample provides strong constraints in the H r 5-8.3 range; for H r > 8.3, where our detected sample drops off, a constant slope in log-space appears plausible for several more magnitudes. We highlight that within the hot component there is evidence for a knee or divot (Shankman et al. 2013(Shankman et al. , 2016;;Alexandersen et al. 2016) feature near H r =8.5 with a similar shallow slope for H r > 9. On the large object end, the cumulative distribution of the intrinsically brightest known cold-classical KBOs fall directly on the OSSOS++ curve and its bright-ward extrapolation. This close match confirms the very steep nature of the H-magnitude distribution of the largest objects, and implies that the inventory of these largest objects is essentially complete. These connections to independent constraints at H r =5 and 9 (with no tuning) gives confidence that in OSSOS++ we have an absolutely calibrated survey and that the H-magnitude distribution's shape between these two ends is correctly represented by our debiased measurement.

The existence of the exponential taper also resolves some literature confusion regarding the measured exponent of the KBO H r distribution. Much of the historical literature fit a single exponential to the apparent magnitude N (< M ) ∝ 10 αM . Because of the finite sky area and the small number of objects detected in any given survey, even fitting a simple exponential was challenging due to a lack of dynamic range. Generally, larger-area surveys were shallower, while deeper surveys made up for their smaller area via the steep H r distribution to end up with comparable (but small) samples. Magnitude distribution shape estimates done in apparent magnitude space (because they lack the precise distance estimates required to translate apparent to absolute magnitude) 'blur' the H-magnitude distribution, resulting in slope estimates that are a function of the depth of the survey; a similar conclusion was reached in Fraser et al. (2014).

Papers that attempted to combine surveys to enlarge the apparent magnitude and dynamic range of observational constraints tended to average this out to an intermediate slope (e.g. Gladman et al. 1998;Fraser et al. 2008). With an even larger apparent magnitude range it became clear that a single exponential could not represent the data, and double/rolling exponent forms were introduced.

Magnitude distribution studies near the solar system's invariable plane (which are thus dominated by low-inclination cold-component objects) exhibited a change of α to shallower values fainter than m r 24-25 (thus H r ∼8-9 in the main belt), when using a rolling (Bernstein et al. 2004), double (Fuentes et al. 2009) or broken power-law (Fraser & Kavelaars 2009) size distributions.

Although computationally convenient, there is no physical motivation for a broken power-law being the correct functional form. Numerical simulations of the streaming instability, however, appear to naturally produce the exponentially tapered form [START_REF] Schäfer | [END_REF]Li et al. 2019). The OSSOS++ sample demonstrates that this functional form is an excellent representation of the cold-component H-magnitude distribution. Much of the discussion of what value of α best matches the actual H r distribution and where to put a break to attempt to mimic the slope evolution, appears to be due to trying to model an exponential taper by combining multiple exponential distributions.

Having such examples of the danger of over-imposing a functional form on reality, we note that although this tapered exponential is clearly impressively similar to the distribution derived from the OSSOS++ sample, numerical simulations (e.g. Li et al. 2019) show that there can be smaller features superposed on this dominant form which, in differential space, manifest themselves as local 'knees'

(broken or double exponents) or 'divots' (Shankman et al. 2013). In particular, the weak relative under-abundance just past H r 7 present in previous data sets (e.g. Adams et al. 2014;Fraser et al. 2014) may be real, in addition to a proposed knee near H r 8.4 in the dynamically hot populations (e.g. Alexandersen et al. 2016;Lawler et al. 2018b).

The SI modeling prediction of the shape of the H-magnitude distribution at small sizes (large

H)
is not yet firmly established. The value of the asymptotic power-law slope seen at small sizes in simulations (Simon et al. 2017(Simon et al. , 2016;;Abod et al. 2019;Rucska & Wadsley 2021) may result from resolution effects. Li et al. (2019) demonstrate that as resolution is improved, what initially appeared to be a roll-over to a 'single-α' asymptotic form, continues to evolve and the transition to the asymptotic form appears to occur at ever smaller sizes. Interestingly, the existing estimates of these asymptotic limits (see Figure 10 in Abod et al. 2019, for example) are not very far from the often-suggested value of α 0.4 ± 0.1 for H r 8 (Bernstein et al. 2004;Fuentes & Holman 2008;Fraser et al. 2010Fraser et al. , 2014)). This allows the possibility that this value of α is then set by the formation size distribution and the similarity of the observed H r distribution slope to the collisional equilibrium value is not evidence of collisional equilibrium having been achieved. The small object size distribution of the cold classical Kuiper Belt appears to be unaltered over the age of the Solar System and to preserve a shape consistent with SI planetesimal formation down to of order kilometer scale.

The comparison with current modeling of SI driven planetesimal formation is not, however, without a significant hurdle. Although the shape of the distribution is compelling, the inferred mass ranges are not a good match. If the current cold component is indeed a relatively unevolved population, then the current surface density may be taken as a proxy of the density at formation. Based on our estimate of the total cold population (Figure 1) we can infer the surface density at the time of the SI process. To estimate this density we convert our H r distribution to a mass distribution by assuming a constant albedo of 0.15 and object bulk density of 500 kg/m 3 (Spencer et al. 2020) and then spread the inferred total mass into a 2 au wide ring centered at 43 au. The resulting estimated primordial surface density in solids is Σ p ∼ 5 × 10 -5 g/cm 2 . This is 100 times lower than the mass scale required for most SI modeling to produce D ∼ 400 km objects (e.g., Abod et al. 2019). The number densities inferred from the current population are far too low to be consistent with SI processes forming the sizes of KBOs in our observed size distribution, which we claim follows the shape seen in SI modeling! Conversely the number density at the time of planetesimal formation implied by SI modeling results appears to be significantly higher than that inferred from the currently observed population.

There appear to be at least the following possible issues to consider:

• The current density could be many factors lower than at the time of planetesimal formation.

This appears unlikely as any process that removed significant mass would very likely have disrupted the binary KBO population we see today. The cratering record seen on Arrokoth is fully consistent with the low total populations reported here (see Greenstreet et al. 2019;Singer et al. 2019b;Abedin et al. 2021) indicating that any period of high number density would have been very rapidly removed which would have implications for the orbit distribution in the cold belt (Gladman & Volk 2021). Thus, it appears unlikely that the surface density at ∼ 43 au was significantly higher than today.

• The SI process is not directly responsible for the production of the objects we see today, but they formed instead after the SI process via particle-particle interactions. A 2 au wide ring would contain many hundreds to thousands of SI cells (e.g. Li et al. 2019), and perhaps the planetesimals from these cells coalesce to form the largest bodies. Here, again, the currently observed density makes this appear improbable as the inferred particle-particle encounter rate would be very low (thus the low numbers of craters) making the rate of planetesimal growth so slow that the largest sized objects would have not yet formed. In addition, there is no reason that the observed H r distribution of the post-growth populations would then so closely resemble that which emerges from the SI processes.

• The models of SI are incomplete and the process happens on scales and at densities that have not yet been fully modelled. This appears unlikely to be the case as the modeling is done in scale free units and then the mass scales are set by imposing a density.

• The classical Kuiper belt was much more tightly confined, radially and azimuthally, at the time of SI driven planetesimal formation and then rapidly dispersed shortly after the planetesimals emerged. This ad hoc solution has some appeal as many planet forming disks exhibit density enhancements (e.g. van der Marel et al. 2013) and the shearing out of particles that might form in such dense regions would be quite rapid compared to collisional time-scales. To achieve the required density enhancement, however, would require a concentration that is 100 times that seen today. Are such density enhancements feasible? Some caution is thus warranted because what is observed today in the cold-classical Kuiper belt is, in this paradigm, the end state of streaming instability plus later accretion and erosion. If the latter two processes are indeed negligible (or could be successfully modeled) then the shape of the cold classical belt H-magnitude distribution becomes a direct measure of the outcome of the planetesimal formation process from the streaming instability. Future numerical work can then in principle constrain the protoplanetary disk's parameters (surface density, viscosity, etc.) as one targets reproducing the observed distribution in the preserved cold classical belt.

Regardless of a possible link with the SI process, our analysis of the OSSOS++ sample has provided a robust high-fidelity measure of the H r distribution of the cold component of the Kuiper belt. This analysis is enabled by the precise characterization of the OSSOS++ surveys. The derived shape exhibits an exponentially tapered form and, from the current evidence, is representative of the initial distribution resulting from planetesimal formation.

Figure 1 .

 1 Figure 1. H r distribution of the cold main Kuiper belt. The grey-dash line represents the distribution of raw detections in the OSSOS++ sample. Red-orange curve (shown as a dotted line for H r > 8.3 where our sample debiasing is less secure) represents the debiased OSSOS++ sample with the shading indicating the Poisson 95% confidence range.

  Gladman et al. (2001);Bernstein et al. (2004);Fraser et al. (2008);Fuentes & Holman (2008);Fraser & Kavelaars (2009) (hereafter G01, B04, Fr08, FH08, and FK09) performed deep 'pencilbeam' surveys to detect faint TNOs and published their detection efficiency functions. The imprecise determination of heliocentric distance at detection in G01, Fr08, FH08 and FK09, however, prevents direct conversion from observed brightness to H r (absolute) magnitudes. We carefully examined the sample of detections from each of these projects. The VLT part of G01 did not yield any detections and we keep only the CFHT component of that project (G01/CFHT). The inclinations in the Blanco part of Fr08 are insufficiently constrained to allow selection between cold and excited objects in that sample, thus we utilized only the CFHT component from that project (Fr08/CFHT). The details provided in FH08 and FK09 are sufficient to allow use of their full sample of detected cold KBOs.

(

  M r ) to an absolute (H r ) value, which requires knowledge of the distance to the sources. As the precise distances of the individual detections are not known we determine a plausible range of H r values by adopting the 95% range of distances of the cold objects from the CFEPS model(Petit et al. 2011; Fraser et al. 2014) (40.3 to 51 au). Using this approach we determine H r (51) and H r (40.3) (using H 51 = M r -17.03 and H 40.3 = M r -16.00) as representative of the range H r values that each survey was sensitive to at the detection limit (M r ).

  of Fraser et al. (2014). The expectation of completeness in the MPC sample is also coherent with the reported detection of large cold objects by Pan-STARRS, which surveys the whole ecliptic. During the period 2010-2014, Pan-STARRS found 7 of the 20 largest objects in the cold belt region defined above. Since 2014 no new large cold objects have been reported by Pan-STARRS, despite continuous operation. An absence of H r < 4 cold objects was also predicted a decade ago based on the CFEPS project (seeSec. 5.1.1 of Petit et al. 2011) whose sample is included in OSSOS++. The total inventory of known cold population objects with H r < 5.5 matches precisely

Figure 2 .

 2 Figure 2. As in Figure 1 the red-orange region represents the debiased OSSOS++ sample. The green shaded area represents the debiased detections from the Deep Ecliptic Survey (DES) (Adams et al. 2014) with the green dashed line their best-fit double exponential. Also shown are the best fits from Bernstein et al. (2004, cyan dotted line) and Fraser et al. (2014, magenta dashed line). The curves have been scaled

Table 1 .

 1 Deep Surveys. If a telescope is listed, the sample is restricted to that particular portion of the study.

	Survey a	η	M r n cold	1/F	H 51 H 40.3	N -	N +
	Fuentes & Holman (2008)	0.88 25.1	30	930 8.1	9.1 22300	45300
	Fraser et al. (2008) CFHT	0.97 25.2	11	1726 8.2	9.2 11000	35000
	Gladman et al. (2001) CFHT 1.00 25.8	3 11700 8.8	9.8 10733	94100
	Fraser & Kavelaars (2009)	0.95 26.5	10 11594 9.5 10.5 67000 224500
	Bernstein et al. (2004)	1.00 28.5	3 115325 11.5 12.5 125700 1011100
	a

We cannot compute i f ree as the semimajor axis of the orbits of the objects are, generally, unknown.

We sometime use the term power-law to refer to the exponential forms such as N (< H) ∝ 10 αH as the underlying mass distribution is, in-fact, a power-law form and H is used as a proxy for that quantity, and referring to exponentially tapered exponential functions becomes exponentially confusing, but does reflect the exponential complexity of reality.

ACKNOWLEDGMENTS JJK and BG acknowledge support of the Natural Sciences and Engineering Research Council of Canada. This work was supported by the Programme National de Planétologie (PNP) of CNRS-INSU co-funded by CNES. This research made use of the Canadian Advanced Network for Astronomy Research (CANFAR) and the facilities of the Canadian Astronomy Data Centre operated by the National Research Council of Canada with the support of the Canadian Space Agency. KV acknowledges support from NASA (grants NNX15AH59G and 80NSSC19K0785) and NSF (grant AST-1824869). Based on observations obtained with MegaPrime/MegaCam, a joint project of CFHT and CEA/DAPNIA, at the Canada-France-Hawaii Telescope (CFHT) which is operated by the National Research Council (NRC) of Canada, the Institut National des Sciences de l'Univers of the Centre National de la Recherche Scientifique (CNRS) of France, and the University of Hawaii. The observations at the CFHT were performed with care and respect from the summit of Maunakea which is a significant cultural and historic site.


Abedin Y Abedin
email: abedin.abedin@nrc-cnrc.gc.ca

J J Kavelaars

Sarah Greenstreet

Jean-Marc Petit

Brett Gladman

Samantha Lawler

Michele Bannister

Mike Alexandersen

Ying-Tung Chen

Stephen Gwyn

Typeset using L A T E X default style in AASTeX63 Collision Probabilities in the Edgeworth-Kuiper belt

Keywords: KBOs, Small Solar System bodies, Collisions, Probabilities

Here, we present results on the intrinsic collision probabilities, P I , and range of collision speeds, V I , as a function of the heliocentric distance, r, in the trans-Neptunian region. The collision speed is one of the parameters, that serves as a proxy to a collisional outcome e.g., complete disruption and scattering of fragments, or formation of crater, where both processes are directly related to the impact energy. We utilize an improved and de-biased model of the trans-Neptunian object (TNO) region from the "Outer Solar System Origins Survey" (OSSOS). It provides a well-defined orbital distribution model of TNOs, based on multiple opposition observations of more than 1000 bodies. In this work we compute collisional probabilities for the OSSOS models of the main classical, resonant, detached+outer and scattering TNO populations. The intrinsic collision probabilities and collision speeds are computed using the Öpik's approach, as revised and modified by Wetherill for non-circular and inclined orbits. The calculations are carried out for each of the dynamical TNO groups, allowing for inter-population collisions as well as collisions within each TNO population, resulting in 28 combinations in total. Our results indicate that collisions in the trans-Neptunian region are possible over a wide range in (r, V I ) phase space. Although collisions are calculated to happen within r ∼ 20-200 AU and V I ∼ 0.1 km/s to as high as V I ∼ 9 km/s, most of the collisions are likely to happen at low relative velocities V I < 1 km/s and are dominated by the main classical belt.

INTRODUCTION

The Edgeworth-Kuiper belt (Kuiper belt hereafter) is a torus-shaped agglomeration of small icy objects, encircling the Solar System. Extending from the orbit of Neptune at 30 AU to more than 200 AU from the Sun, the Kuiper belt consists of more than ∼ 10 5 trans-Neptunian objects (TNOs) with size D > 100 km [START_REF] Petit | [END_REF]. TNOs can be considered to consist of pristine building blocks, leftover from Solar System formation; they are unlikely to have been exposed to temperatures of T > 50 K, (e.g. Jewitt & Luu 2004) over the age of the solar system.

Collisions between bodies inevitably lead to ejecta, regardless of whether the impact is disruptive or crater forming, with the total mass and velocity distribution of liberated material being dictated by the energy of the impact (e.g. [START_REF] Melosh | Impact cratering : a geologic process[END_REF][START_REF] Fujiwara | Asteroids II[END_REF][START_REF] Housen | [END_REF]Benz & Asphaug 1999;[START_REF] Holsapple | Asteroid Impacts: Laboratory Experiments and Scaling Laws[END_REF]. The ejecta size distribution generally follows some power-law, ranging from large boulders to dust particles, with the number of released fragments being inversely proportional to particle size. Understanding the impact velocity distribution and probability will aid in understanding the impact environment of these bodies and the expected rate of dust/debris production. These impacts may very well be responsible for formation of dust rings and satellites around planets [START_REF] Harris | IAU Colloq. 75: Planetary Rings[END_REF][START_REF] Charnoz | [END_REF][START_REF] Charnoz | The Origin of Planetary Ring Systems[END_REF], dusty exospheres around small bodies (e.g., Szalay et al. 2018), serve as source and replenishment of the interplanetary dust flux (e.g. Kuchner & Stark 2010;Vitense et al. 2012;[START_REF] Poppe | [END_REF] or lead to planetary surface cratering (e.g. Greenstreet et al. 2015Greenstreet et al. , 2019)). In particular, dust is dynamically short lived and collisions provide an avenue to refresh that content.

Owing to the fundamental importance of the collisional processing for understanding the Kuiper belt and the outer solar system, a number of very early estimates of collisional environment were made with only the slimmest of observational data to ground them. Stern (1995) argued that TNOs smaller than ∼ 10 km are collisionally evolved, based on collision rate estimation between 35-70 AU, while objects with radii r > 100km may retain their primordial sizes. However, the only available observational statistics at the time was 25 objects of size similar to 1992 QB 1 and the rate of observed Jupiter-Family Comets (JFCs), assuming they originate from the low inclination component of the present Kuiper belt. Davis & Farinella (1997) argued that KBOs with diameters (D < 20 km), are significantly collisionally evolved, producing on average 10 objects per year of size 1-10 kilometers. Durda & Stern (2000) expanded on previous works using an updated number of discovered KBOs and estimated that 2.5 km KBOs have a collisional life-time of 3.5 Gyr, while sub-kilometer objects collide with 100 km TNOs once every 10 7 -10 8 years. In the intervening decades since these early works a much improved picture of the orbital structure and size distribution of the Kuiper belt has emerged.

Surveys for Kuiper belt objects (e.g. Jewitt & Luu 1993, 1995;[START_REF] Gladman | [END_REF]Chiang & Brown 1999;Gladman et al. 2001;Trujillo et al. 2001;[START_REF] Millis | [END_REF]Jones et al. 2006;Kavelaars et al. 2009;Schwamb et al. 2010;[START_REF] Petit | [END_REF]Brown et al. 2015;Alexandersen et al. 2016;Bannister et al. 2016Bannister et al. , 2018) ) have demonstrated that the Kuiper belt contains a rich diversity of orbital sub-components. Contemporary models of the formation of the Kuiper belt imply that Neptune has played a significant role in placing and reshaping the orbits of these primordial bodies (Malhotra 1993(Malhotra , 1995;;Levison & Morbidelli 2003;Gomes et al. 2005;Tsiganis et al. 2005;Nesvorný 2015;Shannon et al. 2019) with the orbital evolution of Neptune "fossilized" into the current state of the Kuiper belt. Furthermore, studies of size-frequency distributions (SFDs) of the TNO sub-populations (e.g. Jewitt & Luu 1995;[START_REF] Gladman | [END_REF]Gladman et al. , 2001;;Trujillo et al. 2001;Bernstein et al. 2004;Fraser & Kavelaars 2008;Fraser et al. 2008;[START_REF] Petit | [END_REF]Fraser et al. 2014;Lawler et al. 2018b) provide clues to the agglomeration and collision history. Studied as a whole, the Kuiper belt provides many insights into the processes of planet formation and planetary system evolution.

With the increased number of discovered TNOs more complete Kuiper belt dynamical and physical models have become possible allowing more precise determination of collision probabilities. Using the CFEPS L7 model [START_REF] Petit | [END_REF]), Dell'Oro et al. (2013) calculated the collision probabilities between different TNO populations and subpopulations, finding intrinsic collision probabilities to be lower than previous estimates. Using refined versions of the L7 model, and informed by various luminosity function studies, Greenstreet et al. (2015Greenstreet et al. ( , 2016) ) made predictions for the collision probabilities and cratering rate onto Pluto and Charon prior to the NASA New Horizons Pluto flyby and later made predictions for New Horizons Kuiper belt Extended Mission target "Arrokoth" (2014 MU 69 ) (Greenstreet et al. 2019).

There are over 1200 KBOs in the OSSOS ensamble sample (Bannister et al. 2018), with high-precision orbital elements and detailed discovery characterization. This sample allows for the construction of a 'de-biased' TNO model, using the OSSOS survey simulator [START_REF] Petit | [END_REF]Lawler et al. 2018a). This improved knowledge of the structure and history the Kuiper belt enables a yet more precise determination of the collisional environment.

We present here our determination of the intrinsic collision probabilities between the orbital sub-populations of the Kuiper belt based on this debiased model. The intrinsic probability is a function of the orbital distribution of objects. The orbital distribution of different dynamical TNO classes is the foundation for computing how frequently these objects collide. When combined with a TNO size frequency distribution model the intrinsic probabilities determine a collision rate. The goal of this work is to use the most up-to-date model of the Kuiper belt orbit distribution to obtain refined collision probabilities and speeds.

In Section 2 we describe the orbital and size distributions used in our probability calculations. Section 3 presents our approach in determining the intrinsic and total collision probabilities in the current Kuiper belt. In Section 4 we present the intrinsic collision probabilities, considering each KBO sub-population in-turn and then combining those into a total collisional probability for objects in the Kuiper belt. In Section 5 we discuss some zeroth-order implications determined from our collisional probabilities.

DYNAMICAL STRUCTURE OF THE KUIPER BELT

Determination of orbital elements for TNOs depends on accurate and multi-opposition astrometry, which is often challenging due to the large distances and faintness of TNOs (e.g. Jones et al. 2006). The apparent daily motion of a TNO, at a distance of r ∼ 40 au, is small (mean motion n ∼ 6.8 × 10 -5 au/day) and hence the observed arc (θ ∼ 0.2 arcmin), even over a month, will be a tiny fraction of the orbit. Often, observations over the course of a number of years are needed before the precise orbit of the TNO can be determined (Bannister et al. 2018). Another difficulty is introduced from the r -4 distance dependence of the reflected light from a KBO, registered by telescopes. Contemporary KBO detection is limited to objects with R-band magnitude brighter than m R 26.5 (Petit et al. 2006;Fraser et al. 2008) (diameter roughly D > 30 km), although there have been reports of detection of smaller D ∼ 1 km KBOs using stellar occultation (Schlichting et al. 2009;Arimatsu et al. 2019). Furthermore, the rate of detection per unit area is a strong function of pointing direction: TNOs moving on highly inclined orbits have a better chance of being discovered away from the ecliptic while those on low-inclination orbits are only detectable when the observed field is near the plane of the TNO's orbit. Similarly, TNOs on eccentric orbits are more difficult to detect, brightening sufficiently only near their perihelia, but spending most of their time near aphelia. Conversely, objects on more circular, moderate inclination orbits are more readily detected. Therefore a dedicated TNO survey is more suitable and desired for extensive follow-up observations, as well as characterization of the above biases. For a review on the biases associated with KBO observations, see [START_REF] Kavelaars | Perspectives on the distribution of orbits of distant Trans-Neptunian objects[END_REF].

Nevertheless, since the discovery of the first KBO (Jewitt & Luu 1993), there have been a number of KBO survey programs, among them Trujillo et al. (2001), Schwamb et al. (2010), Deep Ecliptic Survey [START_REF] Millis | [END_REF], Canada-France Eclipitc Plane Survey (CFEPS) (Jones et al. 2006;Kavelaars et al. 2009;[START_REF] Petit | [END_REF][START_REF] Gladman | [END_REF], Alexandersen et al. (2016) (MA) and the Outer Solar System Origin Survey (OSSOS) (Bannister et al. 2016(Bannister et al. , 2018;;Lawler et al. 2018b). Presently, there are > ∼ 3300 cataloged TNOs in the Minor Planet Center (MPC) https: //minorplanetcenter.net/iau/mpc.html, but only half of them have sufficiently well-constrained orbits to allow for precise long term numerical integration of their orbits and accurate dynamical classification, with the majority of those coming from the OSSOS Ensemble sample which combines the CFEPS, MA and OSSOS samples. Therefore, we choose to utilize the most up-to-date bias-free TNO dynamical model by OSSOS as constrained by the OSSOS Ensemble of over 1000 TNO orbits (Bannister et al. 2018).

For this work we consider the collisional circumstances for different orbit classes of TNOs. We group the TNOs into dynamical groups proposed by [START_REF] Gladman | Nomenclature in the Outer Solar System[END_REF], based on long-term numerical orbital integration. 1. Resonant objects -TNOs that are currently in mean motion resonance (MMR) with Neptune.

2. Scattering objects -TNOs that experience close encounters and "scattering" off Neptune over 10 Myr of numerical orbital integration, 3. Classical and Detached objects -TNOs which do not fall into either of the above categories and is often further sub-divided on semi-major axis and inclination.

(a) Inner belt -classical objects which are interior to the 3:2 MMR with Neptune. (b) Main belt -The main classical belt KBOs, which have semi-major axes between the 3:2 and 2:1 MMR with Neptune. The main belt contains at least two distinct sub-classes: "Hot" and the "Cold" components, where the Cold appears to be constrained to a narrow ring between 42.5 au and 44.5 au and exhibits complex substructure. The Hot component is distinguished by a wide inclination (i) distribution whereas the Cold component's inclinations are tightly confined. (c) Outer belt -Objects with semi-major axes beyond the 2:1 MMR and eccentricities e < 0.24. (d) Detached belt -Objects with semi-major axes beyond the 2:1 MMR and eccentricities e > 0.24.

To compute the collision probabilities in the Kuiper belt we require models that describe the orbital distributions of each of these Kuiper belt classes.

Resonant class

As we are interested in the global collision rates we only consider those resonances with a significant number of known members and a reliable orbit distribution model available in the literature. For the sake of simplifying the computational problem, without significant loss in collisional probability precision we compute the probabilities grouped by semi-major axis:

1. Inner resonances, in order of increasing semi-major axis, N4:3, N3:2 2. Main resonances, (semi-major axis within the main belt resonances), N5:3 and N7:4

3. Outer resonances N2:1, N7:3 and N5:2 Where 'N' stands for a MMR with Neptune and the convention used is p : q, where p and q are the small integers that refer to the number of orbits completed by the inner body (Neptune in our case) and the outer body (TNO) respectively, where p -q is referred to as the "order of resonance". Our collisional modelling does not account for the full resonant dynamics of Kozai-type orbital oscillation nor the confined precision of the argument of pericenter as we treat these objects as having uniformly precessing angular elements. However, Greenstreet et al. (2015) find that accounting for the resonant and Kozai dynamics of the plutinos lead to an increased collision probabilities and rates onto Pluto by ∼50%, though as will be shown later, neglecting the non-uniform precession of the orbital elements of resonant TNOs will not change our overall conclusion.

Scattering

The orbital distribution of the scattering component can be determined via comparison with cosmogonic models of the formation of the Kuiper belt and Oort Cloud. By definition, the orbital elements of the Scattering objects are 'actively' evolving. Shankman et al. (2016) and Lawler et al. (2018b) found that the currently observed orbital distributions of the scattering members of the OSSOS ensemble sample are well represented by the output of simulations in Kaib et al. (2011). We utilize the orbital distributions provided by Lawler et al. (2018b) as our orbit model for the Scattering component.

Classical Kuiper belt

Here, rather than consider parametric distributions of sub-components of the Classical Kuiper belt we utilize an orbit model derived by debiasing a detected sample. In this way the detected Classical belt objects and the discovery survey's efficiency of detecting such orbits are used to invert detected population of KBOs to determine the intrinsic population. For this work we utilized the debiased model described in [START_REF] Petit | OSSOS XXX: The size distribution of the Kuiper belt[END_REF].

Following [START_REF] Petit | OSSOS XXX: The size distribution of the Kuiper belt[END_REF] we split the main classical belt on the basis of orbital inclination; Cold-KBOs with i < 5 • and Hot-KBOs with i ≥ 5 • . Here, "Cold" and "Hot" refer to dynamical state of the orbits, and are not related to the objects' temperature. This somewhat crude splitting is motivated by the fact that, although there is an overlap between Cold and Hot objects in the eccentricity and inclination space, one can see a clear concentration of objects with low (i < 5 • ) inclination (e.g. Brown 2001;[START_REF] Petit | [END_REF]. This inclination cutoff is combined with low (e < 0.08) eccentricity threshold to select out those low-i KBOs that are more strongly interacting with Neptune and group them with the Hot objects for the purpose of considering average collision probabilities. We thus have two distinct TNO populations in the main classical belt, which allows us to test which of these ("Hot" vs. "Cold") are more likely to experience collisions in the Kuiper belt.

Outer/Detached belt

Significantly fewer non-resonant members of the Kuiper belt with a > 47 au, compared to a < 47 au, have been detected, making a direct debiasing statistically less certain. For this region we revert to the CFEPS L7 [START_REF] Petit | [END_REF]) parametric model. The L7 model represents the orbit distribution of this component as having e and i that follow the Hot Main Kuiper belt population but with a following a probability density function (PDF) ∝ a -5 2 and the remaining elements drawn randomly. The a PDF is normalized such that the Outer/Detached component represents a smooth continuation from the Hot Main Kuiper belt. As we will see, collisions among this population are rare and the precise orbital distribution is not likely to significantly modify conclusions about impact rates.

The seven main populations considered here -Cold and Hot Classicals, Inner Resonant, Main Resonant, Outer Resonant, Detached+Outer belt and Scattering disk we have a total of 28 orbit class combinations.

COLLISION PROBABILITIES

We compute the intrinsic collision probability P I and speeds V I between each pair of orbits within, as well as between different dynamical classes. In order to compute the intrinsic collision probability between each and every TNO dynamical class, we randomly select a sample of 10 4 orbits from orbital distributions in each of the sub-population. The intrinsic collision probability is independent of the sample size but larger samples provide somewhat better sampling of the available phase space for the orbit class being considered. For each orbit family pair we computed the intrinsic collision probability between each pair of orbits in each population: the number of calculations roughly scales as ∼ N 2 , where N is the number of individual orbits considered. Each probability computation required an average CPU time of t ∼ 1.5 ms per orbit pair. Performing tests with different sample sizes, we concluded that computation between 10 4 orbits, which is equivalent to ∼ 5 × 10 7 orbit pairs (requiring about 20 hours of computing) provided statistically equivalent behaviour to a 10 5 sample (requiring about 200 hours of computing). We then compute the intrinsic collision probabilities between each of the 28 permutation of our orbit classes.

We use the approach described in Wetherill (1967) to compute the intrinsic collision probability between two orbits:

P I = V I R h 8π 2 sin i | cot α|a 1 a 2 √ p 1 p 2 , p 1 = a 1 (1 -e 2 1 ) p 2 = a 2 (1 -e 2 2 ) R h = p 2 1 + e 2 cos ω cot α = R h ω/ Ṙh (1)
where p 1 and p 2 are the semi-latus rectums of the two orbits, i is the mutual inclination and V I is the relative velocity between the two bodies at the point of orbit intersection. The distance R h is the heliocentric distance of the intersection point (mutual node) of the two orbits and ω is the orientation of perihelion of the impacting orbit and has a meaning of true anomaly of the intersection point. α is the angle between object's rectilinear trajectory and the heliocentric distance vector. The intrinsic collision probability is a function of the orbital elements of the given orbit pair only, without involving their size. Due to gravitational perturbations, the arguments of pericenter and the longitudes of the ascending node of the two orbits will experience secular variations, causing the mutual inclination between the two orbits to process. Eq. 1 does not incorporate the longitude of the ascending node and the argument of pericenter, as the formalism of (Opik 1951) and (Wetherill 1967) assumes that all values for these parameters, during one procession cycle, are equally probable1 .

For the computations we make use of the FORTRAN implementation used by Greenstreet et al. (2015) to calculate the collision probabilities and speeds between each orbit pair within and between different TNO populations. However, we extended this code to provide the intrinsic collision probability as both a function of speed (as in Greenstreet et al. (2015)) and also the heliocentric distance at time of collision. We provide our results as a matrix of 200×135 distance, velocity cells. For each cell we compute the average collision probability within a 1 au, 1 m/s volume at a the heliocentric distance R with an impact speed V I for that cell. The full grid of collision probabilities is provided in electronic form.

The results of the calculated intrinsic collision probability P I and speeds V I , between and within different KBO populations, are summarized in Figures 1 2,3. Each pixel in each subplot is color coded in terms of the intrinsic collision probability P I in [km -2 yr -1 TNO -1 ] occurring with a collision speed V I in [km/s] at a given heliocentric distance R h in [AU]. Given the large range of variation of ∆P I (over seven orders of magnitude), the data in each plot is presented in a logarithmic scale. In addition, we also provide separately the probability density of a collision happening with a given impact speed V I , for the different dynamical TNO populations (Fig. 6). The intrinsic collision probabilities between each KBO population i.e., the average of all pixel values for each subplot are summarized in Table1. The intrinsic probabilities must be multiplied by the size of the populations to determine the collision rate. The probability plots provide a quick visualization of where in the Kuiper belt collisions are most likely to occur and which orbital classes have the largest intrinsic probability of collisions.

Our calculated intrinsic probabilities and speeds indicate that collisions in the Kuiper belt are possible over a wide range of the heliocentric distances R h = 30 -200 AU and speeds V I < 9 km/s, though the majority of these collisions are expected to occur in a narrow range between R h = 30 -55 AU and with much lower speeds (see Fig. 1, 2 and 3). Although, we see collision speeds as high as 9 km/s, there are very unlikely to occur over the age of the Solar System with speed of V I < 6 km/s being more likely (Fig. 6).

The highest intrinsic collision probability occurs within the "cold" sub-component of the classical belt, followed by the "hot" on "cold" impacts, and within the "hot" substructure. This is not surprising given that these populations are more concentrated in the orbit phase-space compared to the other orbit classes considered. It is also evident that the cold classicals collide with median velocity of V I ≈ 0.3 km/s and have a narrow impact speed range (see Fig. 6 (a)). That raises an interesting question as to whether the collisions between these objects are disruptive, crater forming or result in complex structure formation. The expectation of low impact velocities is similar to the result reported in Greenstreet et al. (2019) when considering impacts onto Arrokoth and utilizing the CFEPS L7 orbit model.

We also note a bi-modality in the impact speed distribution for "hot" on "cold" collisions (Fig 6a). The reason for this is the existence of substructures in the semi-major axis distribution of the "cold" and "hot" sub-components of the main belt creating distinct sub-sets of collisional circumstances. A bi-modal velocity distribution is also seen for the scattering population colliding with the "cold" sub-structure of the main belt (Fig. 6d). This structure also appears to be due to the sub-structure in the semi-major axis distribution of the "cold" population.

The intrinsic probabilities we report here are consistent (within an order of magnitude) with those reported in Dell'Oro et al. ( 2013) who performed similar analysis using the CFEPS-L7 model. The impact speeds reported in Dell'Oro et al. ( 2013), however, are somewhat (∼ 70%) higher compared to the values reported here and in Greenstreet et al. (2019) (which were also derived using the CFEPS L7 model).

Combined Intrinsic Probabilities

The P I for the Kuiper belt can now be determined by combining together the P I for the various sub-components determined in the previous section. To do this we must first weight each probability matrix by the relative size of the the components of the Kuiper belt population involved. Frustratingly, the actual relative sizes of these populations depends on the size range of sub-populations being considered as the 'Excited' members of the Kuiper belt (the Hot Classical Kuiper belt, Resonant objects, Scattering disk and Outer/Detached components) appear to have one size distribution while the size-frequency distribution (SFD) of the Cold Classical Kuiper belt appears to be different. In Table 2 we summarize current literature population estimates of the various orbit classes in the Kuiper belt, down to the sizes of objects for which reasonably secure estimates of the populations are available. We use the fractional population estimates of the sub-classes to combine together the P I grids for the individual components and determine a global map of P I for the Kuiper belt region (Figure 4).

COLLISION RATES

The "final" collision probability (yr -1 TNO -1 ) between two objects of radii r and R is obtained by multiplying Eq. 1 with the cross-sectional area of the two objects e.g. P F = P I × A = P I (r + R) 2 . We remind the reader that the factor of "π", in the cross-sectional area, is incorporated in P I .

To obtain the number of collisions N c , per unit time, for an ensemble of objects in a given size range, P F needs to be integrated over a suitable size distributions, representing the number of objects relevant for these populations. For example, the number of collisions N c , per unit time, that a single target of size R = 50 km would experience, with a given KBO population in a given size range r 1 < r < r 2 , scales as the intrinsic collision probability between the target and objects from the given orbit class as well as the size and number of impactors per unit size in that class. In a given impactor size range, the number of collisions onto a single 50 km target will be:

N c = P I r2 r1 (50 + r) 2 dN dr dr (2)
where dN/dr is the number of impactors with size between r and r + dr for the given KBO orbit class. Where as if one is interested in the total number of collisions between all members of two orbit classes then one must integrate over both size distributions:

N c = P I r R (r + R) 2 dN dr dN dR drdR (3)
where dN/dr amd dN/dR are the impactor and target size distributions. If one is interested in the number of collisions over a given time interval, equation 3 is then multiplied by that time interval.

The SFD for Kuiper belt populations can be approximated using exponential distributions of the form dN/dr ∝ r -q . The size distribution does not follow a single power law from the smallest to the largest bodies, but can be described using a small number of different logarithmic slopes, each valid for a given size range (e.g. Bernstein et al. 2004 et al. 2008, 2011;Fraser et al. 2014). For a given impactor KBO population, the number of collisions N ct , a target of R = 50 km, would experiences over a given time T will be the sum over the different logarithmic slopes appropriate to those impactors multiplied by the P I appropriate to the target and impactor orbital classes being considered. That is: where r 1i and r 2i is the range over which the logarithmic slope q i is valid and r oi is the normalization. The total number of collisions can be obtained by adding the contribution from each KBO population:

; Petit (a) (b) (c) (d) (e) (f) (g) (h) (i)
N ct = T P I i r2 i r1 i (50 + r) 2 r r oi -qi dr (4) (a) (b) (c) (d) (e) (f) (g) (h) (i)
N t = j N cj = T j P I j i r2 ij r1 ij (50 + r) 2 r r oij -qij dr (5)
where N cj denotes the contribution of a given KBO population to the number of collisions. Similarly, P I j is the intrinsic collision probability between a given target and KBO population. To provide a sense of the frequency of collisions in the Kuiper belt and to compare our calculated collision probabilities to previous works (e.g Greenstreet et al. 2015Greenstreet et al. , 2016) ) we present an estimate of the total number of impacts on to a R = 50 km target. Precise estimates of the total number of collisions a particular target of R = 50 km would experience depends, naturally, on the orbital elements of the target and the impactor population it could collide with as well as the gravitational focusing and acceleration due to non-negligible mass of the target body. By considering a small target (R = 50 km) we are able to neglect gravitational focusing. For our orbit we use mean collision probabilities for two target populations. We compute the total number of impacts that a mean "inner resonant object" (IRO, our combined 4:3 and 3:2 resonant populations) and a mean "Cold Classical Kuiper Belt Object" (CCKBO) would experience in the current Kuiper belt. The mean impact rates determined here should be roughly similar when compared to more direct computations for specific bodies.

In addition, considering a specific size range of impactors allows for approximate comparison with crater counts on Pluto-Charon system (Singer et al. 2019) and Arrokoth (Spencer et al. 2020). For Charon impactor size scales from crater size as d = 0.07D 1.151 where d and D are the impactor and crater sizes in km (Singer et al. 2019, S2). Using this relation we determine that the 1.4 km < D < 10 km craters measured on the smooth area of Charon's Vulcan Planum (Singer et al. 2019, Table S2) correspond, roughly, to impactors of 100 m d 1 km. For Arrokoth the New Horizons team reported craters in the size range of 200 m < D < 1 km, which translate to impactor sizes 2 m < r < 10 m. This choice of size range maps, roughly, to impact crater sizes (1.4 km < D < 10 km) (Singer et al. 2019) observed on Charon and to (0.2 < D < 1 km) (Spencer et al. 2020) on Arrokoth, using impact cratering scaling laws from (Zahnle et al. 2003), a relationship also adopted by (Greenstreet et al. 2015(Greenstreet et al. , 2016)). To determine the total number of impactors in these size ranges requires knowing the SFD of objects in each impactor group. [START_REF] Petit | [END_REF] provides population estimates for H g < 8, which we convert to H r using the color index (g -r) = 1.0 for the cold and (g -r) = 0.8 for the hot TNO components (e.g. Schwamb et al. 2019), respectively. The impactor sizes considered here are within absolute magnitude range 17 < H r < 32, in the "r"-band with the exact conversion dependent on albedo choice, discussed below. The number of objects between 8.66 < H r < 17 for both, high and low inclination TNOs was then estimated using a single logarithmic slopes of α = 0.4 or (q = 5α + 1 = 3), extrapolating smoothly from the N (H r < 8.66) in Table 2, i.e. we did not consider a divoted size-frequency distribution (Shankman et al. 2016) as was explored in Greenstreet et al. (2016Greenstreet et al. ( , 2019)).

For the range 17 < H r < 32 we use an average slope of α = 0.2(q = 2) based on the slope in the crater size-frequency distribution seen on Charon (q = 1.8 ± 0.2) (Singer et al. 2019) and Arrokoth (q = 2.3 ± 0.6) (Spencer et al. 2020). To convert absolute magnitude to impactor and target size, albedos of 0.14 and 0.08 were assumed for the low and high inclination populations.

Using populations sizes and collisions probabilities for our sub-population groups we then compute the total impact rate onto our nominal IRO and CCKBO objects. Our model prediction of the number of collisions, distributed by velocity and heliocentric distance, that a single target of size R =50 km, belonging to the inner resonant and CCKBO populations would experience over the age of the Solar System are presented in Fig. 5a and Fig 5b. These impact distributions differ, somewhat, from each other and provide an indication of where within the Kuiper belt each subpopulation is experiencing collisions and thus where those populations are most likely to be the source of local dust. In particular we note that for the IRO orbit a significant number of impacts occur interior to 40 au and those impacts inthe 40 au to 45 au zone have speeds of between 0.7 km/s and 3 km/s while for the CCKBO almost no impacts occur interior to 40 au and the bulk of impacts are in the 40-50 au zone with impact speeds well below 1 km/s.

For the case of a IRO, a 50 km spherical body (surface area ∼ 31,000 km 2 ), our model predicts N tot ∼12 impacts over the age of the Solar System, with impactors in the range 0.01< r <1 km (Fig. 5a). The New Horizons science team reported 70 craters with size 1.4 < D < 10 km on Charon's Vulcan Planitia smooth terrain (see supplementary material in (Singer et al. 2019)). Vulcan Planitia roughly subtends an area of A ∼ 35, 000 km 2 (Singer et al. 2019) or ∼ 10% larger than the surface area of 50 km target body. We find that the impact rate onto a typical IRO is dominated (∼ 35%) by cold classical TNO population, followed by the hot (∼ 21%) and detached (∼ 17%) sub-groups. Similar trend was observed by the analysis of (Greenstreet et al. 2015(Greenstreet et al. , 2016)). We find that the inner resonance population accounts for only ∼ 14% of the total impact rate on our nominal IRO. The offset between our predicted number of impacts and observed number of craters on Charon, aside from population estimate and SFD uncertainties is further affected by ignoring the following points:

1. We note that our model is designed to calculate the collision probability within and between different groups of TNO populations, as opposed to using the exact orbit of Charon. 2. In this work we ignore the full resonant and "Kozai-Lidov" dynamics of Charon's orbit, although we do add a correction factor of 1.5 (Greenstreet et al. 2015(Greenstreet et al. , 2016) ) to the collision rates of plutinos onto Charon. Both of these effects will add to non-uniform precession of the orbital elements, as well limit the range the orbital elements can span over a precession cycle, effectively resulting in an increased rate of encounters with other TNOs.

3. Failing to precisely model the relationship between impactor flux and number of observed craters. Impact scalinglaws suggest that, aside from the size of the colliding objects, there is a strong correlation between the impact speed, gravitational acceleration of the target and bulk densities of the target and impactor (e.g. [START_REF] Holsapple | Asteroid Impacts: Laboratory Experiments and Scaling Laws[END_REF]Housen & Holsapple 2003;Zahnle et al. 2003;Housen & Holsapple 2011;Housen et al. 2018). We convert observed impact crater diameters on Charon to impactor size, using an average collision speed of 2 km.s -1 for all populations. However, fig. 6 shows that the distribution of collision speeds, between some TNO populations, exhibit long tails with V I in excess of 3 km.s -1 . The effect of using an average impact speed, as opposed to a detailed distribution, in the crater to impactor size conversion will underestimate the size range of impactors considered in this work, which will translate to reduced total impactor flux.

For a detailed analysis of the impact rates onto a single object, the above shortcomings must be accounted for, as in (Greenstreet et al. 2015(Greenstreet et al. , 2016)). Recognizing the uncertainties, associated with TNO population estimates and in the SFD of objects smaller than 10 km, as well as neglecting the above mentioned points, our results cast confidence in the obtained collision probability calculations between different TNO classes. In a subsequent work, we will explore the actual collision rates in the trans-Neptunian region which urges for consideration of the above shortcomings, associated with our sample collision rate calculations.

For the case of the CCKBO, a target of radius R ∼ 50 km (surface area ∼ 31, 000 km 2 ) is expected to be impacted ∼ 1100 times, with impactors in the size range 2 m < r < 10 m, over the age of the Solar System. New Horizons team has counted 43 craters, over a surface area of 700 km 2 , on the encounter hemisphere of Arrokoth, with crater size range (0.2 < D < 1) km (Spencer et al. 2020). But the surface area, used for crater count on Arrokoth (700 km 2 ), is only 2.3% of that of a spherical 50 km target body, so our modeling predicts N tot = 1100 × 0.023 ≈ 25 craters on Arrokoth's crater count surface. Unlike the IRO case, we find a better agreement with observations. The reason is that the low inclination component of the main classical Kuiper belt is more tightly confined near the ecliptic, effectively resulting in higher collision probability within the cold population and with other TNO groups. Furthermore, Arrokoth is expected to be primarily impacted by the cold classical TNOs, which display a tighter impact speed probability distribution, where the assumed average impact speed V i = 0.3 km s -1 is good assumption for the overall spread of impact speeds. This is further supported by our finding that 56% of the impacts onto our nominal CCKBO are attributed to the cold classicals. We also find that the hot population is the next dominant impactor flux on our typical CCKBO, contributing ∼17% of the total impact rate, followed by the detached ∼10% and inner-resonant ∼8% populations. The population contributing the least to the total impact flux is the scattering population, accounting for ∼2%.

Finally, based on our modelling and analysis we find that most TNO collisions occur in the main classical belt, with the intrinsic collision probability for this sub-component accounting for ∼ 30% (see Tab. 1) of the average total collision probability P I ≈ 4.5 × 10 -21 km -2 yr -1 TNO -1 within the TNO population. Within the main classical belt, collisions are dominated by "Cold" on "Cold" impacts (∼17 %) with median impact speeds V I ∼ 0.3 km s -1 . Next most important source of impacts are the resonant TNOs (taking all resonances together) which account for ∼26 % of the total collision probability. Resonant TNOs experience collisions at speeds in excess of 1 km s -1 and may be more significant sources of dust than the impacts within the "Cold" population. In contrast, we find that the collision probabilities, within the scattering population, are almost four orders of magnitude lower, compared to that in the main classical belt (Tab. 1. See also Fig 1a and Fig. 3i). These results provide a strong proxy as to where in the Kuiper belt the dust production is expected to dominate. Models of dust production and general cratering among TNOs will need to take into account the differences in collisional probabilities between the two major sources of impacts.

DISCUSSION AND CONCLUSIONS

We have calculated the intrinsic collision probabilities and the corresponding impact speeds between different TNO populations, based on a revised dynamical model of the Kuiper belt, derived from the OSSOS ensemble sample, finding that low-velocity impacts at distances 30 < R h < 50 au are most likely.

Our calculations indicate that collisions are concentrated in the main Kuiper belt region (30-55 AU) but impacts over a wide range of distances (30-150 AU), are not excluded. Although our simulations indicate collision speeds in excess of 8 km.s -1 , most of the collisions are likely to happen with speeds well below V I < 5 km/s, over the age of the Solar System.

The highest intrinsic collision probability, P I ∼ 7.5 × 10 -22 km -2 yr -1 TNO -1 , is observed for impacts within the cold main classical population, with these collisions occurring at the lowest median speed V I ∼ 0.3 km.s -1 , and accounting for ∼ 17% of the average total collision probability P I ≈ 4.5 × 10 -21 km -2 yr -1 TNO -1 of all TNO populations. The dominance of impacts in the cold classical belt is not particularly surprising as these objects spend more time physically proximate to each other and close to the ecliptic, compared to the hot classical, scattering and resonant objects whose orbits carry them over large ranges of heliocentric distance, resulting in lower overall space density within these orbital classes. The low impact velocities among objects with the highest impact probabilities results from the narrower semi-major axis, eccentricity and inclination distribution of main classical KBOs, compared to the detached and scattering populations. Moreover, the highest impact speeds between the cold objects do not exceed 1 km.s -1 , raising a fair question as to whether these collisions are catastrophic or even crater forming. The outcome of a collision is a strong proxy as to the amount of material (dust) liberated in these impact events. Perhaps, due to the low impact speeds, collisions within the cold classical population are most likely to lead to crater formation rather than to a complete disruption of colliding bodies. Furthermore, the low impact speeds may also favor slow merging of cold KBOs into a single binary object.

The next highest collision probability (P I = 3.2 × 10 -22 km -2 yr -1 TNO -1 ) is observed between the hot and cold populations, with contribution of ∼ 7% of the total average intrinsic collision probability, though with noticeably higher impact speeds. In a subsequent work, we will investigate if these higher impact speeds are destructive, accounting for specific impact energy involved. Although somewhat lower in probability, collisions between hot and cold classical KBOs may very well dominate the dust production in the Kuiper belt, given their impact speeds.

We see a bi-modality in the velocity distribution for collisions between the scattering and cold populations, which could be a sampling "artifact" caused by the low-fidelity of our input scattering model. A few outliers with high inclination and eccentricity of the scattering model distribution over-contributing to the probability distribution. The gap in probabilities is not evident for all orbital populations and tends to be more prominent for impactors with high and semi-major axis dispersion, which the relative velocity is a strong function of. A few such orbits whose semi-major axis overlap due to numerical resolution can lead to granularity in the computed impact velocity distribution.

The lowest collision probability appears to happen within the scattering and detached populations with median collision speeds in excess of 2 km/s. Collision between these objects could result in complete destruction, depending on the size of colliding bodies. Given the low-intrinsic probability and total number of objects on these orbits, the likelihood of collisions, even after 4 Gyr of evolution in the TNO population, are now very low. Collisions among these populations may have produced significant levels of dust, early in the Solar Systems history when the number objects on such orbits would have been 10s to 100s of times higher.

Similar to Dell'Oro et al. (2013), we find that P I is indeed lower than previously reported values e.g., (Davis & Farinella 1997), though our computed collision speeds slightly differ from Dell 'Oro et al. (2013). We find that the cold main classicals collide on average with speeds V I ∼ 0.3 km/s, about 40% lower than the ones obtained by Dell'Oro et al. (2013). Similarly, our impact speed between the hot component of the main belt from 0.6 -4.4 km/s, with a median value of V I ∼ 2.0 km/s, whereas the value reported by Dell'Oro et al. (2013) is ∼ 20% higher. The most likely reason for that is that we report the median collision speed, based on the Öpik and Wetherill's formalism of collisions between asteroids, whereas the aforementioned authors provide results from best fit to Monte-Carlo simulations.

We also provide an approximate expected number of impacts onto a 50 km sample plutino and CCKBO and compared them to the observed crater densities on Charon and Arrokoth, by the New Horizons spacecraft. We find that, our results display a good agreement with crater densities reported by New Horizons, and also in agreement with the results first predicted by Greenstreet et al. (2015Greenstreet et al. ( , 2016Greenstreet et al. ( , 2019)).

Figure 1 .

 1 Figure 1. Intrinsic collision probabilities P I (colorbar) with a given collision speed V I at a given heliocentric distance R h between different dynamical KBO classes. The intrinsic collision probability in each cell is color coded and has units km -2 yr -1 TNO -1 . The values of the color bars only include values P I > 10 -28 km -2 yr -1 , as smaller values are unlikely to result in collision over the age of the Solar System (see text for details).

Figure 2 .

 2 Figure 2. Intrinsic collision probabilities P I (colorbar) with a given collision speed V I at a given heliocentric distance R h between different dynamical KBO classes. The intrinsic collision probability in each cell is color coded and has units km -2 yr -1 TNO -1 . The values of the color bars only include values P I 10 -28 km -2 yr -1 , as smaller values are unlikely to result in collision over the age of the Solar System (see text for details).

Figure 3 .

 3 Figure 3. Intrinsic collision probabilities P I (color bar) with a given collision speed V I at a given heliocentric distance R h between different dynamical KBO classes. The intrinsic collision probability in each cell is color coded and has units km -2 yr -1 TNO -1 . The values of the color bars only include values P I > 10 -28 km -2 yr -1 , as smaller values are unlikely to result in collision over the age of the Solar System (see text for details).

Figure 4 .

 4 Figure 4. Combined collision probability (color bar) as a function of the heliocentric distance R h and impact speeds VI , between all TNO populations (a) and the main Kuiper belt only (b). Subplot (a) is essentially the sum of all subplots in figures 1, 2 and 3, while subplot (b) presents the intrinsic collision probabilities between "cold-cold", "hot-hot" and "hot-cold" subclasses. Each 2-D pixel, in each subplot, has dimensions of (100 m/s × 1 AU).

Figure 5 .

 5 Figure 5. Number of collisions (color bar) onto an average "plutino" (a) and onto a typical "Cold Classical Kuiper Belt Object"(b). In both cases the target is assumed to have radius R = 50 km. The total number Ntot of impacts are estimated over 4 Gyr, with impactors in the range 100 m < r < 1 km (a) and 2 m < r < 10 m (b). Each 2-D pixel, in each subplot, has dimensions of (100 m/s × 1 AU). The discontinuity and small contribution to high-speed impacts (VI > 8 km.s -1 in both plots are attributed to the "granularity" in our scattering population model.)

Table 1 .

 1 Intrinsic collision probability PI and speeds VI by sub-component.

	Population	PI	Relative PI a	VI (km s -1 )
	Target -Impactor	km -2 yr -1	%	0.05 0.50 0.95
	Cold Classical
	Cold -Cold	7.5e-22	16.6	0.1	0.3	0.7
	Cold -Hot	3.2e-22	7.1	0.6	1.4	3.0
	Cold -Inner res	2.6e-22	5.8	0.5	1.3	2.7
	Cold -Main res	3.0e-22	6.7	0.4	1.0	2.1
	Cold -Outer res	1.6e-22	3.6	0.5	1.2	2.5
	Cold -Detached	7.7e-23	1.7	1.3	2.1	3.4
	Cold -Scattering 4.7e-23	1.0	1.8	2.6	3.5
	Hot Classical
	Hot -Hot	2.9e-22	6.4	0.6	2.0	4.4
	Hot -Inner res	2.3e-22	5.1	0.5	2.5	3.7
	Hot -Main res	2.5e-22	5.5	0.6	1.8	3.8
	Hot -Outer res	1.1e-22	2.4	0.7	1.8	3.8
	Hot -Detached	8.3e-23	1.8	1.1	2.4	4.7
	Hot -Scattering 4.4e-23	1.0	1.8	2.7	4.5
	Resonant
	Inner res -Inner Res. 3.2e-22	7.1	0.8	2.0 4 .3
	Inner res -Main Res. 2.6e-22	5.8	0.7	1.7	3.7
	Inner res -Outer Res. 1.0e-22	2.2	0.9	1.9	3.9
	Main res -Main Res. 2.6e-22	5.8	0.6	1.5	3.1
	Main res -Outer Res. 1.2e-22	2.7	0.8	1.7	3.2
	Outer res -Outer Res. 9.2e-23	2.0	0.7	1.6	3.1
	Scattering
	Scattering -Inner Res. 6.1e-23	1.4	1.6	2.8	4.7
	Scattering -Main Res. 5.5e-23	1.2	1.6	2.7	4.1
	Scattering -Outer Res. 3.3e-23	0.7	1.7	2.8	4.1
	Scattering -Detached	2.3e-23	0.5	1.5	2.6	4.0
	Scattering -Scattering 1.6e-23	0.4	1.7	2.8	4.5
	Detached/Outer
	Detached -Inner Res. 7.8e-23	1.7	1.2	2.2	3.9
	Detached -Main Res. 8.2e-23	1.8	1.2	2.2	4.0
	Detached -Outer Res. 5.2e-23	1.2	1.1	2.2	4.0
	Detached -Detached	3.8e-23	0.8	1.1	2.2	4.4
	a Fraction of total intrinsic probability, not including relative population sizes.

Table 2 .

 2 TNO sub-component population estimates.

	Orbit Class	Fraction	Population
		% of total	N (Hr < 8.66)
	Cold Main Belt	20	61000 ±6000
	Hot Main Belt	12	38000 ±9000
	Detached/Outer	26	80000 ±40000
	Scattering	29	90000 ±20000
	Inner Resonances
	4:3	0.1	400 +800 -120
	3:2	4.2	13000 +6000 -5000
	Main Resonances
	5:3	1.6	5000 +5200 -3000
	7:4	0.9	3000 +4000 -2000
	Outer resonances
	2:1	1.2	3700 +4400 -2400
	7:3	1.3	4000 +8000 -3000
	5:2	3.8	12000 +15000 -8000
	References-Petit et al. (2011),Lawler et al.
	(2018b),Gladman et al. (2012)

We do not consider the effects of resonant dynamics of the colliding objects as these do not strongly alter inter-class collisional probabilities.