Pluto and Charon Impact Crater Populations: Reconciling Different Results

The size–frequency distribution (SFD) of objects across the solar system is influenced by their formation and dynamic evolution. Bodies in different areas of the solar system—especially the inner versus outer—have likely experienced different histories, so measuring that SFD across orders of magnitude of solar distance can elucidate those different pasts. Accurate measurement of the SFD is greatly assisted by impact craters, formed from the smaller bodies that strike much larger ones, and it is often easier to measure impact craters’ SFDs than the smaller impactors themselves. One of the main results from New Horizons’ Pluto–Charon encounter was the crater SFD measurement, which illustrated an unexpected SFD with a distinct change in slope near ≈10–15 km diameter features. However, different methods of reporting these impact craters have resulted in some questions in the community about what the SFD actually is of Pluto's and Charon's crater population and therefore what the SFD is of the impactors that formed them. We have performed new crater population studies on both bodies, compared them with previously published work, and demonstrated that there is no ambiguity with respect to whether there is a transition to a shallower SFD population near ≈10–15 km. We find that this shallower slope is approximately −1.7 ± 0.2, while the steeper slope for larger craters is approximately −3.8 ± 0.6, though both of these uncertainties must be considered absolute minima based on fit sensitivity to the exact parameters and craters used. These values are consistent with previously reported results.


