The DECam Ecliptic Exploration Project (DEEP). V. The Absolute Magnitude Distribution of the Cold Classical Kuiper Belt

Beyond the orbits of the major planets, our solar system hosts a large population of minor bodies known as Kuiper Belt objects (KBOs). In the 30 years since the observational establishment of the Kuiper Belt (Jewitt & Luu 1993), several surveys (e.g., Millis et al. 2002; Bernstein et al. 2004; Petit et al. 2011; Bannister et al. 2018; Bernardinelli et al. 2022) have pushed the inventory of known objects to nearly 4000. These bodies, which are left over from the birth of our planetary system, provide constraints on its formation and dynamical evolution. When taken in aggregate, their dynamics, compositions, and sizes enable us to infer details about the dynamical evolution of the planets, the composition of our solar system's protoplanetary disk, and even the physical processes by which planetesimal formation occurred (see, e.g., Nesvorný 2018; Gladman & Volk 2021).

In particular, the size distribution of the so-called cold classicals (CCs)—which are thought to be relics of the birth of the solar system, relatively untouched and uncontaminated in the ∼4.5 Gyr since their formation—is a sensitive probe of the process of planetesimal formation. If the CCs are truly a quiescent population, a measurement of their size distribution can provide us with a unique opportunity to directly compare a primordial size distribution with predictions made by planetesimal formation models. Such a comparison will enable us to hone our formation models, and better understand the details of the physical processes at play in planetesimal formation, as well as the specific conditions of our own protoplanetary disk.

Recently, the streaming instability (SI; Abod et al. 2019) has begun to emerge as a leading theory of planetesimal formation. Numerical SI simulations predict an exponentially tapered power-law absolute magnitude (H) distribution, enabling a direct comparison between theory and observation. Kavelaars et al. (2021) found that the absolute magnitude distribution of the CCs detected by the Outer Solar System Origins Survey (OSSOS; Bannister et al. 2018) is consistent with an exponentially tapered power law. However, the CCs used by Kavelaars et al. (2021) only went as faint as H_r ∼ 8.3, leaving faint-end consistency with the SI as an open question. Existing literature on measurements of the faint end of the CC absolute magnitude distribution seems to be in weak tension with SI models of planetesimal formation. In particular, Fraser et al. (2014) find a marginally steeper faint-end size distribution than is predicted by SI simulations. However, a dearth of observed objects at the faint end of the distribution, along with the fact that Fraser et al. (2014) did not fit to an exponentially tapered power law, limits the usefulness of such comparisons. We require a larger, deeper set of CC detections by a survey with well-understood biases to thoroughly test any planetesimal formation theory.

The DECam Ecliptic Exploration Project (DEEP) is the first survey with sufficient depth and areal coverage to settle the tension in the shape of the faint end of the H distribution of the CCs. In this paper, we analyze data from 20 DECam fields (comprising an area of approximately 60 deg²), reaching a mean limiting magnitude of m_r ∼ 26.2. We use a subset of our data to reconstruct the luminosity function of the Kuiper Belt as a whole, as well as the luminosity function of the CC population. As the main scientific result of this paper, we use our results to reconstruct the underlying absolute magnitude distribution of the CCs, and find consistency with models of planetesimal formation via the SI (Abod et al. 2019; Kavelaars et al. 2021).

This paper is organized as follows. In Section 2 we outline our observational strategy and discuss the data used in the subsequent analysis. Next we describe our image preprocessing pipeline (Section 3) and the pipeline used to carry out the object search (Section 4). We present an overview of our detections in Section 5. In Section 6 we calculate our detection efficiency using implanted synthetic objects. We compute the luminosity function for our KBOs as a whole in Section 7. In Section 8 we isolate a sample of CCs from our detections. Our main scientific result, presented in Section 9, is a calculation of the absolute magnitude distribution for the CC population. In Section 10 we calculate an estimate of the total mass of the CCs. In Section 11 we test the consistency of our absolute magnitude distributions with the results from Bernstein et al. (2004). The paper concludes in Section 12 with a summary of our results and a discussion of their implications.

DEEP was carried out with the Dark Energy Camera (DECam) on the 4 m Blanco telescope located at Cerro Tololo Inter-American Observatory in Chile from 2019–2023, targeting four patches of sky along the invariable plane (see Trilling et al. 2024; Trujillo et al. 2024; hereafter Papers I and II, respectively, for more details). This paper focuses on the data taken in one of those four patches, our so-called B1 fields, from 2019–2021. These data consist of 20 individual DECam pointings targeting a progressively larger area of sky from 2019–2021, with significant overlap between years in order to enable the tracking of KBOs (see Figure 1). Table 1 gives the pointing for each night of data, as well as the detection efficiency parameters, calculated in Section 6. Note that in order to avoid double counting of single-epoch detections we only use a subset of our fields (indicated in Table 1) to derive the constraints in Sections 7 and 9.

A DEEP exposure sequence, designed with a cadence ideal for a technique called shift-and-stack (Tyson et al. 1992; Gladman & Kavelaars 1997; Allen et al. 2001; Bernstein et al. 2004; Holman et al. 2004; Parker & Kavelaars 2010; Heinze et al. 2015; Whidden et al. 2019), typically consists of ∼100 consecutive 120 s VR-band exposures of the target field. ¹¹ In the shift-and-stack approach, single-epoch images are shifted at the rate of a moving object (rather than at the sidereal rate) so that a moving object appears as a point source in the coadded stack.

There are two primary reasons why the shift-and-stack technique is preferable to long exposures for the discovery of moving objects. First, since the rate and direction of an object's motion are not known a priori, we are able to stack our images in a grid of velocity and direction that spans the space of possible KBO motions. The second benefit is the preservation of the signal-to-noise ratio (S/N) of moving objects. For stationary sources in astronomical images, S/N goes like t^1/2. A source is effectively stationary if its apparent position changes by less than the size of the point-spread function (PSF) over the course of the exposure. This sets an upper limit on the useful exposure time when searching for moving objects, which we will call ${t}_{\max }$. Given DECam's typical VR-band seeing of ∼09 and the typical KBO rate of motion of a few arcseconds per hour, our value of ${t}_{\max }$ is on the order of several minutes. When $t\geqslant {t}_{\max }$, the S/N of stationary sources continues to go like t^1/2, while smearing causes the S/N of moving sources to go like t⁰. Thus moving objects fade into the background while the S/Ns of background sources continue to grow. Because DECam's CCDs have negligible read noise, we lose no sensitivity by adding together many short images, and thus the S/Ns of the moving objects continue to increase like t^1/2.

In this section we describe how our images are processed in preparation for our shift-and-stack pipeline (described in Section 4). The following steps take place after the images have gone through preliminary reductions with the DECam community pipeline (Valdes et al. 2014).

