DE-BIASED POPULATIONS OF KUIPER BELT OBJECTS FROM THE DEEP ECLIPTIC SURVEY

The discovery phase of the DES ended on 2005 May 11. A few KBOs were discovered during subsequent recovery efforts, but are not included in our sample since they were found under different search conditions. (For objects that were discovered during the main survey, we used the latest available orbital information to assign dynamical classes and calculate discovery probabilities.)

A total of 62,392 individual frames are in the DES database, corresponding to one of the eight CCD detectors on the Mosaic camera. Of these, 14,440 were not suitable for de-biasing because these frames were specifically targeted at the recovery of particular objects. Targeted recovery frames were not held to the same criteria for the number and brightness of stars, and have different statistical properties. An additional 1312 frames with solar elongation angle less than 140° were eliminated, because they have a different search efficiency than the bulk of the frames, which were taken near opposition. (Over twice as many objects were found on high solar elongation frames than on low-elongation frames.)

The 46,640 individual frames left were then paired up with matching frames. Each field was observed at least three times: a pair of images taken on the same night a few hours apart, and one image on another night. Sometimes additional frames were taken, for instance if weather interrupted the original observations. Only the pair of frames on the same night was used to search for objects. There were 18,668 pairs of frames that met these criteria. The overlap fraction for each pair is quite high, with a median positional offset of 0026. Each frame is 1065'' × 532'', and the largest mis-registration was less than 80'' in R.A. and less than 90'' in declination.

As of 2012 August, a total of 913 KBOs were listed in the DES internal database. Some objects have been excluded from this analysis because they were found outside the main survey, during recovery or other follow-up observations (101 objects). We also exclude eight objects found on low solar elongation fields (<140°). Several additional objects, typically those that were lost shortly after discovery, are excluded because of missing information: objects lack discovery field coordinates (36 objects), discovery distances (104 objects), or discovery magnitudes (1 object). This leaves 663 objects discovered during the main DES survey.

Of these objects, 478 objects had sufficient observations post-discovery to receive preliminary designations from the MPC (72%). Dynamical classification typically requires additional observations to achieve low enough errors, and 316 objects (48%) have been classified (see Section 2.3). Almost all (304) of these classifications are considered to be secure (quality 3 as defined in E05). For our analyses, we used the 304 securely classified DES objects.

We note that there is a potential for recovery bias in the particular subset of objects that were both discovered and classified by the DES. One could imagine that objects in unusual orbits could be missed more frequently than objects that were closer to their predicted location when follow up observations were taken a few months later. However, we do not think that this recovery bias introduces a major problem for the current work for the following reasons. (1) The single biggest reason for failure to recover an object was insufficient telescope time to re-observe an object within three months of discovery, after which observations were typically not attempted; this affected objects in all orbits alike. It was more likely that an object was lost because we did not have enough clear weather nights than because we looked in the wrong place and did not find it. (2) When the estimated orbital parameters for the lost objects were compared to those that were ultimately designated and classified, Elliot et al. (2005) found that the only main difference was that the lost objects were fainter (and thus would be harder to recover on an average frame). In this work we find that the rate of recovery was twice as high for objects with brighter discovery magnitudes (M_disc ⩽ 23) than for fainter objects, and we only do analyses using the better-recovered subset of brighter objects. Future work is planned to deal with the fainter objects, including the improved magnitude calibration, and we will also more fully explore the recovery bias at that time.

Over the years, several different schema for dynamically classifying KBOs have been proposed (e.g., Elliot et al. 2005; Lykawka & Mukai 2007; Gladman et al. 2008). In this paper, we use the current DES classification scheme, a modified version of the scheme laid out in E05 which incorporates the Scattering concept from Gladman et al. (2008). Since differences in the precise definitions of classes can affect the membership of objects and consequently the parameters derived for the class, we describe our classification scheme in full below.

To test for the dynamical class of an object, we use both the current observed orbital parameters as well as a 10 Myr integration of the orbit forward in time, and two additional 10 Myr integrations for the ±3σ values for a, e, and i. An object is said to be securely classified (quality 3) only if all three integrations agree on the same dynamical class. (In quality 2 objects, the nominal classification only agrees with one of the +3σ or −3σ integrations, while in quality 1 it agrees with neither.)

An object is tested for membership in each of the following dynamical classes, in order. (1) Resonant objects occupy a mean-motion resonance with Neptune (resonances up to ninth order are identified from libration in the resonant angle). (2) Centaurs have osculating perihelia that reach values less than the osculating semimajor axis of Neptune. (3) Scattered Near objects have a values that vary by more than an arbitrary amount, Δa ⩾ 0.02, over the 10 Myr integration (where Δa = (a_max − a_min)/a_mean). This is along the same lines as the Scattering class of Gladman et al. (2008), which uses a_max − a_min > 1.5 AU. An excitation statistic, $s=\sqrt{e^2 + \sin (i)^2}$, is applied to the remaining objects, which are sent to two additional classes: (4) Scattered Extended objects have s ⩾ 0.25, and (5) Classical objects have s < 0.25. Note that the exact criteria used can have an important impact on class membership: for instance, what might in other works be considered Hot Classical objects (low e, large i) with i > 15° would be classified as Scattered Extended under this scheme.

For the purposes of this paper, we have re-grouped a few objects in order to be more broadly consistent with the classes used by CFEPS and elsewhere. First, we restrict the population of Centaurs to only those with semimajor axes less than Neptune's (resulting in a population of smaller objects than found in any other class, due to the magnitude bias toward finding closer objects). Second, we group all of the dynamically excited objects with a_Nep ⩽ a ⩽ 80 AU into a single Scattered class, which contains objects from the Centaur, Scattered-Near, and Scattered-Extended classes. Ten objects with a = 82–739 AU are not considered in this analysis, because they are too sparse in parameter space to adequately define a dynamical class; for a brief discussion of such rare objects see Section 5.3.

To compute the probability that an object would have been detected by the DES at some time during the survey, we extend the methods described in E05 to include observational biases introduced by an object's H, a, e, and i values. We refer all inclinations to the Kuiper Belt plane (KBP), as derived in E05 (i = 174 ± 023, Ω = 992 ± 66), which is consistent with the invariable plane of the solar system at the 1σ level. Following the approach used in E05, we consider a set of N_O objects (j = 1, ..., N_O) discovered in a set of N_F search fields (k = 1, ..., N_F). For the jth object discovered by the survey, we denote its H magnitude, semimajor axis, eccentricity, and inclination by the symbols H_j, a_j, e_j, and i_j. We assume that the H magnitude of an object is not correlated with its orbital elements, and that angles that describe the orientation of the orbit in three-dimensional space average out over time. (For resonant objects, this assumption is not quite true, and is addressed in a separate bias factor discussed below.)

Each of the N_F fields of the survey data set is characterized by (1) the magnitude for which the detection efficiency has dropped to 1/2, which we denote by m_{1/2, k}, (2) the Kuiper Belt latitude and orientation, or tilt, of the frame (β_K and θ_K), used to calculate which orbits pass through a field, (3) the Kuiper Belt longitude of the frame (λ_K), and (4) the solid angle, Ω_k, subtended by the search field. Below we describe how biases related to all four of these characteristics have been removed.