Introduction
NASAʼs New Horizons mission to the Pluto-Charon system and the Kuiper Belt provided the first well-resolved images of objects that spend most or all of their orbits within the Kuiper Belt. These images provided the first observations of the crater populations on these worlds without the complication of planetocentric impactors driven into moons around the gas and ice giant planets of the outer solar system. Robbins et al. (2017) and Singer et al. (2019) found that Plutoʼs and Charonʼs crater populations demonstrate a distinctly different power-law slope for craters smaller than ≈10-15 km in diameter (D), where the transition point varies somewhat on different terrains of Pluto and Charon. Such an observation was not generally predicted, and it has potentially wide-reaching implications for planetesimal formation and collisional evolution of the Kuiper Belt (e.g., Greenstreet et al. 2015Greenstreet et al. , 2016Greenstreet et al. , 2019Nesvorný et al. 2018Nesvorný et al. , 2019. Robbins et al. (2017) reported a break in slope but to a power-law exponent that was not nearly as shallow as that of Singer et al. (2019). Unfortunately, this led some researchers to think that there was a disagreement and state, "The reason for this discrepancy within the New Horizons team is not clear" (Morbidelli & Nesvorý 2020).
Given that this change in slope is incredibly important for understanding the dynamical evolution of the Kuiper Belt, another look at the crater data is warranted. In this work, we discuss a reassessment of the crater data, and we review how the Robbins et al. (2017) catalog was generated and how the slopes were reported therein. This reconciliation is supported by a new crater mapping effort by Robbins, reported herein, and we found that there is a distinct change in the crater sizefrequency slope at D ≈ 10-15 km.

Methods
A substantial reason for the apparent mismatch between the 2017 and 2019 results is due to how the data were gathered and reported, and the new work described herein was generated differently. Here we detail the methods used in each effort. Robbins et al. (2017) Despite the publication date of 2017, fully two years after New Horizons' encounter with the Pluto-Charon system, work on a crater database began immediately following resolved images being downlinked from the spacecraft. The goal of this effort was a global crater database, or at least a crater database covering all visible, resolved terrain. Crater rims were traced on raw, minimally calibrated images from the spacecraft (e.g., radiometric corrections were applied), and later those rim vertices, in units of pixels, were converted to units of decimal degrees on each body through the CAMPT program in the United States Geologic Surveyʼs Integrated Software for Imagers and Spectrometers (ISIS; Anderson & Sides 2004). Doing this conversion required accurate pointing information, which evolved over the time from the initial downlinks to when the data were published, and it continued to change after. CAMPT also required a proper camera model, and it was discovered roughly a year after the encounter that the framing camera model distortion was inverted in ISIS, meaning further changes had to be made. It was thought that tracing rims on raw images would obviate the need to continuously modify the rims with each new map projection and coordinate change. In addition to those traces, other researchers mapped craters on individual-but map-projected-images and mosaics that were produced by the New Horizons team in later months as derived map products were made.

Methods from
All of those crater identifications were merged together through a density-based spatial clustering algorithm (specifically, the density-based spatial clustering of applications with noise, or DBSCAN, from Ester et al. 1996) adapted for craters (Robbins et al. 2014). Any cluster analysis has free variables that must be tuned to the problem at hand, and there will almost always be problems (i.e., features that are merged but should not be, and features that are not merged that should be). The clustering was necessary because multiple images covered the same terrain, so craters were duplicated, and multiple researchers identified the same features.
After the clustering algorithm was applied, all clusters were then manually examined on the global mosaics-which had been assembled by that time-and manually edited if needed based on whether there appeared to be errors with features clustered that should not have been or unclustered features that should have been. From the final clustering, average locations and diameters were output, along with a subjective confidence of whether the feature was a crater or not, which was also an average of the individual identifications' confidences that went into each cluster. All slope values reported in Robbins et al. (2017) were based on craters with a confidence of 3, 4, or 5, which goes from a "50/50 chance" up to "certain." In contrast, the comparisons presented later in Section 4 all use the highest confidence value (5) unless otherwise stated.

Methods from Singer et al. (2019)
Singer et al. (2019) mapped craters on specific image sequences (e.g., a mosaic of images taken at one point in the encounter) or image scans as map-projected products produced by the New Horizons team. This produced one crater data set per image sequence, or scan, that was at an approximately consistent image resolution (because the spacecraft moved during scans and sequences, slight image scale differences are present). Then, each data set was further subdivided by geomorphic province. Poorly imaged areas were excluded. When the version 1 base maps were released Planetary Data System (PDS; Schenk et al. 2018aSchenk et al. , 2018b, all craters in each data set were checked and shifted as necessary so that they would align with the published base map in latitude and longitude.

Methods Used in This Work
Robbins mapped craters on the PDS-released version 1 300 m pixel −1 base maps and higher-resolution internal products with the goal of mapping the well-imaged encounter region of Pluto and the well-defined equatorial crater population of Charon. Rims were traced on these mosaics in units of decimal degrees with circles and ellipses fit following the extensive discussion in Robbins (2019). Robbins et al. (2017) Often, especially when observing something new for the first time, there is a strong desire to extract the greatest amount of information as soon as possible. Unfortunately, this desire must be tempered with caution to ensure that errors do not build up from a potentially rapidly changing source. Despite best attempts, issues crept into the 2017 work.

Sources of Error from
First, perhaps the most obvious potential issue is that the clustering did not work as well as it has in the past, for it is sensitive to the clustering criteria (the criteria being the distance between two craters scaled by the mean diameter, and the size difference between those two craters relative to their mean diameter). Because the three primary crater analysts used very different products, craters could project to very different locations, several crater diameters apart from each other. The philosophy with the final database was to trust the initial identifications, meaning that two high-confidence features ∼5D away from each other would remain as separate craters in the database rather than be merged together, even if there was not an obvious corresponding feature for one of them. In hindsight, this did not work well given the large disparities in positions as the map products and control networks progressed.
Second, the changing coordinate system during the process meant that craters might not align with the current datum. The differences were up to tens of kilometers. That fed into some of the reasoning behind the first issue. Indeed, even when Singer et al. (2019) were constructing their database, the coordinate system changed yet again, and craters were manually shifted to the latest base map near the time of publication.
That leads to a third issue where, because Robbins' rim traces were fit to circles after converting from pixels to decimal degrees, if the projection was off by a few degrees, then the fit could return a crater of a very different size. For example, if a crater 10 pixels wide at 30°latitude is instead projected to 50°l atitude (some initial shifts were, indeed, that significant), it will be only 75% as large, simply from the projection change. Even if averaging with other crater annotations from the clustering, it will still offset the crater from its true diameter.
A fourth problem is the use of confidence values in the database, how the values were reported, and how they have been subsequently used (or not used) by other researchers. If the confidence thresholds are different, such as raised to only accept values of 5 when performing slope analyses, the slopes from Robbins et al. (2017) will be different: they will be closer to those from Singer et al. (2019). This is especially important on Pluto, where geologic context is critical to interpreting what features may or may not be impact craters. Singer et al. (2019) had the benefit of at least two additional years of the communityʼs interpretation of these terrains.
Fifth, Robbins et al. (2017) reported on slopes from terrains and image scales that varied within single units. While care was taken to try to identify only single geologic units, the benefit of several more years and dedicated mapping efforts White et al. 2020) allowed for more specific terrains to be used in Singer et al.ʼs work. Merging of multiple terrains that might have experienced different erosion histories will confuse the crater population and lead to erroneous results. Similarly, combining craters from multiple image scales and not accounting for a single conservative crater completeness level will act to shallow small craters, which will affect the overall slopes that are calculated.
All of these effects combined are not intended to fully negate the analysis of Robbins et al. (2017) or indicate a faulty crater database, but they do mean that much more caution is needed when using that work than one might otherwise think. Additionally, these issues should be considered "lessons learned" for other crater analysts and mission teams in the future.

Regions for Comparison and Study
To compare the data between the three works, we must compare the same regions. Singer et al. (2019) was the more specific of the two published works, using distinct images or image sequences from the New Horizons Pluto-Charon encounter (which eliminates the first issue with Robbins et al. 2017, discussed above) and more specific mapping (which eliminates the fifth issue discussed above), while Robbins et al. (2017) used more vague shapes based on very early mapping and nomenclature boundaries. Figure 1 shows monochromatic base maps for Pluto and Charon with regions of interest outlined for this work. For Charon, the geologic map from Robbins et al. (2019) was used to identify the Vulcan Planitia boundary, which broadly corresponds with that unitʼs boundary used in Singer et al. (2019), though it eliminates the eastern area mapped as "Rough Terrain" in the 2019 geomorphologic map. The higher-resolution framing camera strip at ≈160 m pixel −1 through the middle of that geographic province was used as the second region, and we also removed a heavily pitted region in the southern half of the high-resolution strip to create a third region for direct comparison of values from Singer et al. (2019).
The larger Cthulhu region on Pluto crosses several different image sequences, but it is primarily made from the line scan rather than framing camera due to its better signal-to-noise ratio properties in practice. It cuts off in the south due to the solar incidence angle, east because of switching to the LORRI framing camera images and complications from the vast glacier covering the Sputnik basin, north due to a geologic unit change, and west due to unfavorable emission angles; the boundary herein is similar to but contracted from the boundaries used in both previous works. Terrain near the eastern boundary is particularly fraught with uncertainty in crater identification due to both the chaotic nature of the terrain and the extreme brightness differences across the terrain making crater identification difficult, similar to the crater detection issues on Iapetus's bright/dark transitional terrain. The three other regions are better-scale LORRI strips, taken at ≈75-235 m pixel −1 . Each of these three strips, however, was further cropped along approximate northeast-southwest trending diagonals to bracket the regions of solar incidence (to the northwest) or transition to Sputnik (to the southeast). For the Pluto comparisons discussed in the next section, the crater data from Singer et al.ʼs data set labeled MPAN that intersected each region in Figure 1 were used. The "MPAN" data correspond to the PDS "Request ID" field, PEMV_P_MPAN1, and they are also further labeled by geographic subsets (see Singer et al. 2019, for details). For the Cthulhu region, we used data from the subset "MPAN_Cthulhu" only; for the LEISA ride-along, we used data from the "MPAN_Cthulhu," "MPAN_Fretted," and "MPAN_Mid-lat" subsets (resulting in a slight lack of overlap; see Figure 1 and note 4 in Table 1); for the mid-latitude LORRI strip, data from the "MPAN_Burney," "MPAN_Fretted," and "MPAN_Mid-lat" subsets were used; and for the northernmost LORRI strip, data from "MPAN_Burney" and "MPAN_North" were included. Robbins mapped craters on the MVIC data, as well as high-resolution LORRI data (PELR_P_LORRI_MVIC_CA, PELR_P_MPAN_1, and PELR_P_LEISA_HIRES) for this work but was not as specific for cutoffs and data provenance.
For the Charon comparisons discussed in the next sections, data from the subset "LORRI_CA_VP" (Request ID "PELR_C_LOR-RI_MVIC_CA") were used for the two smaller-area comparisons, and data from "MVIC_CA_VP" (Request ID "PEMV_C_LORRI_ MVIC_CA") were used for the larger area that covers the majority of Vulcan Planitia.

Size-Frequency Distribution Comparison
One of the primary methods to compare crater populations is through the venerable size-frequency distribution (SFD), first codified in a 1978 report that was formalized in the Crater Analysis Techniques Working Group (1979). The graphing methods are variations on the theme of plotting crater diameter on the abscissa and crater number (in some form) on the ordinate axis. The two recommended versions are the cumulative (CSFD) and relative (RSFD), which plot the sum of craters from large to small as the ordinate value, or the differential number normalized to a D −3 slope as the ordinate value, respectively. The latter form has the benefit of seeing how small-crater populations compare without being masked by the larger craters, and it has the benefit of making better use of graph space by removing the common ∼D −3 slope. We use the updated format from Robbins et al. (2018), which represents each crater as a probability density function (PDF) based on the idea that diameters are not exact or perfectly repeatable, and it sums those PDFs to make an empirical distribution function (EDF). The abbreviations are thus CSFD EDF and RSFD EDF . Figure 2 shows CSFD EDF (top) and RSFD EDF (bottom) for the four regions on Pluto, while Figure 3 follows the same Note. See the text for exact images used. Only the power-law slope (exponent) is shown. Note (1)-To illustrate the effects of fitting to a truncated Pareto distribution, we report parenthetical slopes when using D max = 2× the largest crater in the region in that data set. Using D max = 1× the largest crater shallows all slopes by up to 1. Note (2) format for the two regions on Charon. Each shows a comparison between the highest-confidence data from Robbins et al. (2017), the data from Singer et al. (2019), and the highestconfidence data from this new work, all with 95% confidence intervals (choosing a tighter confidence threshold in this new work does not alter any conclusions). Importantly, all three SFDs for each region are extremely similar to within the confidence regions, with the Charon data more similar than the Pluto data (except for Cthulhu, possibly due to the improved signal-to-noise ratio of the MVIC scan). On Pluto, the largest differences between the three crater databases are in the LEISA ride-along plot and the mid-latitude LORRI strip plot, and those differences are only in the Robbins et al. (2017) data relative to the two newer data sets. For the LEISA ride-along, this can be interpreted as the earlier work going to smaller crater diameters than either of the more recent ones, for the SFD slope is fairly stable down to D ≈ 2 km. In the mid-latitude LORRI strip, the 2017 database shows a deficit of craters D ≈ 3-5 km compared with the two more recent sets. Looking at the images with the data overplotted does not reveal an obvious reason for this difference.
On Charon, the primary difference between the three data sets is on the larger Vulcan Planitia province (as opposed to the LORRI high-resolution strip), where the new database shows more craters in the ∼2-3 km range than the previous two. This work used 300 m pixel −1 base maps, limiting reasonable crater detection to ≈2.4-3 km (≈8-10 pixels; Robbins et al. 2014). Truncating at 3-4 km removes the deviation.

Slope Comparison
One of the most important results from the New Horizons mission was the unpredicted, distinct, unambiguous change in SFD slope near D ≈ 10-15 km. Since SFDs are plotted on log 10 -log 10 axes and crater populations tend to follow a power law, they appear as a straight line on most SFDs. The exponent in the power law appears as a linear slope on these graphs. By convention, the slope refers to the fit to a differential SFD. A D −3 slope would be a flat line on an RSFD, a more negative slope (such as D −4 ) is referred to as "steeper," and a more positive slope (such as D −2 ) is referred to as "shallower." The slopes can be fit using a maximum-likelihood estimate (MLE) fit to a Pareto distribution (for large craters) that is only bounded at the small end, and the data can be fit with an MLE fit to a truncated Pareto distribution (for small craters) that is bounded on both ends (Aban et al. 2006; a Pareto distribution is the same as a power law). All reported uncertainties in this work are 1σ. An important point here is that, in the limit that D max → ∞ on the truncated fit, the truncated and nontruncated fits will return the same results. However, one cannot have an infinitely large crater on a surface, so some reasonable limit must be placed. Unfortunately, using the largest crater as that limit precludes the fit from considering the absence of larger craters, which, if that absence is considered, steepens the fit (i.e., evidence of absence when those features could have been observed if they existed is still information that will, and should, affect the result; this is an important distinction from much more common leastsquares fitting methods, which can only use the existing data). While we found no case where the truncated versus nontruncated fits produced results different by >1σ, this should be kept in mind when interpreting the large-crater slopes.
The break to a shallower small-diameter slope was strongly suggested by Robbins et al. (2017), who identified smaller-crater slopes near −1.4 to −2.5 (depending on body and terrain) and large-crater slopes of −2.9 to −4.4 (or as steep as −8.4, which is likely spurious due to Sun angle limitations). Singer et al. (2019Singer et al. ( , 2021 focused more on determination of slopes per image data set and region, reporting a typical value near −3 for larger craters on both Pluto and Charon (some were steeper, some were shallower) and −1.7 for smaller craters on Vulcan Planitia when explicitly excluding a heavily pitted region near the terminator (see Table S2 of Singer et al. 2019 for all regions and slopes). Singer et al. (2019) also looked specifically at the data sets for Vulcan Planitia and determined that the break in diameter for this region was close to D ≈ 13 km.
Those results could be interpreted as different, despite reported 1σ uncertainties of ≈±0.2-0.8 in the former and ≈±0.1-1.0 in the latter. However, taking the steepest value from the 2017 work for smaller craters (−2.5 ± 0.2) versus the shallowest value from the 2017 work for larger craters (−2.9 ± 0.2) across all terrains on either body, one could interpret those as the same, thus negating the slope change and calling the Singer et al. (2019) results that identified that inflection point into question.
We recalculate the slopes using the Singer et al. (2019) data and new data in this work, giving margin to a nominal ∼15 km transition; the only variable is the minimum diameter for small craters due to differences in estimated minimum completeness. We provide two sets of small-crater fits, the first with a D max cutoff of 10 km and the second with a D max cutoff of 12 km, with the goal to remove effects from some of the ambiguity of where the transition lies; our D min for larger craters in this work is 18 km, thoroughly avoiding the transition diameter. This is an important difference from both previous works, where Robbins et al. (2017) "eyeballed" different diameter cutoffs for fitting, and Singer et al. (2019) mostly used a 10 km break diameter (in two cases using a 15 km break diameter), where all craters smaller than the break were fit and then all craters larger than the break were fit. The new fits are reported in Table 1 and illustrated in Figure 4. Table 1 and Figure 4 demonstrate a large range of fit values, and there is a suggestion that slopes from craters in this work are slightly steeper than those from the 2019 work (in 15 of 18 cases). However, there are two important conclusions from the table: (1) the slopes are consistently, statistically, significantly shallower for small-diameter craters than large-diameter craters; and (2) the slopes from both Singer et al.ʼs (2019) data and this new work are all within 1σ of each other when comparing the same diameter ranges for the same surfaces in the two crater databases, despite the preference for Robbins to measure craters that yield slightly steeper SFD slopes.
With that in mind, and points 1 and 2 above robustly determined, it should be noted that slightly modifying the diameter range of craters that are fit can alter the fit values. For example, despite smaller numbers, the LORRI high-resolution strip across Charon has the best-defined crater population in this work and the smallest diameter completeness of ∼2 km (Figure 2; note that Singer et al. 2019 used a minimum diameter of 1.4 km for this areaʼs nonpitted region). In Table 1, the fits reported for D = 2-10 km are −1.6 ± 0.3 and −1.9 ± 0.2  and this work, respectively). Changing the fit to D = 2-12 km yields slopes of −1.6 ± 0.2 and −1.8 ± 0.2, respectively, which are within each otherʼs 1σ range again. However, changing the fit to D = 4-12 km yields −1.4 ± 0.5 and −1.8 ± 0.5 for the two, which is a nontrivial (though still within 1σ) shallowing for the former dataʼs SFD. Using D = 4-10 km yields −1.5 ± 0.7 and −1.9 ± 0.7, respectively. These changes are due primarily to the relatively few craters fit (on the order of tens), the uncertainty of the completeness limit, and the break in slope near D ≈ 10-15 km. Also, while these are all within each others' 1σ confidence values, it is still an important point of the approximate nature of these slopes that their associated uncertainties must be considered when using them for any derivative work, such as trying to understand the impactor population.
For completeness and as direct a comparison as possible, we can similarly look at the fits in Table S2 of Singer et al. (2019) and the almost exact overlap in areas for the LORRI_CA_VP data (the Charon high-resolution data in our Table 1). They found 16 craters of D 10 km, as did Robbins in this new effort. They fit −2.8 ± 0.8 to the data using an MLE fit. To return that fit value, Singer et al. (2019) used a truncated Pareto fit, bounded as (10 km, D max measured ). If instead, we use an upper bound of 2D max measured , as in Table 1 of this work, the fitted slope is −3.4 ± 0.6, while an unbounded fit (D max = ∞) is −3.6 ± 0.6. Performing the same two bounded and one unbounded fits with Robbins' 16 craters in this work yields −2.8 ± 0.8, −3.5 ± 0.7, and −3.6 ± 0.6, respectively. To these two significant figures, the values are identical, except for the middle one.
As a second direct comparison, we can look at Singer et al.ʼs "smooth area" subset. The mapping is very slightly different in this work, as they found N = 70 craters, while our new subset from their work is 71 craters. They fit [1.4 km, 10 km] as the range and found −1.8 ± 0.2. With the slightly revised crater subset, fitting the same range, we find −1.8 ± 0.2 from their data and −1.9 ± 0.2 from the new Robbins data. Incidentally, when trying to match the mapping between the 2019 work and this, we initially had a slightly larger area to the south that included one extra, D ≈ 9.0 km crater (N = 72 from Singer et al.ʼs work). The inclusion of that one crater shifted the fit to −1.7 ± 0.2 (to three significant figures, it was −1.73 ± 0.22; without that crater, the fit is −1.78 ± 0.22). This is a further example that the exact craters that are fit can have the ability to appear to nontrivially alter the fit, even though these are again clearly within each others' 1σ and even when it might seem as though N is large (∼70 in this case).

Summary and Discussion
The population of impact craters on a planetary surface can tell us a lot about the processes that formed and modified them. Plutoʼs and Charonʼs surfaces offer the first relatively isolated look at the crater population in the outer solar system that is uncomplicated by planetocentric impactors around the gas and ice giants. Robbins et al. (2017) published the first global  databases for each body, while Singer et al. (2019) focused on the SFD slopes and implications for Kuiper Belt object formation and evolution. However, the slopes and graphs reported in each paper can appear to disagree, resulting in problems for researchers who want to use the results for their own work and potentially casting some doubt on whether the break in slope is real.
In this new work, Robbins undertook a new crater identification and measurement effort using a more stable image base and having the benefit of half a decade of work on understanding Plutoʼs and Charonʼs geologic context. As illustrated in Figures 2 and 3, this new work, compared with the highest-confidence craters from Robbins et al. (2017) and the Singer et al. (2019) data, shows remarkable overlap in the SFDs. Not only that, but the slopes are comparable when fitting over the same diameter ranges. This means that the results of a break in slope to a shallower population for D  10-15 km is certainly robust. As reported in Singer et al. (2021), we would need to undercount ∼3 km impacts by a factor of ∼10 in order to maintain the large-crater slope.
However, what exactly that slope value is is somewhat uncertain. With on the order of tens of craters, slight shifts in the diameters that are fit can change the results by what can appear to be nontrivial amounts. Singer et al. (2019) discussed complications from the geology of Pluto, and both prior works addressed issues with Sun angle, signal-to-noise ratio, and extreme reflectivity variations potentially masking impacts. Therefore, neither the 2019 work nor this one posit that the crater SFD slopes on Pluto necessarily represent an original, unmodified, and fully detected crater-and therefore impactor-population.
Therefore, Charonʼs is likely both a less modified and more easily measured population than Plutoʼs. The better number statistics across Vulcan Planitia, rather than the smaller-area subset in the high-resolution LORRI strip, are tempting and would indicate that the shallow smaller craters' slope is roughly −1.9 ± 0.2, and it indicates that the steep larger craters' slope is roughly −3.8 ± 0.6. (These values use a weighted mean of the fits carried out to more significant figures than shown in Table 1, then truncated to two significant figures; the uncertainty is an average uncertainty rather than added in quadrature, which would reduce the overall uncertainty.) However, Charonʼs high-resolution strip across Vulcan Planitia is likely to be more accurate for smaller craters, given that better-resolved features make crater identification more certain, and the better pixel scales allow reasonably confident measurements to smaller diameters. We therefore think that it is most representative of the small impactor flux of Kuiper Belt objects, as discussed in Singer et al. (2019). As discussed in Section 4.3, the slopes calculated are variable based on just changing the diameter cutoffs by ∼2 km at either end and inclusion or exclusion of even a single crater. But the slopes tend to center around a slope of ≈−1.7 ± 0.2 (using a weighted mean of D max = 10 or 12 km gives the same results to two significant figures). This −1.7 is the value preferred by Singer et al. (2019), and it is substantially shallower than any value stated in Robbins et al. (2017) for reasons described above.
Independently, we can look to Nixʼs craters to see if they provide a supporting constraint. Nix was the only small satellite of Plutoʼs four to be imaged with sufficient resolution to reliably detect numerous impact craters. Robbins et al. (2017), which included data from Weaver et al. (2016), reported craters on two images of Nix, and here we use those from the MVIC scan, which had more favorable Sun for crater detection. Using only the highest-confidence craters, fitting a power-law distribution with only a lower bound of 3 km (N = 9 craters) yields a −2.3 ± 0.5 slope, while a truncated fit with D max = 2× the largest identified crater (14.1 km) yields −1.9 ± 0.7. Assuming completeness to the next-smallest crater, D min = 2.4 km, yields slopes −2.3 ± 0.5 and −2.0 ± 0.6. Using confidence 4-5 craters (N = 16), the four slopes are −2.6 ± 0.5, −2.4 ± 0.6, −2.4 ± 0.4, and −2.2 ± 0.5. Using only data from Weaver et al. (2016), the fit truncated with D min = 2.4 km yields −2.0 ± 0.5 bounded and −2.3 ± 0.4 unbounded when using all-confidence data; when using just high-confidence data, the two are identical (their lower-confidence craters are all <2.4 km).
The purpose of these eight Nix fits is to again emphasize that there is variation with small N when choosing what to fit, and one can practically cherry-pick the value one desires and be able to reasonably justify it. Attempting to not do that, we simply use the Nix experiment to demonstrate that Nixʼs craters are generally consistent with a shallow slope for the small impact craters in the system.
In summary, the data from 2017 are generally consistent-when using high-confidence features-with the data from 2019 and the new data gathered here. Plutoʼs impact crater population is affected by various geologic processes and practical imaging effects that make it difficult to discern the small-crater population. However, Charonʼs craters in the Vulcan Planitia region are not affected by either of those problems, and small, kilometer-scale craters would need to be consistently undercounted by more than an order of magnitude by multiple researchers for there to not be a break in SFD slope for D ≈ 10-15 km. Therefore, we reaffirm this slope change. However, exactly what the differential SFD slope of craters larger than D ≈ 10-15 km is, and what the slope of craters smaller than D ≈ 10-15 km is, is somewhat uncertain. We can robustly say that they are somewhere around −3.8 ± 0.6 and −1.7 ± 0.2, respectively, but given the changes in these values based on exactly how the data are fit, these error bars must be understood to be conservative. These SFD results are important for understanding the population of the Kuiper Belt, how it formed, and how it has evolved (e.g., Gladman et al. 2001;Zahnle et al. 2003;Kenyon & Bromley 2004;Bierhaus & Dones 2015;Greenstreet et al. 2015Greenstreet et al. , 2016Greenstreet et al. , 2019Singer et al. 2019;Morbidelli et al. 2021;Parker 2021), and any modeling must take into account the uncertainty in the observational crater data.