3.1. Synthetic TNOs

3.2. Flux Calibration and Synthetic Source Injection

3.3. Difference Imaging

4.1. The Grid

To carry out our search, we developed a novel method of generating the shift-and-stack grid. We begin by computing the grid bounds. A unique grid is required for each field, as its exact shape depends on the epoch and sky position of a pointing. We must also select the range of topocentric distances (Δ) of interest. For the search described in this paper, we consider Δ ∈ [35, 1000] au. ¹⁶ With R.A., decl., and Δ fixed, an object's position vector in the topocentric frame, x _T , is uniquely determined. We then change the origin from the topocentric to the barycentric frame so that we have x _B . ¹⁷ After changing the origin to the barycentric frame, we assign a velocity v_bound such that an object at position x _B would be barely bound to the solar system (specifically, we use the speed appropriate for semimajor axis a = 200,000 au). We then uniformly sample a collection of velocity vectors v _B on the surface of the sphere of radius v_bound. With the velocity vectors specified, the state vectors are fully determined, allowing for the computation of the corresponding R.A. rates and decl. rates as observed in the topocentric frame. ¹⁸ We repeat this process at several discrete values of Δ, and then use the Python package alphashape to draw a concave hull bounding the computed rates. This hull encompasses the full region of physically possible sky motions for the distances of interest.

Zoom In Zoom Out Reset image size

Once we have computed the boundary of the concave hull for a given pointing, we choose a finite set of rates at which to stack our data. Toward this end, we employ a new method in which we fill the hull with a large sample of random points, and then use a K-means clustering algorithm to divide the region into N clusters. One can use dimensional analysis to determine that N ∝ A_hull/², where A_hull is the area of the hull and is the desired grid spacing to minimize the maximum source trailing (roughly determined by the PSF width divided by the duration of the exposure sequence). We take the centroid of each cluster as a grid point. When compared to a rectangular grid, this method allows us to use ∼20% fewer grid points, and simultaneously reduce both the mean and maximum distance of any given point in the hull to the nearest grid point. As a result, we achieve more even coverage of physically possible rates, while minimizing the computational cost and opportunities for false positives. Although the process is random by nature, a sufficient number of random samples makes it nearly deterministic. We show an example of a grid for a 4 hr exposure sequence with 1'' seeing in Figure 2.

Zoom In Zoom Out Reset image size

4.2. The Shift-and-Stack Procedure

After computing the shift-and-stack rates, we proceed with the stacking. We designate the first exposure in an exposure sequence as the reference image. We then compute the R.A. and decl. at the center of the reference image, and use the R.A. rate and decl. rate to compute the amount by which we have to shift each image to match with the reference image's center. In other words, we stack the pixels along the path taken by the center of the reference image for the given R.A. and decl. rates. We do not consider variations in the focal plane geometry across the chip, as the solid angle of a single chip is rather small (we thus assume that the chips are locally flat), and all of the images are SWarped (Bertin et al. 2002) for the difference imaging process. We perform separate stacks for both the weighted signal (i.e., the signal multiplied by the weight) and weighted images, and then obtain the full stack by dividing the stacked weighted signal image by the stacked and weighted image.

After each stack we use the Python package sep (Bertin & Arnouts 1996; Barbary 2016) to extract all sources with at least 3 pixels with values above 1.5σ. ¹⁹ Each stack is contaminated with of order a few thousand spurious sources mostly consisting of cosmic rays, dead pixels, oversaturated pixels close to bright stars, and residuals from poorly subtracted stars and galaxies. Since we use approximately 100 stacks per chip, the total number of spurious sources per chip is close to 10⁵. Because the vast majority of spurious sources are not PSF-like, we trained convolutional neural networks (CNNs) using tensorflow (Abadi et al. 2015) to reject them automatically. We trained one CNN on synthetic sources superimposed on a background made from DEEP difference images, and another on the autoscan training set that was used to train a random forest algorithm for background rejection in the DES supernova search (Goldstein et al. 2015). Both CNNs retain nearly all of the signal and fail to reject different types of background, thus enabling a significant performance gain by requiring a source to be classified as good by both CNNs. This procedure cuts the number of sources per stack by three orders of magnitude, down to a few hundred. After all of the stacks are completed for a given chip, we consider the complication that most objects are bright enough to be detected in adjacent stack rates. To eliminate this redundancy, we employ the DBSCAN (Ester et al. 1996) clustering algorithm in pixel space to group detections associated with the same object.

The grid spacing in the initial shift-and-stack is good enough for source detection, but is too coarse to provide the best values for the position and rate for a given source. To refine the parameters, we use a Markov Chain Monte Carlo (MCMC) approach in which we perform targeted stacks on our candidates to maximize the S/N. These targeted stacks are still restricted to the parameter space of bound orbits, but are now continuous in R.A. and decl. rates. This procedure enables us to obtain refined R.A. and decl. rates with uncertainties, while simultaneously optimizing the measured R.A. and decl. of the source. We use the uncertainties in the R.A. rate and decl. rate (which are typically about 01 hr⁻¹ in each dimension) to classify our detections probabilistically in Section 8. After refining the parameters of our detections, we discard all sources with rates slower than 3 pixels hr^–1 (079 hr⁻¹, or distance ≳ 150 au), as such slow rates tend to accumulate false positives due to subtraction artifacts much more quickly than faster rates. We feed our remaining candidates through a final CNN that reduces the number of sources per chip to ∼10. The images we show to the CNN are similar to the right panel in Figure 3, and contain more information than the cutouts we show to the first CNN. Good sources tend to show a characteristic radial pattern, while false sources do not.

Zoom In Zoom Out Reset image size

Once we have refined our sources' rates and positions, we compute their fluxes and flux uncertainties using sep. We use these values to calibrate the magnitudes of our detections against the known magnitudes of our implanted sources, as well as obtaining magnitude uncertainties.

Finally we do a reverse stack on our data, in which we repeat the above procedure with negated R.A. and decl. rates. Because no physical KBO would appear as a point source when stacked at these rates, all sources that result from this stack are false positives. This reverse stack enables the critical step of accounting for false positives in our detections. In Figure 4 we show differential histograms of the number of sources resulting from the forward and reverse shift-and-stack as a function of S/N, both before and after applying weights and various cuts. In Figure 5 we show a scatter plot of the R.A. and decl. rates of all detections from the forward and reverse shift-and-stack, prior to human vetting.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.3. Candidate Vetting