• 1.  
  We assume that all search fields have the same maximum efficiency at bright magnitudes, _max, and the same characteristic range, σ_m, over which the detection efficiency drops from _max to 0; as in E05, we have fixed the values of σ_m = 0.58 and _max = 1. For the kth search field, we employ the functional form used by Trujillo et al. (2001) and in E05 (their Equation (19)) for the detection efficiency as a function of magnitude, (m, m_{1/2, k}):

  $$\begin{equation} \epsilon (m, m_{1/2,k}) = \frac{\epsilon _{{\rm max}}}{2} \left[ 1+\tanh { \left(\frac{ m_{1/2,k} - m}{\sigma _m} \right) } \right]. \end{equation} \tag{ 1 }$$

  Due to its elliptical orbit, the magnitude of the jth object varies from m_{min, j} to m_{max, j} (considering solely distance variations and ignoring physical changes). Between these extremes, the relative time that the object is at magnitude m is given by the probability distribution function p_j(m) for the jth object. Hence the likelihood, ξ_{mag, j, k}, for the jth object being discovered in the kth search field due solely to its magnitude (we shall consider other factors later) is given by

  $$\begin{equation} \xi _{{\rm mag},j,k} = \int _{m_{{\rm min},j}}^{m_{{\rm max},j}} p_j(m) \epsilon (m, m_{1/2,k}) dm. \end{equation} \tag{ 2 }$$

  Our strategy for evaluating the integral in Equation (2) involves transforming it from an integral over magnitude to an integral over heliocentric distance. To do this we make several simplifying assumptions. First we assume that all objects are discovered at opposition, and we set the Earth–Sun distance to 1 AU when computing the variation in the opposition magnitude of an object throughout its orbit. We define R_{d, j} as the heliocentric discovery distance for the jth object and m_{d, j} as its discovery magnitude. We neglect photometric variability due to changing phase angle and/or rotational light curve, since the mean amplitude of variability is only 0.1 mag, with only 15% of all KBOs varying by more than 0.15 mag (Thirouin et al. 2010). Then the magnitude of the jth object at a heliocentric distance R is given by the equation,

  $$\begin{eqnarray} m_j(R) &=& m_{d,j} + 5 \log _{10}(R/R_{d,j}) + 5 \log _{10} \nonumber\\ && \times [ (R-1 \,{\rm AU}) / (R_{d,j} - 1 \,{\rm AU}) ]. \end{eqnarray} \tag{ 3 }$$

  For the jth object, the minimum and maximum heliocentric distances are R_{min, j} and R_{max, j} respectively:

  $$\begin{equation} R_{{\rm min},j} = a_j (1 - e_j), \end{equation} \tag{ 4 }$$

  $$\begin{equation} R_{{\rm max},j} = a_j (1 + e_j). \end{equation} \tag{ 5 }$$

  To find the amount of time that the jth object spends at a given heliocentric distance, R, it is possible to directly numerically integrate a two-body orbit and then calculate the amount of time spent in each distance bin. However, this approach is far too slow for an equation that is called frequently, and there are no good analytical approximations that work near perihelion and aphelion, where the probability of finding an object peaks sharply. We thus calculated a grid of numerical integrations using e = 0.001 and e = 0.01–0.99 (Δ_e = 0.01). For each eccentricity, we fit a piecewise function composed of five line segments (f₁ to f₅) at predefined break points (p₁ to p₄). The coordinates for the slopes and offsets and break points of each line are stored for each eccentricity in a python executable file, and may be quickly called during numeric integrations. The function, which has been normalized over the allowed R range, is scaled to the appropriate a value when called. Tests indicate that the results are within ±15% of the values obtained using numeric integration (for objects with a < =80 AU, which is the cutoff in this paper for the Scattered class).

  $$\begin{equation} p_j(R) =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}0& R < R_{{\rm min},j} \\ f_1(e,R) & R_{{\rm min},j} \le R \le p_1 \\ f_2(e,R) & p_1 \le R \le p_2 \\ f_3(e,R) & p_2 \le R \le p_3 \\ f_4(e,R) & p_3 \le R \le p_4 \\ f_5(e,R) & p_4 \le R \le R_{{\rm max},j} \\ 0 & R > R_{{\rm max},j}. \end{array}\right.. \end{equation} \tag{ 6 }$$

  We are now in a position to rewrite Equation (2) in terms of an integral over heliocentric distance. Since our approximations reduce the apparent magnitude of an object with a given H magnitude to only the variation that depends on its semimajor axis and eccentricity, we emphasize this by using the symbol ξ_{j, k}(H, a, e) for the likelihood factor instead of ζ_{MAG, j, k} used in E05 (see the Glossary in their Appendix C). Substituting p_j(R)dR for p_j(m)dm in the integral and changing the limits of integration from magnitudes to the corresponding heliocentric distances, we find

  $$\begin{equation} \xi _{j,k}(H,a,e) = \int _{R_{{\rm min},j}}^{R_{{\rm max},j}} p_j(R) \epsilon (m_j(R), m_{1/2,k}) dR. \end{equation} \tag{ 7 }$$

  Equation (7) can be evaluated with substitutions from Equations (1), (3), and (6), using Equations (4) and (5) to set the limits of integration.
• 2.  
  We now consider the bias factor for inclination, ξ_{j, k}(i), using the same approach as in the inclination distribution analysis of Gulbis et al. (2010, hereafter G10). First we write an expression for the conditional probability of finding an object at latitude β with an inclination i, which is Equation (9) of G10:

  $$\begin{equation} p(\beta | i) =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{\cos \beta }{\pi \sqrt{\sin ^2 i - sin^2 \beta }}& \sin {i} > \vert \sin \beta \vert \\ 0 & \sin {i} \le \vert \sin \beta \vert \ \&\ i \ne \beta \ne 0\\ 1& i = \beta =0 \end{array}\right.. \end{equation} \tag{ 8 }$$

  (Note that in G10 and this work all inclinations and latitudes are relative to the KBP.)We define ξ_{j, k}(i) as the likelihood (based on the inclination alone) of detecting the jth object in the kth search field, which we determine by integrating the conditional probability for the inclination, i_j, over the range of KBP coordinates covered by the kth search field:

  $$\begin{equation} \xi _{j,k}(i) = \int _{\beta _{{\rm min},k}}^{\beta _{{\rm max},k}} \Delta \lambda (\beta ^{\prime }_k, \theta _k) p(\beta ^{\prime }_k | i_j) \,d\beta ^{\prime }_k, \end{equation} \tag{ 9 }$$

  where Δλ is the range of longitudes at each latitude on the frame (a geometric parameter derived in G10, their Equations (11)–(14)), and θ_K is the tilt of the frame with respect to the KBP.
• 3.  
  Resonant objects, by definition, oscillate around the resonant angle with a certain amplitude. Observationally, this has the effect that they are observed preferentially at some longitude relative to Neptune but not at others. Using the libration amplitudes found in a 10 Myr integration of the orbit of each object, we can determine a longitude bias, where the fields in the allowed longitude range, λ_{j, min} ⩽ λ ⩽ λ_{j, max}, are upweighted and all other longitudes are set to zero (note that for non-resonant objects ξ = 1):

  $$\begin{equation} \xi _{{\rm long},j} =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}360 / (\lambda _{j,{\rm max}} - \lambda _{j,{\rm min}})& \lambda _{j,{\rm min}} \le \lambda \le \lambda _{j,{\rm max}} \\ 0& {\rm otherwise} \end{array}\right.. \end{equation} \tag{ 10 }$$
