Pluto and Charon Impact Crater Populations: Reconciling Different Results

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.

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.

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. 2018a, 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.

Robbins mapped craters on the PDS-released version 1 300 m pixel⁻¹ 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).

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⁻¹ 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).

Zoom In Zoom Out Reset image size

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⁻¹. 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_LORRI_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.

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⁻³ 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⁻³ 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 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 highest-confidence 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).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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⁻¹ 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.

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₁₀–log₁₀ 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⁻³ slope would be a flat line on an RSFD, a more negative slope (such as D⁻⁴) is referred to as "steeper," and a more positive slope (such as D⁻²) 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 least-squares 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. (2019, 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 (Singer et al. 2019 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).