The final step in the discovery pipeline is the visual inspection of the grids that were labeled as good by the final CNN. This step includes the implanted synthetic sources and the sources from the reverse shift-and-stack. To do this visual inspection we developed a website where users vote yes, no, or maybe on a candidate and have their votes recorded in a database. The images are presented to the vetters in a blind manner, meaning that the vetter has no indication of whether an object is an implanted source, a source from the reverse shift-and-stack, or a true candidate. By blindly vetting all sources, we can reliably compute a voter's true and false positive rates for yes, no, and maybe votes. We require votes from three unique vetters for each object, and then combine the votes into a probability that a source is "real" using the following framework.

The odds (i.e., betting odds) of a source being good given a vote are

$$\begin{eqnarray}&&O(+|{\rm{vote}})=O(+)\displaystyle \frac{P\left({\rm{vote}}|+\right)}{P\left({\rm{vote}}|-\right)},\end{eqnarray} \tag{ 1 }$$

where O represents odds and P represents probability. The + symbol means that the object is truly a good source, and the − symbol means that the object is truly a bad source. We calculate the prior O(+) using the excess in the number of sources in the forward shift-and-stack over the number of sources in the reverse shift-and-stack. The fraction on the right-hand side of Equation (1) is the Bayes factor

$$\begin{eqnarray}&&B\equiv \displaystyle \frac{P(\mathrm{vote}| +)}{P(\mathrm{vote}| -)}.\end{eqnarray} \tag{ 2 }$$

The quantity P(vote∣+) is calculated as the probability of assigning a given vote to an implanted source. Similarly, the quantity P(vote∣−) is calculated as the probability assigning a given vote to a source from the reverse shift-and-stack.

Given multiple votes, we can simply update the information by taking the product of Bayes factors as

$$\begin{eqnarray}&&O(+|{\rm{votes}})=O(+)\displaystyle \prod _{i}{B}_{i},\end{eqnarray} \tag{ 3 }$$

where B_i is the Bayes factor of the ith voter. We can then convert the odds from Equation (3) into the probability that a source is real as

$$\begin{eqnarray}&&P(+| \mathrm{votes})=\displaystyle \frac{O(+| \mathrm{votes})}{1+O(+| \mathrm{votes})}.\end{eqnarray} \tag{ 4 }$$

We assign the values calculated by Equation (4) as a weight (w) for each of our detections. This treatment allows us to take a probabilistic approach in studying our detections in Sections 6–9. ²⁰

In principle, P(vote∣+), P(vote∣−), and O(+) can all vary with magnitude and rate of motion. In practice, however, we found that P(vote∣+) and P(vote∣−) did not vary much over the range of brightnesses and rates that we considered for our fits in Sections 7 and 9, so we considered them to be constants. Similarly, parameterizing O(+) had little effect on our detections' weights after human inspection, so we chose to treat it as constant. ²¹ See Table 2 for a tabulation of our vetters' Bayes factors, and see Figure 6 for correlations between the vetters' votes.

Zoom In Zoom Out Reset image size

Table 2. Anonymized Bayes Factors of Each of Our Vetters

Person	B(yes)	B(no)	B(maybe)
Vetter 0	17.50	0.12	0.49
Vetter 1	76.96	0.10	1.40
Vetter 2	849.50	0.10	3.17
Vetter 3	42.03	0.19	1.09
Vetter 4	23.54	0.35	1.76

Download table as:  ASCIITypeset image

In this section we qualitatively analyze our detections; we do more thorough quantitative analyses in Sections 6–9. In our 20 nights of data we detected a weighted sum of 2297.9 objects with weights greater than 0.01, corresponding to 2896 unique sources. We have elected to omit 3698 sources with weights less than 0.01, as such detections are rather unlikely to be real, and their omission does not change the results of our analysis. While the majority of our remaining sources have weights close to 1, there are some more ambiguous cases. We show a mosaic of all detections with weights ≥ 0.4 in Figure 7.

Zoom In Zoom Out Reset image size

Based on our distribution of observed magnitudes (shown in Figure 8), it appears that our detection efficiency begins to fall off at magnitudes fainter than r ∼ 26 (though we explore this in detail in Section 6).

Zoom In Zoom Out Reset image size

It is also informative to examine the sky moving rates of our detections, which we display in Figure 9. In this figure, the size of each marker is proportional to its weight. The most apparent feature here is a large population of objects moving at approximately 3'' hr⁻¹, mostly corresponding to CCs (see Section 8). Based on the increased density of small points at slow rates of motion, we note a propensity for slow-moving false positives to pass our CNNs, but to be later given low weights after human inspection (there is a noticeable overdensity of points at slow rates in Figure 5 which is absent in Figure 9). At faster rates, points with low weights appear to be evenly spread. The presence of such features provides further evidence that a great deal of care in avoiding false positives is required when making statistical use of single-night detections. When such detections can be linked to several epochs, single-night false positives will be less problematic.

Zoom In Zoom Out Reset image size

In order to make use of our detections, we must understand the efficiency with which we recover our implanted synthetic sources. We parameterize the detection efficiency as a function of apparent magnitude using a single hyperbolic tangent function, given by the following equation

$$\begin{eqnarray}&&\eta (m)=\displaystyle \frac{{\eta }_{0}}{2}\left(1-\tanh \left(\displaystyle \frac{m-{m}_{50}}{\sigma }\right)\right),\end{eqnarray} \tag{ 5 }$$

where η₀ is the peak detection efficiency, m₅₀ is the magnitude at which the detection efficiency drops to η₀/2, and σ is the width of the hyperbolic tangent function. ²² We weight each of the detections in our fit using Equation (4), and maximize the likelihood given by

$$\begin{eqnarray}&&\mathrm{ln}{ \mathcal L }(\theta )=\displaystyle \sum _{i}{w}_{i}\mathrm{ln}\left(\eta ({m}_{i}| \theta )\right)+(1-{w}_{i})\mathrm{ln}\left(1-\eta ({m}_{i}| \theta )\right),\end{eqnarray} \tag{ 6 }$$

where θ is the vector of function parameters, and undetected fakes receive w = 0. We display the best fit for each night in Figure 10, and list the fit parameters in Table 1.

Zoom In Zoom Out Reset image size

We use our characterized detections to compute the differential sky density Σ of the Kuiper Belt as a function of apparent magnitude. ²³ We reiterate that for this analysis we are using only single-epoch data (i.e., we have not linked these detections across multiple epochs), so we must treat each night as an independent survey. However, because our survey was designed to detect objects multiple times, our nights are not all statistically independent. As such, we selected a subset of our data consisting of the statistically independent set of nights {B1b 20201015, B1c 20201016, B1e 20201017, B1a 20201018, B1f 20201020, B1d 20201021} to do our fits. This subset offers the best combination of survey area and depth among all possible subsets of our data. ²⁴ Note that when fitting our data, we truncate our detection efficiency at m₅₀, and ignore all fainter detections.