• 4.  
  A minor correction is included for the solid angle of each set of fields that is available for object searches. The effective search area for a pair of frames varies slightly due to mis-registration (typically a few arcseconds) and obscuration by other objects on the field. Defining Ω_s as the solid angle of a full CCD and Ω_k as the net solid angle for the kth CCD, the solid angle component of the likelihood factor (following E05) is

  $$\begin{equation} \xi _{{\rm ang},k} =\Omega _k / \Omega _s. \end{equation} \tag{ 11 }$$

  Neglecting any correlation among a, e, i, and H, the combined likelihood, ζ_{j, k}(H, a, e, i), for detecting the jth object in the kth search field is the product of the four separate likelihood components in Equations (7), (9), (10), and (11):

  $$\begin{equation} \zeta _{j,k}(H,a,e,i) = \xi _{j,k} (H,a,e) \xi _{j,k} (i) \xi _{{\rm long},j} \xi _{{\rm ang},k}. \end{equation} \tag{ 12 }$$

  If we had observed all the search fields simultaneously, then the probability for having detected the jth object for the N_F search fields would be equal to the sum of the probabilities for each search field. This would be the case for a single night, or even a single lunation. However, we observe the search fields over a period of years, with the possibility of an object moving between search fields, so instead we take the efficiency of detecting the objects to be 1 minus the product of the likelihoods of having not detected it in each search field. We refer to this quantity as the detection probability of the survey for the jth object, q_det(H_j, a_j, e_j, i_j):

  $$\begin{equation} q_{{\rm det}}(H_j, a_j, e_j, i_j) = 1 - \prod _{k=1}^{N_F} (1 - \zeta _{j,k} (H,a,e,i)). \end{equation} \tag{ 13 }$$

  From the probability of detecting each individual object and other assumptions, we can characterize the distribution of unbiased orbital parameters and estimate the numbers of objects in each of the dynamical classes, as described in the next section.

For each class of objects, we want to determine the distribution functions for the parameters defining the class (a, e, i, and H). Here we present the functional forms used.

The simplest differential H magnitude distribution function, p_H(α, H), is an exponential function with exponent α. (We later present evidence for a double exponential, with two different power law indices following a break point, but since the break point lies near or outside of the range of the distribution functions for the classes we can fit, we only use a single exponent here.) All distribution functions are normalized over the full range of parameters considered, where c_H is a normalization constant so that the integral of the distribution function between H_min to H_max is equal to one. Thus the form of the H distribution is modeled using:

$$\begin{equation} p_H(\alpha ,H) =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}c_H 10^{\alpha H}& H_{{\rm min}} \le H \le H_{{\rm max}}, \\ 0& H < H_{{\rm min}}\ {\rm or}\ H_{{\rm max}} < H, \end{array}\right. \end{equation} \tag{ 14 }$$

where

$$\begin{equation} c_H = \frac{ \alpha \ln 10}{10^{\alpha H_{{\rm max}}} - 10^{\alpha H_{{\rm min}}} }. \end{equation} \tag{ 15 }$$

Next we consider the differential distribution function for semimajor axes. Between a_min and a_max for each class, we model this distribution function, p_a(a_o, Δa, a), as a Lorentzian with its central peak offset from zero by a_o and FWHM of 2Δ_a:

$$\begin{equation} p_a(a_o, \Delta a, a) =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{c_a}{1+\left[(a-a_o)/ \Delta _a \right] ^2}& a_{{\rm min}} \le a \le a_{{\rm max}}, \\ 0& a<a_{{\rm min}}\ {\rm or}\ a_{{\rm max}} < a, \end{array}\right. \end{equation} \tag{ 16 }$$

where the normalizing constant, c_a, is given by

$$\begin{equation} c_a = \frac{1}{ \Delta _a \left[ \arctan \left((a_{{\rm max}} - a_o) / \Delta _a \right) - \arctan \left((a_{{\rm min}} - a_o) / \Delta _a \right) \right]}. \end{equation} \tag{ 17 }$$

The distribution of eccentricities is similarly modeled as a single Lorentzian with parameters e_o, Δ_e, and c_e, using analogous versions of Equations (16) and (17).

The inclination distribution is also represented as a single Lorentzian, using i_o, Δ_i, and c_i. However, we note that other models have been used for the Classical distribution, notably a double Gaussian that represents two populations of objects, the "core" and "halo" populations (Brown 2001; Elliot et al. 2005; Gulbis et al. 2006), or alternately the "Hot" and "Cold" populations. To be specific, the first Gaussian component is the narrower one (core) and has a characteristic width (the standard deviation in the Gaussian expression) of Δ_i1. The second Gaussian component is broader (halo) with a characteristic width Δ_i2 > Δ_i1. We normalize each component separately, and b is a number between zero and one, representing the fraction of objects in the narrower Gaussian component. The inclination distribution function, p_i(b, Δ_i1, Δ_i2, i), can be written as:

$$\begin{eqnarray} && p_i(b,\Delta _{i1}, \Delta _{i2}, i) \nonumber\\ && \quad = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\frac{b c_{i1}}{\sqrt{2\pi \Delta _{i1}^2} } \exp \frac{-i^2}{2\Delta _{i1}^2} + \frac{(1-b) c_{i2}}{\sqrt{2\pi \Delta _{i2}^2} } \exp \frac{-i^2}{2\Delta _{i2}^2}& i_{{\rm min}} \le i \le i_{{\rm max}} \\ 0& i < i_{{\rm min}}\ {\rm or}\ i_{{\rm max}} < i \end{array}\right.. \nonumber\\ \end{eqnarray} \tag{ 18 }$$

The normalizing constants, c_i1 and c_i2, are given by:

$$\begin{eqnarray} c_{i1} &=& \frac{1}{I(\Delta _{i1}, i_{{\rm max}}) - I(\Delta _{i1}, i_{{\rm min}}) }\ {\rm and} \ \nonumber\\ c_{i2} &=& \frac{1}{I(\Delta _{i2}, i_{{\rm max}}) - I(\Delta _{i2}, i_{{\rm min}}) }, \end{eqnarray} \tag{ 19 }$$

where the integral, I(Δ, i) is given by the error function (erf):

$$\begin{equation} I(\Delta, i) = \frac{1}{2} {\rm erf }\left(\frac{i}{\sqrt{2} \Delta }\right). \end{equation} \tag{ 20 }$$

Using the current definition of the Classical class, however, this double Gaussian model is not required (see Section 4.2.2).

3.2.1. Range of Normalization

Using the functional forms defined in the previous section, we can now estimate the unbiased distribution of objects in a, e, i, and H, as well as the total number of objects per class. Both tasks are accomplished using an ML method.

3.3.1. Distribution Function Parameters

For each dynamical class, we fit for the parameters that best fit the unbiased distribution functions. We do this using a likelihood function, L, which is given by a product over every discovered object, j, of the probability density functions, raised to the inverse power of the detection probability for that object:

$$\begin{eqnarray} L &=& \prod _{j=1}^{N_0} [ p_H(\alpha ,H_j) p_a(a_o,\Delta _a,a_j) \nonumber\\ && \times p_e(e_o,\Delta _e,e_j) p_i(b, \Delta _{i1}, \Delta _{i2}, i_j) ]^{1/q_{{\rm det}}(H_j,a_j,e_j,i_j)}. \quad\quad \end{eqnarray} \tag{ 21 }$$

Each class has a different set of parameters, and each distribution has been normalized over the range considered (see Section 3.2.1) so that information from fields that did not detect objects is included in this fit.

Computationally, it is easier to maximize the natural logarithm of L, defined as M:

$$\begin{equation} M = \sum \limits_{j=1}^{N_0} \displaystyle\frac{ \ln \left[ p_H(\alpha ,H_j) p_a(a_o,\Delta _a,a_j) p_e(e_o,\Delta _e,e_j) p_i(b, \Delta _{i1}, \Delta _{i2}, i_j) \right]}{q_{{\rm det}}(H_j,a_j,e_j,i_j)}. \end{equation} \tag{ 22 }$$

To calculate the errors, we assume that the number of objects in the sample is great enough so that the formal errors on the ML estimates approximately follow a Gaussian distribution, using the method of Crooke et al. (1999). We define the matrix X as a j × n matrix, where there are j objects and n parameters (α, a_o, Δ_a, e_o, Δ_e, b, Δ_i1, and Δ_i2 if a double-Gaussian inclination model). Each element X_{j, n} is given by the expression:

$$\begin{equation} X_{j, n} = \left. \frac{\delta M}{\delta x_n} \right|_{H=H_j, a=a_j, e=e_j, i=i_j }. \end{equation} \tag{ 23 }$$