For a given probability distribution Σ(m), the expected number of detections by a survey is given by

$$\begin{eqnarray}&&\bar{N}={\rm{\Omega }}\int \eta (m){\rm{\Sigma }}(m|\theta ){dm},\end{eqnarray} \tag{ 7 }$$

where Ω is the survey's areal coverage and η(m) is its detection efficiency. The variable θ is a vector of function parameters. Next, the probability of randomly drawing an object with magnitude m from Σ(m) is given by

$$\begin{eqnarray}&&P(m|\theta )=\int {\rm{\Sigma }}({m}^{{\rm{{\prime} }}}|\theta )\epsilon ({m}^{{\rm{{\prime} }}}|m,\delta m){{dm}}^{{\rm{{\prime} }}},\end{eqnarray} \tag{ 8 }$$

where is a functional representation of the magnitude uncertainty for which we have adopted a Gaussian centered at m, with a width of δ m. ²⁵

We calculate the underlying luminosity function (Σ) of the Kuiper Belt by maximizing the likelihood (${ \mathcal L }$) given by

$$\begin{eqnarray}&&{ \mathcal L }(\theta )=\displaystyle \prod _{k=1}^{n}{e}^{-{\bar{N}}_{k}}\displaystyle \prod _{j=1}^{{N}_{k}}P{\left({m}_{j,k}|\theta \right)}^{{w}_{j,k}},\end{eqnarray} \tag{ 9 }$$

(see, e.g., Loredo 2004; Fraser et al. 2014) where the index k runs through each night of data and the index j runs through the survey's detections. The value w_j,k denotes the weight of the jth object detected by the kth survey, as calculated by Equation (4). ²⁶

Several previous works have studied the form of the luminosity function of KBOs (Bernstein et al. 2004; Petit et al. 2006; Fraser et al. 2008; Fraser & Kavelaars 2009; Fuentes et al. 2009). Following the example of these studies, we fit our data with functional forms of varying complexity. We first try a single power law given by

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{single}}(m)={10}^{\alpha (m-m0)},\end{eqnarray} \tag{ 10 }$$

where m₀ is the magnitude at which the density of objects is one per square degree and α is the power-law slope. We next try a rolling power law given by

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{rolling}}(m)={{\rm{\Sigma }}}_{23}{10}^{{\alpha }_{1}(m-23)+{\alpha }_{2}{\left(m-23\right)}^{2}},\end{eqnarray} \tag{ 11 }$$

where Σ₂₃ is the number of objects with m_r = 23 per square degree, while α₁ and α₂ control the shape of the function. We fit a broken power law given by