The correlation matrix, C, can be expressed in terms of X:

$$\begin{equation} \boldsymbol{C} = [\boldsymbol{X}^T\boldsymbol{X}^{-1}]. \end{equation} \tag{ 24 }$$

The variance of the parameter x_n is given by the appropriate diagonal element, C_{n, n}, of the correlation matrix, and the standard deviation is found by taking the square root:

$$\begin{equation} \sigma _n = \sqrt{C_{n,n}}. \end{equation} \tag{ 25 }$$

3.3.2. Inferred Number of Objects

Once we have determined the parameters for all of the unbiased distribution functions, we can estimate the number of objects in the class, again using ML. Within the range of H, a, e, and i values that the ML distribution function parameters are valid, we assume that the total population for a dynamical class is well described by the distribution functions. Given that assumption, we calculate the fraction of the total population detected by the DES, f, where 0 ⩽ f ⩽ 1, by integrating the distribution functions over the detection probability of the range of object parameters:

$$\begin{eqnarray} f &=& \int\!\!\int\!\!\int\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\int \limits ^{H_{{\rm max}}, a_{{\rm max}}, e_{{\rm max}}, i_{{\rm max}}}_{H_{{\rm min}}, a_{{\rm min}}, e_{{\rm min}}, i_{{\rm min}}} \;\;q_{{\rm det}}(H,a,e,i) p_H(\alpha , H) p_a(a_o,\Delta _a,a) \nonumber\\ && \times p_e(e_o, \Delta _e, e) p_i(b, \Delta _{i1}, \Delta _{i2}) \,dH \,da \,de \,di. \end{eqnarray} \tag{ 26 }$$

To turn the discovered fraction of objects into an estimated total number of actual objects, with error bars, we use the generalized binomial distribution, along with another ML function, where the observed number is N_o, the actual number of objects is $N_o^{\prime }$, and Γ is the gamma function:

$$\begin{eqnarray} L(N_o^{\prime }) &=& \frac{N_o^{\prime }!}{N_o! (N_o^{\prime } - N_o)!} f^{N_o} (1-f)^{N_o^{\prime }-N_o} \nonumber\\ &=& \frac{ \Gamma (N_o^{\prime } + 1) }{\Gamma (N_o+1) \Gamma (N_o^{\prime } - N_o + 1)} f^{N_o} (1-f)^{N_o^{\prime }-N_o} \quad \end{eqnarray} \tag{ 27 }$$

and thus

$$\begin{eqnarray} M(N_o^{\prime }) &=& \ln \left[ \frac{ \Gamma (N_o^{\prime } + 1) }{\Gamma (N_o+1) \Gamma (N_o^{\prime } - N_o + 1)} \right] \nonumber\\ && + N_o \ln f + {N_o^{\prime }-N_o} {\rm ln} (1-f). \end{eqnarray} \tag{ 28 }$$

The number of objects, $N_o^{\prime }$, is found by maximizing Equation (28), while the error is found by numerically calculating the total derivative:

$$\begin{equation} \sigma (N_o^{\prime }) = 1/ \sqrt{\left\vert \frac{d M(N_o^{\prime })}{d N_o^{\prime }}\right\vert }. \end{equation} \tag{ 29 }$$

Numerical implementation was done using Mathematica 8.0 and Python 2.7.3. Equation (13) gives the probability of detecting an object with an orbit described by a, e, and i, an absolute H magnitude, and which was discovered with magnitude m_d at distance R_d. Values for a, e, i, and H listed in the DES database in 2012 March were used in these calculations.

Although 304 DES objects have the highest-quality orbital classification, we restrict our sample to the 273 objects discovered with a VR magnitude of 23 or brighter. The nominal 50% detection efficiency of the DES determined by E05 was VR = 22.5, with VR = 23 corresponding to an efficiency of 15%, although the actual efficiency varies from frame to frame. The original photometric calibration of the DES (used in this and previous work) was based on the USNO-B1.0 catalog, which is however known to have magnitude uncertainties up to 0.5 mag. Recent recalibration of the DES (Buie et al. 2011) has been completed but not yet reapplied to the survey results. Future work using the revised photometry will have much more reliable magnitudes, particularly for objects detected at the faint edge of the DES's range. Until the re-calibrated data is available, we have calculated probabilities for all 304 objects, but restricted the de-biased class analysis to objects with m_{d, VR} ⩽ 23. Table 1 shows the number of objects the DES found in each dynamical class, as well as the breakdown for objects with m_d ⩽ 23 and additional H magnitude constraints used in subsequent fits to single exponential functions. Three non-resonant classes (Centaur, Classical, and Scattered) and 16 resonances (five with ⩾5 objects) have been identified for the following analysis.

The detection probabilities (Equation (13)) have been calculated for all objects with secure orbits, as shown in Column 2 in Table 2. The probabilities range from a maximum of 0.32 for the 3:2 resonant object 2002 GF₃₂ to a minimum of 1.5 × 10⁻⁷ for the faint Scattered object 2001FU₁₈₅. (Similar low probability objects are discussed in Section 5.3.)

Modelers work in physical sizes, while observers deal with apparent magnitudes, and properly comparing the results of the two requires some conversions and assumptions. The absolute, or H magnitude, which is easily derived from the observations, is the preferred variable where both types of researchers can meet. Very few KBOs have actual measured diameters (e.g., Elliot et al. 2010), so observers would have to make assumptions about the object's physical properties, such as the albedo (widely variable among KBOs; Stansberry et al. 2008) to produce sizes. Using the apparent magnitude, as many observers and modelers have unfortunately done, is not recommended, because it conflates detection biases and leads to errors. We recommend that modelers report sizes in kilometers or else convert to H magnitudes after noting the albedo assumed.

With the data in a uniform framework it is now possible to compare whether the size (or rather, H magnitude) distribution is the same across different dynamical classes. The implications of differences for the objects' collisional histories is discussed in Section 6.

4.2.1. Exponential Fits from Probability-weighted CDF

Only classes with a relatively large number of objects can use the ML method above to find distribution functions for all parameters. However, even classes with as few as five objects can be fit for an H magnitude distribution. A simple method to estimate the number of objects in a class is to convert the probability of detection into a predicted number of actual objects. A probability-weighted cumulative distribution of the estimated number of objects, N, brighter than a particular H magnitude, where there are N_o observed objects, is calculated by summing the inverse of the detection probability, q_det:

$$\begin{equation} N ({\le} H) = \sum _{j=1}^{N_o} \frac{1}{q_{{\rm det},j}}. \end{equation} \tag{ 30 }$$

For a given class x and some normalization scale C_x, the cumulative distribution function (CDF) can be fit with an exponential function of the form

$$\begin{equation} N_x ({\le} H) = C_x 10^{\alpha _x H}. \end{equation} \tag{ 31 }$$

A single exponent α has been shown to break down when the range of H magnitudes is large enough (Jewitt et al. 1998; Gladman et al. 2001; Bernstein et al. 2004; Fuentes & Holman 2008; Fraser & Kavelaars 2008; Shankman et al. 2013). This is evident in a turnover to a shallower slope, and can be mathematically described by introducing two slopes, α₁ brighter than a certain H magnitude, referred to as the break, H_b, and α₂ for fainter objects:

$$\begin{equation} N_x ({\le} H) =\left\lbrace \begin{array}{@{}l@{\quad }l@{}}C_x 10^{\alpha _1 H}& H \le H_b \\ C_x 10^{\alpha _1 H_b} 10^{\alpha _2(H - H_b)}& H_{b} < H \end{array}\right.. \end{equation} \tag{ 32 }$$

A plot of the probability-weighted CDF for each class is shown in Figure 1. We examined eight classes: Classical, Scattered, Centaur, 2:1, 3:2, 5:2, 5:3, and 7:4. We required that each class have N_o ⩾ 5 objects with m_d ⩽ 23, to avoid problems with sparse sampling at faint magnitudes (note that similar results are obtained using m_d ⩽ 23.5). Similarly, a few objects in the Scattered class have such large semimajor axes (hundreds of AU) that they distort any population fit due to sparse sampling and much lower probabilities of discovery; all objects with a > 80 AU are excluded from these analyses.

Zoom In Zoom Out Reset image size

A single exponential cannot explain the full distribution of objects. The location of the turnover to a shallower slope was estimated by examining the most populous class, the Classical objects. The single exponential fit with the lowest relative error has a break in slope at H ⩽ 7.2. A similar break is observed in the next two most populous classes, the 3:2 and Scattered objects, though the location of the break is not as well defined as in the Classical objects, since the 3:2 objects have more scatter at both high- and low-H objects, and both the Scattered and 3:2 have a likely coincidental lack of objects right around the putative break point. The slope for objects with H < 7.2 is α = 1.02 ± 0.01 from the Classical data; fits to several break points from 6.4 to 7.4 are shown in Figure 2.

Zoom In Zoom Out Reset image size

We note that break points are commonly seen at the faint edge of observational data, although the observed slope after the break varies by group. Fuentes & Holman (2008) observed objects spanning the break, which they placed at R = 24.3 with a slope of $\alpha _1=0.7^{+0.2}_{-0.1}$ before and $\alpha _2=0.3^{+0.2}_{-0.3}$ after. This corresponds to a diameter of D = 118(p/0.05)^−0.5 km (where albedo p = 0.05), or an H magnitude of about 8.5. The break seen in the DES exponential for Classical objects, at H = 7.2, could thus perhaps be a limiting magnitude issue, if the DES survey efficiency was being mis-estimated at the faintest, most difficult to detect magnitudes. To test the hypothesis that the limiting magnitude was causing an artificial turnover, we tightened our discovery magnitude threshold from VR = 23 to VR = 22.5. The same exponent was found, but the turnover in the Classical population appeared at H ⩽ 6.7 instead.

However, recent work by Fraser et al. (2014) also finds break points at a similarly low H, with their Cold population breaking at $H_{b,{\rm Cold}}=6.9^{+0.1}_{0.2}$ and the Hot population at $H_{b,{\rm Hot}}=7.7^{+1.0}_{0.5}$ in r'. Since their method is different (see Section 5.2) and their data is independent from what is presented here, it lends credence to the break point we have identified around H_b = 7.2 being real.

In addition, we have one population that clearly samples the size distribution after the break in slope: the Centaurs. With seven objects spanning a large range of H magnitudes after the purported break (7.5 < H < 11), the Centaurs are well fit by a single, shallower exponential, with α = 0.42 ± 0.2. (We note that similar slopes can be fit to the few Classical and 3:2 objects with H > 7.2, though the exact slopes are strongly subject to the exact choice of the artificial break point and have a much smaller lever arm spanning only 0.5–1 mag.)

Our proposed model for the size distribution of the Kuiper Belt is to assume that brighter objects follow the Classical H magnitude distribution for 5 ⩽ H ⩽ 7.2, and then shifts to a shallower slope for smaller objects as derived from the Centaur distribution. (There are hints in the 3:2 distribution that brighter objects may follow a shallower slope, but we have very few objects with H < 5, and none are Classical objects.) Combining the bright slope fit to Classical objects with H ⩽ 7.2 with the faint slope fit to Centaurs with H ⩾ 7.5, we can construct a double exponent: α₁ = 1.02 ± 0.01, α₂ = 0.42 ± 0.02, and H_b = 7.2 (see Table 3). Note that the exact location of the break point may be an artifact of the limiting magnitudes of our fields, as very few objects (especially Classical) were found fainter than M > 23, and so our CDFs suffer from incompleteness.

We check the compatibility of each of the eight classes of KBOs with this double exponent using a Kolmogorov–Smirnov (K-S) test, with one additional variable: although the values for α₁, α₂, and H_b are kept the same, the absolute number of objects is allowed to vary. The reasoning behind this choice is that objects in a resonance may have the same overall size distribution, but, due to varying efficiencies of resonance capture, may differ in absolute abundance. At the 5% level, all eight dynamical classes are compatible with this double exponent. The scaling constants required in the region following α₁, relative to the number of Classical objects, N_C = 1.6 × 10⁻⁴, are as follows: N_Scattered = 1.33N_C, N_{3: 2} = 0.27N_C, N_{2: 1} = 0.19N_C, N_{5: 2} = 0.13N_C, N_{7: 4} = 0.03N_C, and N_{5: 3} = 0.03N_C. The Centaurs, which follow α₂, have a different scale factor, N_Centaur = 0.02. In addition, the Hot and Cold populations are also consistent on a K-S test with the Classical distribution, with N_Hot = 2.56N_C and N_Cold = 0.72N_C. Visually, this agreement can be seen in Figure 1, which shows the scaled Classical fit with the data for each class.

We can estimate the total number of objects in each class less than a particular H magnitude. Table 4 shows the estimated total number of objects in each class from H = 4–9, using two models: (1) a scaled broken exponential with H_b = 7.2 for all dynamic classes except the Centaurs, and (2) a single exponential for the Centaurs. The errors are calculated by taking the difference between the nominal number and the numbers derived from the ±1σ errors. Note that the Centaurs are estimated using their own independent fit and errors.

Table 4. Estimated Number of Objects Less than H

⩽H	Classical	Scattered	3:2	2:1	5:2	7:4	5:3	Centaur
(1) α1 = 1.02, Hb = 7.2, α2 = 0.42, NC = 1.6 × 10−4								(2) αCen = 0.42
								NCen = 0.02
4	2 ± 0	2 ± 1	0 ± 1	0 ± 0	0 ± 0	0 ± 0	0 ± 0	1 ± 0
5	19 ± 2	26 ± 2	5 ± 1	4 ± 0	2 ± 1	1 ± 0	1 ± 0	2 ± 0
6	200 ± 30	270 ± 30	54 ± 7	38 ± 5	26 ± 3	7 ± 0	6 ± 0	5 ± 2
7	2100 ± 300	2800 ± 400	570 ± 80	400 ± 50	270 ± 40	69 ± 9	60 ± 8	13 ± 5
7.2	3400 ± 500	4500 ± 600	910 ± 90	640 ± 90	440 ± 60	110 ± 20	95 ± 15	16 ± 6
7.5	4500 ± 700	6000 ± 900	1200 ± 200	850 ± 130	580 ± 100	150 ± 20	130 ± 20	21 ± 9
8	7300 ± 1300	9600 ± 1400	2000 ± 300	1400 ± 200	940 ± 160	240 ± 40	200 ± 40	35 ± 15
9	19000 ± 5000	25000 ± 6000	5100 ± 1200	3600 ± 800	2500 ± 600	620 ± 150	530 ± 130	90 ± 50
10	49000 ± 15000	66000 ± 19000	13000 ± 4000	9300 ± 2700	6400 ± 1900	1600 ± 500	1400 ± 400	240 ± 130
Scale	1.0	1.33	0.27	0.19	0.13	0.03	0.03	NA

Download table as:  ASCIITypeset image

4.2.2. Estimating Class Distribution Parameters

For those objects with more members, we can more fully characterize the class properties. This approach uses an ML method to model the distribution functions for the parameters (a, e, i, and H) that define the class. These distributions are then used to calculate the detected fraction, f, from Equation (26), to find the total number of objects within our parameter range. A downside to ML is that more objects are needed in a class to fully characterize the four-dimensional parameter space. Only the Classical, Scattered, and 3:2 classes have sufficient objects to attempt an ML fit, and the Scattered fit is marginal. This is why we recommend using the values derived from the CDFs in Section 4.2.1 for comparison between many classes of objects.