$$\begin{eqnarray}\begin{array}{c}{{\rm{\Sigma }}}_{{\rm{broken}}}(m)=\left\{\begin{array}{cc}{10}^{\left.{\alpha }_{1}(m-{m}_{0}\right)} & \,m\lt {m}_{B},\\ {10}^{\left.{\alpha }_{2}(m-{m}_{0})+({\alpha }_{1}-{\alpha }_{2})({m}_{B}-{m}_{0}\right)} & \,m\geqslant {m}_{B},\end{array}\right.\end{array}\end{eqnarray} \tag{ 12 }$$

were m₀ is a normalization parameter, m_B is the magnitude at which the break occurs, and α₁ and α₂ are the bright and faint slopes, respectively. Finally we fit an exponentially tapered power law with the functional form

$$\begin{eqnarray}&&\begin{array}{c}{{\rm{\Sigma }}}_{{\rm{t}}{\rm{a}}{\rm{p}}{\rm{e}}{\rm{r}}}(m)=\exp \left[-{10}^{\left.\beta ({m}_{B}-m\right)}\right]{10}^{\left(\alpha -\beta )m-({m}_{0}\alpha \right)}\\ \,\times \text{}\left[\alpha {10}^{m\beta }+\beta {10}^{{m}_{B}\beta }\right],\end{array}\end{eqnarray} \tag{ 13 }$$

where α is the faint-end power-law slope, β is the strength of the exponential taper, m₀ is a normalization parameter, and m_B is the magnitude at which the exponential taper begins to dominate.

In each case, we obtain a best fit using an optimizer, and then estimate our uncertainties by running an MCMC. We then use a survey simulation technique to test the quality of our fits. We randomly sample a population of objects from our best-fit Σ(m∣θ), and then impose the detection criteria of our survey to simulate detections. After we simulate our detections, we construct a simulated empirical distribution, E(m). By repeating this process many times, we end up with an ensemble of simulated empirical distribution functions, {E_i (m)}. This ensemble is representative of the distribution of possible outcomes for our survey under the assumption that our best-fit Σ(m∣θ) is the underlying truth. Next we calculate the upper and lower limits for the central 95th percentile of {E_i (m)} as a function of m, which we call u(m) and l(m), respectively. We use these values to define a test statistic that measures the fraction of the interval m over which E_i (m) is an outlier at the 95% level. The statistic, which we call S, is given by

$$\begin{eqnarray}&&S=\int f(m)\,{dm},\end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray}\begin{array}{c}f(m)=\left\{\begin{array}{cc}0 & {\rm{if}}\,l(m)\lt E(m)\lt u(m),\\ 1 & {\rm{otherwise}}.\end{array}\right.\end{array}\end{eqnarray} \tag{ 15 }$$

We calculate S for each simulated survey to assemble a set {S_i }, and then for our actual detections, S_d . We then compute the quantile of S_d among {S_i }, which we call Q_outlier. Values of Q_outlier > 0.95 indicate a poor fit. Finally, we compute the Bayes information criterion (BIC) for each of our fits. In general, lower values of BIC indicate a preferred model. We demonstrate our results in Figure 11, and summarize them in Table 3.

Zoom In Zoom Out Reset image size

Table 3. Best-fit Parameters and Statistics for Each of the Distributions Tested on the Full KBO Sample

Functional Form	Best-fit Parameters	BIC	Qoutlier
Single Power Law	$\langle \alpha ,{m}_{0}\rangle =\langle {0.60}_{-0.02}^{+0.02},{23.31}_{-0.10}^{+0.09}\rangle$	5.2	0.99
Rolling Power Law	$\langle {\alpha }_{1},{\alpha }_{2},{{\rm{\Sigma }}}_{23}\rangle =\langle {0.89}_{-0.10}^{+0.11},-{0.07}_{-0.03}^{+0.02},{0.38}_{-0.09}^{+0.11}\rangle$	0	0.19
Broken Power Law	$\langle {\alpha }_{1},{\alpha }_{2},{m}_{0},{m}_{B}\rangle =\langle {0.90}_{-0.12}^{+0.16},{0.49}_{-0.05}^{+0.04},{23.59}_{-0.1a}^{+0.09},{24.59}_{-0.37}^{+0.26}\rangle$	2.1	0.18
Tapered Power Law	$\langle \alpha ,\beta ,{m}_{0},{m}_{B}\rangle =\langle {0.21}_{-0.08}^{+0.10},{0.14}_{-0.04}^{+0.05},{13.95}_{-11.40}^{+5.33},{28.74}_{-2.18}^{+2.71}\rangle$	7.9	0.83

Note. The BICs have been normalized such that the minimum value among the distributions is 0.

Download table as:  ASCIITypeset image

First we note that like several studies before us, we strongly rule out the single power law. The rolling, broken, and exponentially tapered power laws all provide acceptable fits. Judging by the BIC, the rolling power law is marginally preferable to the broken and exponentially tapered power laws, but we contend that a dearth of detections brightward of m_r = 24 limits the usefulness of such comparisons. There is also no overwhelming physical motivation for choosing between the distributions. While Kavelaars et al. (2021) showed that the absolute magnitude distribution of the CCs is well described by an exponentially tapered power law, there is no reason a priori that it should fit the full Kuiper Belt luminosity function. In fact, we might expect the full Kuiper Belt luminosity function to be fit poorly by an exponentially tapered power law, because it is likely a mix of multiple distributions. However, as we show in Section 8, our detections are dominated by CCs. It is therefore unsurprising that the exponentially tapered power law yields a reasonable fit. Given these considerations, we opt not to choose a preferred model, and instead claim for the time being that all three distributions provide acceptable fits to the DEEP detections.

While our single-night detections do not provide nearly enough information for secure dynamical classification, we can use the detections' rates of motion to make assumptions about their dynamical classes. Consider the space of allowable sky motions for bound orbits shown in Figure 12. As before, the large black outline contains the space of possible rates of motion for objects on bound orbits at topocentric distances of 35–1000 au. The blue region is the allowable parameter space for CCs (42.4 au < a < 47.7 au, 0 < e < 0.2, and 0° < i < 4°), and the red region is the allowable parameter space for a somewhat arbitrary definition of hot classicals (HCs) that we are only using for the purpose of demonstration (42.4 au < a < 47.7 au, 0 < e < 0.2, and 4° < i < 45°). Note that the exact shape and orientation of these regions vary as functions of sky position and epoch.

Zoom In Zoom Out Reset image size

In practice, there is some overlap in the rates of motion between different dynamical classes, which is accentuated by the uncertainties in our rate measurements. We use survey simulations to isolate the CCs in our detections. Using the OSSOS++ Kuiper Belt model, ²⁷ we calculate which objects our survey would have detected, and apply a smear in the R.A. and decl. rates consistent with our real detections. For each simulation we use a kernel density estimator to draw a region containing 95% of the CCs. We find that we can isolate a fairly pure sample (on average ∼70% purity; see Table 4) of CCs from our single-night data. This purity is somewhat serendipitous, as the DEEP B1 fields happen to be in a patch of sky that is relatively uncontaminated by objects in resonance with Neptune. We speculate that the uncertainties inherent to the distribution of objects in the OSSOS++ model are larger than the purely statistical uncertainties listed in Table 4. Nevertheless, these values should provide a realistic estimate of the purity of the CC samples that we attempt to isolate for the analysis in Section 9.

Finally, we obtain a distance estimate and uncertainty for each of our CC detections. These estimates rely on the fact that in the regime of the CCs, the relationship between an object's heliocentric distance and the inverse of its apparent rate of motion can be well approximated as linear. For each night of data, we project our CC population model (obtained from the OSSOS++ model) into the space of R.A. rate and decl. rate, and fit a line relating distance to the inverse of the apparent sky motion. We then use the resulting relationship to compute the heliocentric distance of each of our detections. By sampling from the covariance of our detections' rates, we obtain Gaussian distance uncertainties of at most 1–2 au, which end up being precise enough to fit an absolute magnitude distribution (Section 9).

Determining the absolute magnitude (H) of an object requires simultaneous knowledge of its apparent magnitude and its heliocentric distance. While our detections have reasonably well-constrained apparent magnitudes (∼±0.1 mag), their distances are not as well constrained from single-night detections a priori. However, as we found in Section 8, any true CCs are within a narrow range of r, with uncertainties of only 1–2 au, leading to uncertainty of only 0.1–0.2 mag in H.

For a given probability distribution Σ(H), the expected number of detections by a survey is given by

$$\begin{eqnarray}&&\bar{N}={\rm{\Omega }}\iint \eta (m){\rm{\Sigma }}(H(r,m)| \theta ){\rm{\Gamma }}(r)\,{dmdr},\end{eqnarray} \tag{ 16 }$$

where Ω is the survey's areal coverage and Γ is the underlying radial distribution of objects in the survey's field of view. Note that we use the OSSOS++ Kuiper Belt model (Kavelaars et al. 2021) to compute a kernel density for Γ. Since our fields are all observed near opposition, it is a good approximation to use

$$\begin{eqnarray}&&H=m-5{\mathrm{log}}_{10}(r(r-1)).\end{eqnarray} \tag{ 17 }$$

Next, the probability of randomly drawing an object with magnitude m and heliocentric distance r from the distribution Σ(H) is given by

$$\begin{eqnarray}&&\begin{array}{c}P(m,r|\theta )=\iint {\rm{\Sigma }}(H({r}^{{\rm{{\prime} }}},{m}^{{\rm{{\prime} }}})|\theta )\epsilon ({m}^{{\rm{{\prime} }}}|m,\delta m){\rm{\Gamma }}({r}^{{\rm{{\prime} }}})\\ \,\times \text{}\gamma ({r}^{{\rm{{\prime} }}}|r,\delta r){dm}^{{\rm{{\prime} }}}{dr}^{{\rm{{\prime} }}},\end{array}\end{eqnarray} \tag{ 18 }$$

where is a functional representation of the magnitude uncertainty and γ is a functional representation of the uncertainty in the heliocentric distance. Note that we have taken to be a Gaussian centered at m, with a width of δ m, and γ to be a Gaussian centered at r, with a width of δ r.

We again use a modified version of the method described by Fraser et al. (2014) in which we maximize the likelihood given by

$$\begin{eqnarray}&&{ \mathcal L }\left(\theta \right)=\displaystyle \prod _{k=1}^{n}{e}^{-{\bar{N}}_{k}}\displaystyle \prod _{j=1}^{{N}_{k}}P{\left({m}_{j,k},{r}_{j,k}|\theta \right)}^{{w}_{j,k}},\end{eqnarray} \tag{ 19 }$$

where the index k runs through the surveys (in our case individual nights of data), and the index j runs through the survey's detections.

We again fit single, rolling, broken, and exponentially tapered power laws. Note, however, that we have changed the definition of the rolling power law to

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{{\rm{rolling}}}(H)={{\rm{\Sigma }}}_{8}{10}^{{\alpha }_{1}(H-8)+{\alpha }_{2}{\left(H-8\right)}^{2}},\end{eqnarray} \tag{ 20 }$$

where Σ₈ is the number of objects per square degree per magnitude at H_r = 8.

For these fits, we must also account for the contamination in our sample from non-CCs. We do this by sampling from each night's detections that have rates consistent with being a CC, accepting each detection with a probability given by the CC fraction column in Table 4. Note that we omit the B1f field due to the projected low purity of the isolated CC sample. We then do the maximum likelihood fit, and compute our uncertainties as in the previous sections. We repeat this procedure 100 times, and then aggregate all of our all samples from all MCMCs, resulting in a combined posterior. We then calculate the best-fit parameters and uncertainties from the aggregation. For the fitting statistics, we now report the means of BIC and Q_outlier from our simulations as 〈BIC〉 and 〈Q_outlier〉, respectively. We show our fits in Figure 13, and present our results in tabular form in Table 5.

Zoom In Zoom Out Reset image size

Table 5. Best-fit Parameters and Statistics for Each of the Distributions Tested for the Absolute Magnitude Distribution of the CC Subsample of Our Detections

Functional Form	Best-fit Parameters	〈BIC〉	〈Qoutlier〉
Single Power Law	$\langle \alpha ,{H}_{0}\rangle =\langle {0.60}_{-0.04}^{+0.04},{7.43}_{-0.13}^{+0.12}\rangle$	4.0	0.86
Rolling Power Law	$\langle {\alpha }_{1},{\alpha }_{2},{{\rm{\Sigma }}}_{8}\rangle =\langle {0.77}_{-0.08}^{+0.10},-{0.10}_{-0.05}^{+0.04},{2.26}_{-0.33}^{+0.35}\rangle$	0	0.20
Broken Power Law	$\langle {\alpha }_{1},{\alpha }_{2},{H}_{0},{H}_{B}\rangle =\langle {1.02}_{-0.23}^{+0.48},{0.50}_{-0.09}^{+0.07},{7.57}_{-0.19}^{+0.14},{8.20}_{-0.67}^{+0.51}\rangle$	5.2	0.66
Tapered Power Law	$\langle \alpha ,\beta ,{H}_{0},{H}_{B}\rangle =\langle {0.26}_{-0.07}^{+0.08},{0.19}_{-0.06}^{+0.08},{2.22}_{-3.72}^{+2.42},{10.62}_{-1.14}^{+1.86}\rangle$	5.9	0.56

Note. The BIC values have each been rescaled such that the minimum value among the distributions is 0.

Download table as:  ASCIITypeset image

While the BIC favors the rolling power law, we do not believe that we can give preference to any distribution. This lack of constraining power is largely due to our dearth of detections at the bright end of the H distribution (we have no CC detections brighter than H_r ∼ 6.3). In contrast, Kavelaars et al. (2021) found that the OSSOS++ detections were inconsistent with rolling and broken power laws due to a taper at the bright end of the H distribution. The OSSOS++ model sensitivity was only possible due to the high purity of their sample and its well-characterized orbits. Our inability to rule out the rolling and broken power laws may be due to some combination of a relatively small survey area and contamination in our CC sample. On the other hand, the OSSOS++ sample only went as deep as H_r ∼ 8.3 (at least two magnitudes brighter than DEEP), requiring the assumption of a faint-end power law. In contrast, we constrain the faint end of the H distribution rather well. Our measured faint-end slope, $\alpha ={0.26}_{-0.07}^{+0.08}$, is consistent with SI simulations (Abod et al. 2019). In Figure 14 we show a comparison of our exponentially tapered power-law fit with that of Kavelaars et al. (2021). The two results are consistent (i.e., the OSSOS result is within our 95% confidence region), and both surveys are consistent with the result of Bernstein et al. (2004), which is the deepest KBO survey to date.

Zoom In Zoom Out Reset image size

The consistency of our detections with a rolling power law warrants careful consideration. First, we note that for appropriate parameters, a rolling power law is functionally equivalent to a Gaussian. All rolling power-law fits we have presented in this work satisfy these criteria, and can therefore be equally well represented as normal distributions. If the CCs are normally distributed in H, it would imply that they have a characteristic size, after which the H distribution turns over. While we currently have no reason to suspect that the size distribution turns over, our data cannot rule it out as a possibility. Furthermore, the only survey in the literature with data beyond our limit (Bernstein et al. 2004) cannot distinguish between the models, as it is consistent with extrapolations of both the rolling and exponentially tapered power laws. While more precise measurements of the CC H distribution from forthcoming DEEP data will help to bring the true distribution into clearer focus, an even deeper targeted CC survey may be necessary to distinguish between models.

Given our absolute magnitude distribution, we can calculate the total mass of CCs as

$$\begin{eqnarray}&&{M}_{{CC}}=\displaystyle \frac{1}{F}\int {\rm{\Sigma }}(H|\theta )M(H|p,\rho )\,{dH},\end{eqnarray} \tag{ 21 }$$

where F is an estimate of the average fraction of the total CC population per square degree in one of our fields (we calculate F ≈ 3.27 × 10⁻⁴ using the OSSOS++ model), and M is the mass of a body as a function of H given its geometric albedo, p, and mass density, ρ. If we assume that ρ and p are constant among the population of CCs, we can manipulate Equation (21) to pull the ρ and p dependence out of the integral, as

$$\begin{eqnarray}&&{M}_{{CC}}=\displaystyle \frac{\rho }{F}\displaystyle \frac{{1329}^{3}\pi }{6{p}^{3/2}}\int {\rm{\Sigma }}(H| \theta ){10}^{-0.6H}\,{dH}.\end{eqnarray} \tag{ 22 }$$

Note that the quantity

$$\begin{eqnarray*}&&\displaystyle \frac{{1329}^{3}\pi }{6{p}^{3/2}}\int {\rm{\Sigma }}(H|\theta ){10}^{-0.6H}\,{dH},\end{eqnarray*}$$

gives the volume of CCs per square degree in units of km³ deg⁻². We find that M_CC (H_r < 10.5) (i.e., the mass of the CCs up to our detection limit) is (at 95% CL)

$$\begin{eqnarray}&&\begin{array}{c}{M}_{CC}({H}_{r}\lt 10.5)={0.0014}_{-0.0003}^{+0.0008}{M}_{\oplus }\times {\left(\displaystyle \frac{p}{0.15}\right)}^{-3/2}\\ \times \text{}\left(\displaystyle \frac{\rho }{1.0\text{}{\rm{g}}\text{}{{\rm{c}}{\rm{m}}}^{-3}}\right)\times \text{}{\left(\displaystyle \frac{F}{3.27\times {10}^{-4}{\deg }^{-2}}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 23 }$$

If we extrapolate our fits our to H_r = 12, we find

$$\begin{eqnarray}&&\begin{array}{l}{M}_{{CC}}({H}_{r}\lt 12)={0.0017}_{-0.0004}^{+0.0010}{M}_{\oplus }\times {\left(\displaystyle \frac{p}{0.15}\right)}^{-3/2}\\ \ \ \times \ \left(\displaystyle \frac{\rho }{1.0\ {\rm{g}}\ {\mathrm{cm}}^{-3}}\right)\times \ {\left(\displaystyle \frac{F}{3.27\times {10}^{-4}\ {\deg }^{-2}}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 24 }$$

To facilitate a comparison, we can modify the form of the mass estimate reported by Bernstein et al. (2004) to

$$\begin{eqnarray}&&\begin{array}{c}{M}_{CC-B04}({H}_{r}\lesssim 14)={0.0016}_{-0.0003}^{+0.0003}{M}_{\oplus }\times {\left(\displaystyle \frac{p}{0.15}\right)}^{-3/2}\\ \,\times \text{}\left(\displaystyle \frac{\rho }{1.0\text{}{\rm{g}}\text{}{{\rm{c}}{\rm{m}}}^{-3}}\right)\times {\left(\displaystyle \frac{r}{43.8{\rm{a}}{\rm{u}}}\right)}^{6},\end{array}\end{eqnarray} \tag{ 25 }$$

where r is the average heliocentric distance of a CC. ²⁸ While the mass estimate from Bernstein et al. (2004) goes slightly deeper than our extrapolation, the extra mass is negligible, so the two mass estimates are in excellent agreement. Finally, we note that since our H distribution is consistent with that found by Kavelaars et al. (2021), our mass estimates must also be in agreement.

The only survey in the literature that is significantly deeper than the present work is that of Bernstein et al. (2004), which reached m₅₀ = 29.02 ²⁹ over a search area of 0.019 deg². Although the survey area is quite small, we can use it as a powerful lever arm to determine whether our fits remain valid down to H_r ∼ 12. To do so, we simulate the survey of Bernstein et al. (2004), using our Σ(H) fits (and their uncertainties) as the true underlying CC H distribution. For each H form, we simulate the survey 10³ times, and then ask whether the true number of detections made by the survey (N = 3) is commensurate with the suite of simulations. In particular, we calculate P(≤N), the probability that the survey would have made fewer than or exactly three detections. For the broken power law we find P(≤3) < 0.02, for the exponentially tapered power law we find P(≤3) = 0.16, and for the rolling power law we find P(≤3) = 0.65.

First, we note that our best-fit exponentially tapered and rolling power laws are both consistent with the Bernstein et al. (2004) detections. We are hesitant to rule out the broken power law because we find that dropping the m₅₀ of Bernstein et al. (2004) from r ∼ 29.02 to r ∼ 28.82 yields P(≤3) = 0.05 for the broken power law. While this result is still indicative of a marginal fit, it is not strong enough to rule out the broken power law altogether. Although the Bernstein et al. (2004) survey had sensitivity as faint as r ∼ 29.02, its faintest detection was r ∼ 28.23, so it is possible that the reported m₅₀ was overestimated. If, on the other hand, the reported efficiency of Bernstein et al. (2004) is accurate, it is possible that the H distribution of the CCs flattens out (and possibly turns over) somewhere between the limits of DEEP and Bernstein et al. (2004). While DEEP will not be able to resolve this issue, a very deep targeted CC survey should be able to answer the question (Stansberry et al. 2021).

In this paper we have presented our single-night detections from 20 nights of data in the DEEP B1 field. By using a shift-and-stack technique we were able to achieve an r-band depth of ∼26.2 over approximately 60 deg² of sky. Our data yielded 2297.9 single-epoch candidate detections, including 1849.8 detections fainter than m_r ∼ 25—the most detections fainter than m_r ∼ 25 ever reported in a single survey by more than an order of magnitude.

Our claim of fractional discoveries is a first for KBO science, and as such may seem peculiar. However, as we have shown in Section 4.2, a weighted treatment of our detections allows us to properly account for false positives. By accounting for false positives, we are able to make full use of the data from even our deepest nights, whose faintest detections cannot be recovered in another epoch. Additionally, our statistics remain reliable near the detection limit where false positives tend to accumulate.

Using 554.4 single-night detections from six unique DECam pointings, we computed the luminosity function of the Kuiper Belt as a whole (Section 7). Then using ∼280.0 CC detections from five unique DECam pointings we calculated the luminosity function of the CC population (Appendix 12), down to m_r ≳ 26.5. In both cases, we were able to confirm that a single power law is not suitable to describe the underlying distribution. We found that rolling, broken, and exponentially tapered power laws yielded acceptable fits, though low discrimination power at the bright end prevented us from giving overwhelming preference to any model.

The most significant scientific result of this work is a measurement of the absolute magnitude distribution of the CCs down to H_r ∼ 10.5. While a dearth of bright objects limits our constraining power at the bright end, our plethora of faint detections enable us to tightly constrain the faint end of the distribution. Our detections are consistent with an exponentially tapered power law with a faint-end slope $\alpha ={0.26}_{-0.07}^{+0.08}$. This faint-end slope is marginally shallower than previous measurements in the literature, but is consistent with simulations of planetesimal formation via the SI. This is of particular interest because the theory predicts a physically motivated size distribution (as opposed to other size distributions which have simply been engineered to adequately describe observations) of CCs that matches observations over the full range of sizes probed to date. However, we urge that these conclusions must be approached with caution. Of particular consequence to the conclusions from this work, the limited resolution in SI simulations leaves the theoretical faint-end slope somewhat uncertain (see Kavelaars et al. 2021 for an in-depth discussion of the outstanding problems with the SI as a complete theory of planetesimal formation). While the exponentially tapered power-law H distribution is a good fit to our CC detections, we cannot rule out rolling or broken power laws. Both the exponentially tapered and rolling power law distributions are consistent with the results of Bernstein et al. (2004), implying that the H distribution continues to flatten out beyond the DEEP limit. In contrast, our broken power law may be in tension with the the results of Bernstein et al. (2004). Interestingly, our rolling power law fit would imply a characteristic size for the CCs, beyond which the probability density rolls over.

Finally, we note some limitations of this work that can be improved upon in future studies. The most obvious limitation is our use of single-night detections. Although we developed robust new techniques to account for false positives, linked detections with well-constrained orbits are still preferable. Since we used single-night detections in this work, our selection of a CC subsample of our detections was rather rough, relying heavily on the OSSOS++ solar system model and the serendipity that Neptune was near these fields, meaning that they were necessarily relatively devoid of resonant objects. Linked orbits will enable proper dynamical classification, thus reducing the uncertainties in our studies of individual dynamical classes. In forthcoming works, we will analyze three additional fields similar to the B1 field. We will link our detections, yielding well-determined orbits for the smallest known KBOs. Our catalog of discoveries will enable studies of Kuiper Belt populations to unprecedented depth, providing deep insight to the formation and evolution of our planetary system.

We thank Gary Bernstein and an anonymous referee for insightful and engaged reviews that have improved the clarity and fidelity of this paper.This work is based in part on observations at Cerro Tololo Inter-American Observatory at NSF's NOIRLab (NOIRLab Prop. ID 2019A-0337; PI: D. Trilling), which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation.This work is supported by the National Aeronautics and Space Administration under grant Nos. NNX17AF21G and 80NSSC20K0670 issued through the SSO Planetary Astronomy Program and by the National Science Foundation under grant Nos. AST-2009096 and AST-2107800. This research was supported in part through computational resources and services provided by Advanced Research Computing at the University of Michigan, Ann Arbor. This work used the Extreme Science and Engineering Discovery Environment (Towns et al. 2014), which is supported by National Science Foundation grant number ACI-1548562. This work used the XSEDE Bridges GPU and Bridges-2 GPU-AI at the Pittsburgh Supercomputing Center through allocation TG-AST200009.K. J. Napier is supported by the Eric and Wendy Schmidt AI in Science Postdoctoral Fellowship, a Schmidt Futures program.H. Smotherman acknowledges support by NASA under grant No. 80NSSC21K1528 (FINESST). H. Smotherman, M. Jurić P. Bernardinelli, and C. Chandler acknowledge the support from the University of Washington College of Arts and Sciences, Department of Astronomy, and the DiRAC Institute. The DiRAC Institute is supported through generous gifts from the Charles and Lisa Simonyi Fund for Arts and Sciences and the Washington Research Foundation. M. Jurić wishes to acknowledge the support of the Washington Research Foundation Data Science Term Chair fund, and the University of Washington Provost's Initiative in Data-Intensive Discovery. This project used data obtained with the Dark Energy Camera (DECam), which was constructed by the Dark Energy Survey (DES) collaboration. Funding for the DES Projects has been provided by the US Department of Energy, the US National Science Foundation, the Ministry of Science and Education of Spain, the Science and Technology Facilities Council of the United Kingdom, the Higher Education Funding Council for England, the National Center for Supercomputing Applications at the University of Illinois at Urbana-Champaign, the Kavli Institute for Cosmological Physics at the University of Chicago, Center for Cosmology and Astro-Particle Physics at the Ohio State University, the Mitchell Institute for Fundamental Physics and Astronomy at Texas A&M University, Financiadora de Estudos e Projetos, Fundação Carlos Chagas Filho de Amparo á Pesquisa do Estado do Rio de Janeiro, Conselho Nacional de Desenvolvimento Científico e Tecnológico and the Ministério da Ciência, Tecnologia e Inovação, the Deutsche Forschungsgemeinschaft and the Collaborating Institutions in the Dark Energy Survey. The Collaborating Institutions are Argonne National Laboratory, the University of California at Santa Cruz, the University of California at Los Angeles, the University of Cambridge, Centro de Investigaciones Enérgeticas, Medioambientales y Tecnológicas-Madrid, the University of Chicago, University College London, the DES-Brazil Consortium, the University of Edinburgh, the Eidgenössische Technische Hochschule (ETH) Zürich, Fermi National Accelerator Laboratory, the University of Illinois at Urbana-Champaign, the Institut de Ciències de l'Espai (IEEC/CSIC), the Institut de Física d'Altes Energies, Lawrence Berkeley National Laboratory, the Ludwig-Maximilians Universität München and the associated Excellence Cluster Universe, the University of Michigan, NSF's NOIRLab, the University of Nottingham, the Ohio State University, the OzDES Membership Consortium, the University of Pennsylvania, the University of Portsmouth, SLAC National Accelerator Laboratory, Stanford University, the University of Sussex, and Texas A&M University.C.F. acknowledges support from the BASAL Centro de Astrofísica y Tecnologías Afines (CATA) ANID BASAL project FB210003 and from the Ministry of Economy, Development, and Tourism's Millennium Science Initiative through grant IC120009, awarded to The Millennium Institute of Astrophysics (MAS).

Here we calculate the luminosity function the subset of our detections that have rates consistent with being CCs (see Section 8). For these fits, we follow the same general sampling procedure as Section 9. We summarize our results in Table 6 and show the fits in Figure 15.

Zoom In Zoom Out Reset image size

Table 6. Best-fit Parameters and Statistics for Each of the Distributions Tested for the Luminosity Function of the CC Subsample of Our Detections

Functional Form	Best-fit Parameters	〈BIC〉	〈Qoutlier 〉
Single Power Law	$\langle \alpha ,{m}_{0}\rangle =\langle {0.59}_{-0.04}^{+0.04},{23.72}_{-0.13}^{+0.12}\rangle$	1.9	0.93
Rolling Power Law	$\langle {\alpha }_{1},{\alpha }_{2},{{\rm{\Sigma }}}_{23}\rangle =\langle {0.97}_{-0.17}^{+0.20},-{0.09}_{-0.04}^{+0.04},{0.19}_{-0.07}^{+0.10}\rangle$	0	0.21
Broken Power Law	$\langle {\alpha }_{1},{\alpha }_{2},{m}_{0},{m}_{B}\rangle =\langle {0.99}_{-0.21}^{+0.43},{0.48}_{-0.08}^{+0.07},{23.86}_{-0.18}^{+0.14},{24.55}_{-0.66}^{+0.43}\rangle$	3.6	0.55
Tapered Power Law	$\langle \alpha ,\beta ,{m}_{0},{m}_{B}\rangle =\langle {0.16}_{-0.05}^{+0.09},{0.17}_{-0.05}^{+0.07},{12.23}_{-8.05}^{+5.73},{28.02}_{-1.50}^{+2.47}\rangle$	7.9	0.80

Note. The BIC values have each been rescaled such that the minimum value among the distributions is 0.

Download table as:  ASCIITypeset image

Here a single power law provides a very marginal fit. As with the full set of KBOs, we cannot rule out the rolling, broken, or exponentially tapered power laws.