We used somewhat expanded subsets of objects for the Scattered and 3:2 classes compared to the CDF analysis. For the ML analyses, a well-sampled parameter space is important, and we examined cutoffs for H magnitude from H = 6.5–8.5. For Classical objects, setting the break point between H = 6.7–7.2 does not substantially change the fitted value for α (e.g., α_{H ⩽ 6.7} = 1.05 ± 0.12), but including fainter objects does (α_{H ⩽ 7.5} = 0.60 ± 0.07). The distributions for a, e, and i, however, are resilient to the choice of H_b. Choosing a lower (brighter) cutoff means decreasing the sample size, and since there were few enough objects to start with, the cutoff was set at H = 7.5 for the Scattered and 3:2 classes, while Classical objects use the previously determined CDF cutoff of H = 7.2.

For all three classes we chose to model the a, e, and i distributions using a single Lorentzian (Equation (16)), following the preference for Lorentzians over Gaussians of E05 and G10. The Classical inclination distribution is also well fit by a double Gaussian (as in G10), but since the resulting distributions are not meaningfully different, the simpler model is preferred in this work. (Note: whether a single or double Gaussian is preferred depends critically on the exact classification scheme used. Four high-inclination objects that were deemed Classical objects in Gulbis et al. (2010) are classified as Scattered-Extended objects in this work, and that appears to explain why the double-Gaussian fit is not required here. Choice of classification scheme matters, particularly when comparing results!) The distributions and histograms of the biased and de-biased detections are shown in Figures 3 (Classical), 4 (3:2), and 5 (Scattered).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We list the best fit values for the Classical, 3:2, and Scattered distributions in Table 5. The values for α agree within their errors with the slopes found for the Classical objects in the CDF analysis. The distributions were normalized over the observed range of parameters for a and H,⁹ and from 0 ⩽ e ⩽ 0.95 and 0 ⩽ i ⩽ 170° (the full range has been truncated to avoid numerical errors near boundaries). We note that the spread in the eccentricity distribution for the Scattered objects is artificially narrow due to a small number of low-probability objects at e = 0.18 (an effect that would presumably go away with more objects).

Table 5. Maximum Likelihood Fitted Distribution Parameters^a

Param	Classicalb	3:2c	Scatteredd
α	1.01 ± 0.07	0.95 ± 0.16	1.04 ± 0.19
ao	43.94 ± 0.09	39.36 ± 0.05	45.87 ± 0.24
as	0.61 ± 0.09	0.2 ± 0.07	0.67 ± 0.28
eo	0.06 ± 0.01	0.24 ± 0.01	0.18 ± 0.01
es	0.05 ± 0.01	0.03 ± 0.01	0.01 ± 0.004
io	1.78 ± 0.1612	14.37 ± 1.52	21.99 ± 0.77
is	1.31 ± 0.24	4.14 ± 1.68	1.69 ± 0.56
H ⩽ 7.0	(2800)	(840)	(3100)
H ⩽ 7.2	4500 ± 400	(1300)	(5000)
H ⩽ 7.5	(9000)	2500 ± 450	10000 ± 2100

Notes. ^aParentheses indicate derived quantities. ^bObjects fit: quality 3, m_d ⩽ 23, and H ⩽ 7.2. ^cObjects fit: quality 3, m_d ⩽ 23, and 5 ⩽ H ⩽ 7.5. ^dObjects fit: quality 3, m_d ⩽ 23, and H ⩽ 7.5.

Download table as:  ASCIITypeset image

Using the calculated distribution functions, we evaluate Equation (26) along a grid of 4–7 points each in a–e–i–H space to find the detected fraction, f. The grid was constructed assuming that each grid object was discovered at its median distance from the Sun (i.e., half its time is spent closer to the Sun, and half its time is further out) and with a hypothetical discovery magnitude at that distance calculated using:

$$\begin{equation} m_d = H + 5 \log _{10}(R_d / 1) + 5 \log _{10}(R_d - 1/ 1). \end{equation} \tag{ 33 }$$

For the Classical objects, the fraction of total objects detected f_Classical = 0.027 for H ⩽ 7.2. The detected fraction of 3:2 objects is half as big, f_{3: 2} = 0.012, while the fraction of Scattered objects detected is ten times smaller, with f_Scattered = 0.0023 (both for H ⩽ 7.5).

Note that the assumed discovery distance in the grid has an effect on the fraction calculated. If we instead assumed that objects were always found at perihelion, instead of at the median distance, the detected fractions for the Scattered and 3:2 objects are eight times smaller. In the DES data, Classical objects were found, on average (median), close to the median distance (+1.6 AU from perihelion, and −0.6 AU from the median distance), while the 3:2 (+3.5 AU from perihelion and −4.3 AU from the median distance) and Scattered (+2.6 AU from perihelion and −5.9 AU from the median distance) objects were discovered somewhere in between. We chose to use the median distance rather than the perihelion (or the average distance of the discovered objects) to establish a minimum number of objects per class, and to avoid having the values of the grid depend on the properties of the objects found.

The corresponding total number of objects (from maximizing Equation (28) for $N_o^{\prime }$) is found to be N_<7.2 = 4500 ± 400 Classical objects, N_<7.5 = 2500  ±  450 3:2 objects, and N_<7.5 = 10, 000 ± 2100 Scattered objects (Table 5). Compared to the CDF fit values (Table 4), the ML method finds a larger total number of objects for all classes. At H = 7.2, there are 10% more Scattered, 30% more Classical, and 40% more 3:2 objects in the ML estimate than in the CDF estimate.

Comparisons between the results of the two largest independent surveys for KBOs, CFEPS and DES, need to keep in mind three differences in observations and analysis.

• 1.  
  Differences in observational strategies. The specific strategies used have been described in detail elsewhere (Millis et al. 2002; Elliot et al. 2005; Petit et al. 2011; Gladman et al. 2012), and may affect populations of observed and recovered objects. For instance, one should note that they used different filters (g' for CFEPS and VR for DES) and each cite results in H magnitudes relative to those filters, although the offset between them is small, roughly g' − VR = 0.1 mag.One area of potential concern relates to the recovery bias. The DES found a much larger number of objects initially, but due to limitations in telescope resources was unable to track and classify all of them, resulting in a smaller percentage recovered and classified compared to CFEPS. The breakdown of which types of DES objects were lost after initial discovery was reported in Elliot et al. (2005), and lost objects tended to be fainter (requiring better follow-up conditions) and faster-moving (allowing for less time to recover an object before errors accumulated) than objects that were successfully tracked. Future work will explore how important the recovery bias is. We note that the advent of all-sky surveys like LSST and Pan-STARRS will mean that many of these objects will eventually be recovered and linked back to the original observations (with high quality orbits due to the long observational baselines), mitigating any current recovery bias.
• 2.  
  Differences in class definitions. Although mean-motion resonances are straightforward to identify, the precise meaning of a "Classical" or "Scattered" object can differ significantly between groups. The properties of smaller groups of objects can be particularly affected. The CFEPS Scattering population (11 objects; Shankman et al. 2013) draws members from what the DES would classify as Scattered and Centaur groups, and as a result has a very different H magnitude distribution than the DES Scattered class. Also, as noted earlier, differences in definitions result in high-inclination CFEPS Classical objects being classified as DES Scattered objects. Special care should be taken when comparing results between different groups and models to make sure the same objects are being analyzed.
• 3.  
  Differences in statistical approach. CFEPS follows a recursive backward-modeling approach, making synthetic populations of objects (drawn from distribution functions for H–a–e–i as well as orbital longitude and libration amplitude for resonant objects). These synthetic populations are run through the CFEPS Survey Simulator to find which ones would have been detected. The properties of the initial distributions are adjusted until they match the recovered populations.By contrast, the main DES approach is a forward-modeling method, where we calculate the biases of discovering each object at a randomized point in its orbit, and turn that into a discovery probability. Two methods are used to extrapolate from individual probabilities to class populations. The H magnitude CDF method has each known object stand in for the (1/probability) objects that we didn't find with similar parameters. The ML method fits for distributions of H-a-e-i, assuming functional forms such as exponentials or Lorentzians, for each class. Although the ML method results in larger absolute numbers, it produces results consistent with the CDF fits at the 1–2σ level. Given the inability to use the ML method for most dynamical classes, we use the CDF results to compare with the size distributions found by other groups.

Despite the differences noted above, the CFEPS and DES results do agree on the slopes for most classes. CFEPS has carved the main belt into several sections ("hot," "stirred," and "kernel"), with slopes of 0.8, 1.2, and 1.2, respectively (Petit et al. 2011); when combined together, as in Figure 1, the overall CFEPS distribution is similar in slope to the DES Classical slope α_C = 1.02  ±  0.01 (CDF method), with similar absolute numbers as well. (The "hot" and "stirred"/"kernel" slopes, meanwhile, are consistent with the results when the full DES sample is divided into "Hot" and "Cold" populations; see Section 5.2.)

The absolute numbers of CFEPS objects are generally the same order of magnitude as the DES results, but not entirely consistent, with CFEPS generally citing lower abundances. After accounting for the slight differences between results quoted in H_g (Petit et al. 2011) and Gladman et al. (2012) and H_VR (this work), we find that the CFEPS Classical objects must be scaled upward by a factor of 1.7 to be consistent with the absolute number of DES objects (at the 5% level using a K-S test), while the 3:2 objects require a factor of 1.5. By contrast, the CFEPS Scattered objects are low by a factor of 14, although this discrepancy probably has more to do with the different definitions of "Scattered"/"Scattering." The other, smaller resonant classes (2:1, 7:4, 5:2, and 5:3), which have few objects in either survey (and hence higher uncertainty) are all consistent without any rescaling.

The ratios between classes also differ between the DES and CFEPS. Even ignoring the problematically defined Scattered classes, the ratio of 3:2 to Classical objects is nearly 4 for the DES but around 6 for CFEPS. The ratio for 2:1 objects, meanwhile, is either a factor of 5 (DES) or 21 (CFEPS). The 5:2 abundances, on the other hand, are approximately the same in both surveys (seven times as many Classicals as 5:2s). Since resonance ratios are used to distinguish between planet formation scenarios, future work is need to resolve these discrepancies.

There are many ways of subdividing populations of KBOs, and one long-standing division has been to divide the Classical belt between "Hot" and "Cold" objects, based on the orbital inclinations. We have noted before that the precise choice of class definition strongly impacts the values derived: for instance, several of the highest-inclination Classical objects in previous work (e.g., Gulbis et al. 2010) have been classified here as Scattered objects, removing the need for a double Gaussian (or "Hot Classical" population) inclination distribution. It is useful therefore to make comparisons wherever possible using the same classification scheme.

The definitions of "Hot" and "Cold" classes of Fraser et al. (2014) for all objects (not just Classicals) have simple inclination and distance cuts as follows: 38 ⩽ R ⩽ 48 and i < 5 are "Cold," while those from 30 ⩽ R ⩽ 150 and 5 < i < 90 are "Hot." We examined the DES sample using the same criteria (plus our restrictions to objects with M ⩽ 23 and a ⩽ 80 AU) and found that the best-fit slopes (H ⩽ 7.2) were α_Hot = 0.90 ± 0.02 and α_Cold = 1.32 ± 0.02. Despite using completely different objects and de-biasing methods, these numbers are very similar to Fraser's values: $\alpha _{{\rm Hot,Fraser}} = 0.87^{+0.07}_{-0.2}$ with a break at $H=7.7^{+1.0}_{-0.5}$, and $\alpha _{{\rm Cold,Fraser}} = 1.5^{+0.4}_{-0.2}$ with a break at $H=6.9^{+0.1}_{-0.2}$ (Fraser et al. 2014). The value after the break point for Fraser's Cold population (which is better constrained than the faint Hot population), is $\alpha _{2,{\rm Fraser}} = 0.38^{+0.05}_{-0.09}$, entirely consistent with our Centaur slope estimate of α = 0.42 ± 0.02.

The lowest calculated probabilities fall into two general categories: objects with discovery magnitudes well below our average detection efficiency, and very distant objects. Seven of the ten most improbable objects found by the DES have M_disc ⩾ 24, and the remaining three have a ⩾ 200. We have already noted that our current analysis is restricted to objects brighter than M_disc = 23 (future work using the magnitude recalibration may be able to extend that limit to fainter objects). We have also excluded the most distant objects (a ⩽ 80 AU) after finding that a few objects widely scattered in semimajor-axis space leads to huge distortions in the distributions derived for all parameters. Only nine objects that otherwise met our selection criteria (m_d ⩽ 23 and quality 3 orbits) were excluded, with 82 ⩽ a ⩽ 653 AU and classifications of Scattered-Near, Scattered-Extended, or Centaur. The CFEPS sample also has few objects in this range, with two scattering objects (2005 RH₅₂ and 2003 HB₅₇) around 150 AU (Petit et al. 2011) and one (2003 YQ₁₇₉) reported in a 5:1 resonance at 88 AU.

The lowest detection probability for any object is 10⁻⁷ for 2001FU₁₈₅, a very faint Scattered object (M_disc = 25.3). The lowest probability with M_disc ⩽ 23 is 10⁻⁵ for 87269, an extremely eccentric Centaur with a = 653 AU and e = 0.97. There are likely thousands of similarly far-flung objects, which are only detectable during a small fraction of their orbit, but characterizing their distribution parameters will require a larger sample and is beyond the scope of this work.

The existence of objects in mean-motion resonances is evidence that Neptune migrated outward from where it formed (Malhotra 1993). Many early numerical simulations (e.g., Hahn & Malhotra 2005) made the assumption that this outward migration happened more or less smoothly, with the giant planets maintaining generally low eccentricity, and that the population of objects that Neptune migrated through was dynamically "cold" (low e and i). However, when such models attempted to reproduce the observed structure of the Kuiper Belt, they had difficulties explaining the observed eccentricity and inclination distributions. For instance, a pre-existing dynamically "hot" (〈e = 0.1〉) disc allows for more efficient resonance capture into higher-order resonances, and Hahn & Malhotra (2005) found that objects would be expected in many "exotic" resonances that the DES has in fact found objects in, including the 5:2, 9:4, 7:3, 11:6, 12:5, 3:1, 7:2, and 4:1 (Table 1).

However, even with an initially excited disc, problems remain with smooth migration models. If the resonances were swept out of a Hot Classical disc, the inclination distribution of the 3:2 resonance should resemble that of the Classical population it formed from, and this is not the case (compare Figures 3 and 4). Instead, the 3:2 inclination distribution is much closer to the Scattered inclination distribution (Figure 5), indicating a different formation model (see Section 6.2). The smooth model is similarly argued against by the similarity in colors of the 3:2 and Scattered, but not Classical, populations (Sheppard 2012).

Comparing absolute numbers of objects, we find that the smooth migration model of Hahn & Malhotra (2005) overestimates the number of Classical objects by a factor of seven (130,000 predicted H ⩽ 9, compared to the DES estimate of 19,000), and finds the opposite ratios of 3:2 and 2:1 objects (twice as many 2:1, whereas we find 0.7 times as many) and Classical to Scattered (five times as many Classical, whereas we find 30% more Scattered). Other problems noted with smooth migration include the inability to explain objects with i > 15° or to populate the extended scattered disc (perihelia beyond 40 AU), where objects such as Sedna and several DES objects have been found. For all of these reasons, smooth migration cannot fully explain the evolution of the outer solar system.

Current evidence suggests that the early history of the solar system was much more chaotic than previously thought, as proposed and later refined in a series of papers collectively referred to as the "Nice Model" (Gomes et al. 2005; Tsiganis et al. 2005). In its current incarnation (Morbidelli & Crida 2007; Levison et al. 2011; Nesvorný & Morbidelli 2012), the model posits that a number of giant planets (up to six) were originally configured so that Jupiter and Saturn were near a major resonance, likely a mutual 3:2 resonance, and dramatic changes in planetary orbits ensued when the resonance was crossed. In simulations with more than four giants, one or two planets similar to Uranus or Neptune were ejected from the early solar system; a major resonance crossing has also been invoked to explain the late heavy bombardment of the inner solar system. All of these events left a dramatic mark on the Kuiper Belt.

In general, the Nice Model has been much more successful than previous attempts at explaining Kuiper Belt structure. Regardless of precisely how objects are partitioned between Classical and Scattered designations, there exists a bi-modal inclination distribution (Gulbis et al. 2010) which the Nice Model can reproduce. Additionally, the Nice Model predicted that some of the primordial population of objects was captured to make up Jupiter's Trojan population at the same time as the giant planets were sculpting the resonant and scattered populations of the Kuiper Belt (Morbidelli et al. 2005; Nesvorný et al. 2013), although observations at the time seemed to contradict the prediction that such objects would have similar size distribution slopes. A recent reanalysis of the Trojan size distribution by Fraser et al. (2014) has found that, in fact, the Trojan population has a similar two-slope fit to that work's results for the Hot population (α_{1, Trojan} = 1.0 ± 0.2, α_{2, Trojan} = 0.36 ± 0.01, and $H_{b,{\rm Trojan}} = 8.4^{+0.2}_{-0.1}$), which are also in agreement with this work (except for some difference in the break point, which is likely due in part to H magnitude conversions and differing average albedos).

Nonetheless, some important details still need to be worked out. As noted in Levison et al. (2008), the observed 3:2 distribution had higher eccentricities than the simulations could produce; as seen in Figure 4, the debiased eccentricity distribution is even more highly eccentric. The simulations also produced a Classical distribution that was too eccentric, although it has been noted by Batygin et al. (2011) that this may be less due to a failure of the model and more to the assumption that the space Neptune was pushing objects out into was previously unpopulated. Parker & Kavelaars (2010) find that a separate, primordial Cold population is also supported by the high wide binary fraction, which is too large for the mechanism of Levison et al. (2008).

Another interesting feature is that Cold Classical objects (defined by inclination) are less eccentric than the stability limit requires (Dawson & Murray-Clay 2012). For the least-excited (i < 2°) objects, there is a "wedge" in e−i space where no objects exist despite orbits being stable in that regime. Morbidelli et al. (2014) identified the area as just inside the 7:4 mean motion resonance. The wedge has been argued to support either a primordial cold belt of objects (Dawson & Murray-Clay 2012), or to be a natural result of resonance interactions during smooth migration at the tail end of Neptune's outward journey (Morbidelli et al. 2014). Clear evidence of this wedge is seen in the DES Classical objects in the same region (between the 5:3 and 7:4 mean motion resonances), where there are 43 objects. We can do CDF and ML fits as above, and find that the CDF slope of the population just inside the 7:4 resonance (42.3 ⩽ a ⩽ 43.6) agrees with the full Cold population, with α_{inside7: 4} = 1.25 ± 0.03, which is not surprising given the ML fit to the inclinations (i_o = 1.7, i_s = 1.7), while the eccentricities are indeed well below the stability limit of e = 0.1 (e_o = 0.043, e_s = 0.012). (We note that we find 10%–15% as many 7:4 as Classical objects in this region, or about half what the non-debiased results and model simulations of Morbidelli et al. (2014) both found; the difference is probably due to the slightly further distance, which adds 0.2 mag, as well as to the higher inclinations and eccentricities of 7:4 objects, all of which slightly bias against discovery even of this relatively close class.)

As more detailed model predictions of fine regions of phase space are made, it will become increasingly necessary that (1) all parties involved are using the same dynamical definitions, and (2) that observers make debiased data available and that modelers compare to debiased data wherever possible.

There has long been evidence that the size distribution of KBOs cannot follow a single power law. A break in the size distribution power law has been proposed for objects smaller than about 100 km diameter (Gladman et al. 2001; Bernstein et al. 2004; Kenyon & Bromley 2004), although the exact location varies by model. Destructive collisions have been evoked to explain the existence of a break at roughly the same size (e.g., 40 km; Pan & Sari 2005). The exact location of the size break, however, has not yet been observed because of the faint magnitudes of objects in this size regime: 100 km with an albedo of 0.05 roughly corresponds to H = 9.

Collisional evolution modeling by Fraser (2009) presented theoretical support for a depletion region from D = 20–40 km, with reduced collision rates for larger objects (D = 50–100 km), allowing the steep power law from the brighter Kuiper Belt to continue until that point. A "divot" is expected in the power law at around that point. Recently, there has been a claimed detection of that "divot" around H = 9 (Shankman et al. 2013), based on 11 Scattering objects discovered and de-biased by CFEPS using their Survey Simulator. These objects range in brightness from H = 7.1–10, with only two objects (L4k09 and L4v11) fainter than the supposed "divot" point (H = 9.5 and H = 10.0, respectively). The lack of additional objects is itself used as a constraint on the steepness of the slope.

Collisional grinding has come under fire recently, however, since it would have disrupted the many wide binaries found in the outer solar system (Parker & Kavelaars 2012; Fraser et al. 2014). Fraser et al. (2014) has compiled detections up to H = 9 and finds a much shallower turn-over at H ∼ 7, in line with what we have found in this independent work. After accounting for the albedo distribution, that paper concluded that the turn-over diameter is at D = 136 km ± 8. The picture that has emerged is of a Cold population (roughly the DES Classical group) that formed in situ around 40 AU, while the hot objects (DES Scattered plus resonances) formed around 15–35 AU and were scattered outward. The cold belt cannot have endured order-of-magnitude collisional evolution, and must have always had a surface density similar to the current density, while that of the hot objects must have been 10⁵ times greater. These five orders of magnitudes imply a growth time rate that is 10⁵ times slower in the outer solar system, and this creates intractable problems for object formation time scales and the development of a broken size distribution (Fraser et al. 2014). Another theoretical explanation for the existence of a break is required.

Locating the break observationally also is challenging. Most observations do not have numerous objects on both sides, although meta-analyses such as Fraser et al. (2014) do a good job of pulling together disparate observations. The only DES class with objects fainter than about H = 8.5 is the Centaurs, which approach closer than Neptune and are consequently often much brighter at discovery. Although Centaurs are now located inside the inner edge of the Kuiper Belt, they are thought to be the transitional stage between the scattered disc and Jupiter family comets (JFCs). This is supported by both dynamical evolution calculations (Levison & Duncan 1997; Tiscareno & Malhotra 2003) and the similarity of the Centaur inclination distribution (Gulbis et al. 2010) to both the Scattered objects and JFCs (and not to the main Kuiper Belt or resonant objects).

The slope observed for the Centaur objects, α₂ = 0.42  ±  0.02, is consistent with other populations of small objects reported in the literature, such as the JFCs (α = 0.49 ± 0.05; Solontoi et al. 2012). It is also the same power law as the faint end of the Jupiter Trojans as measured by Solontoi et al. (2012; α = 0.44  ±  0.05) and Fraser et al. (2014; α_{2, Trojan} = 0.36  ±  0.01). With seven DES objects between H = 7.5–11 and m_d ⩽ 23, we see no sign of the purported "divot" of Shankman et al. (2013), which most of our objects are fainter than (Figure 6). No sign of turnover is seen out to H = 11 (D = 40 km assuming albedo = 0.05).

Zoom In Zoom Out Reset image size