Observational Completion Limit of Minor Planets from the Asteroid Belt to Jupiter Trojans

The distributions of orbital elements and the physical characteristics of the minor planets between Mars and Jupiter are of great scientific interest because they are thought to retain important clues about the formation of the terrestrial planets and the origin and evolution of the solar system. As of 2019 September, more than 800,000 of these minor planets have been discovered, offering a statistically large sample to permit deep analyses for dynamical structures and their correlations with their physical properties. However, the discovered sample is subject to observational biases that must be understood for the proper interpretation of their demographics. Because these minor bodies are usually detected by their reflected sunlight, the observed sample suffers from the selection effect that brighter objects are more easily detectable than fainter objects. This bias manifests as an artificial peak in the distribution of the absolute magnitude, H, even when the intrinsic population of fainter/smaller objects grows rapidly with H. Therefore, in attempting demographic statistical analyses, it is often important to identify an observationally complete sample.

In the previous literature, the absolute magnitude at which a sample becomes observationally incomplete, ${H}_{\mathrm{lim}}$, has typically been determined using one of two methods. The first, and not uncommon method, is for a single value to be given without explicit justification, or to be sourced from an older work. The second approach, which we refer to as the exponential-fitting method, makes the assumption that the intrinsic population's size distribution is expected to increase as a power law as one considers objects of smaller and smaller diameters. Assuming constant albedo, H then follows an exponential function, and the H value at which a decline in the population of observed objects is seen relative to a fitted exponential function is taken to be ${H}_{\mathrm{lim}}$. Use of the exponential-fitting method requires several assumptions. First, it must be assumed that the asteroids follow a single power-law size distribution. Second, a range of H must be chosen to apply a linear regression. Third, the magnitude at which the deviation from linear is extreme enough to define ${H}_{\mathrm{lim}}$ must be chosen. The latter two require subjective decision making. Alternatively, the peak of the H distribution may be taken as the completeness limit (e.g., Ryan et al. 2015), a method that we adopt and expand in the present work.

Additionally, a common approach in previous studies has been to adopt a single value of ${H}_{\mathrm{lim}}$ for the entire asteroid belt or for large subsets, such as the inner, middle, and outer belts. For example, Gladman et al. (2009) and Malhotra & Wang (2017) adopt ${H}_{\mathrm{lim}}$ = 15, and Michtchenko et al. (2016) and Cambioni & Malhotra (2018) adopt ${H}_{\mathrm{lim}}$ = 15.5. Dermott et al. (2018) adopted ${H}_{\mathrm{lim}}$ = 16.5, while noting that the observational completion limit within the inner main belt varies from ${H}_{\mathrm{lim}}$ = 17 at 2.1 au to ${H}_{\mathrm{lim}}$ = 16.5 at 2.5 au.

To better assess the demographics of the asteroids at higher resolution in semimajor axis, we found a need to model the completion limit as a function of semimajor axis, ${H}_{\mathrm{lim}}$(a), and to ascertain that its value is contemporary with our sample population because ${H}_{\mathrm{lim}}$(a) evolves with time as observational surveys go deeper with technological advances.

We have developed a method to measure ${H}_{\mathrm{lim}}$(a) from observations, and to describe it with a simple, physically motivated, parametric model. We think that our approach will provide researchers with completion limit models for the asteroids and Jupiter Trojans that can be used to more accurately and consistently simulate observations, debias observations for statistical analyses, and model/synthesize minor-planet populations. This approach will lead to more reproducible results, and allow for the completion limit to be updated consistently over time as more minor bodies are discovered.

For our analysis, we use the Minor Planet Center Orbit Database (MPCORB). ¹ The MPCORB database provides, among other things, the osculating orbital elements of observed solar system minor bodies. As of 2019 September 29, this database contains data for 845,049 minor planets of which 729,567 have had multiple observations near opposition, and therefore have orbital elements and absolute magnitudes that are reliably known. In this work, we utilize two parameters from MPCORB: semimajor axis and absolute magnitude (a, H).

We investigated using proper elements from the AstDyS database (Knežević & Milani 2000), and decided that (a) because the present day orbital parameters are more appropriate for understanding observational bias, and (b) because AstDyS does not provide proper elements for all objects, the MPCORB database was more appropriate for this work.

3.1.  H_lim in Radial Bins

The upper panel of Figure 1 is a 2D density plot (i.e., heatmap or 2D histogram) showing the distribution of absolute magnitude (H) as a function of semimajor axis, a, in the range 1.5–5.5 au, for objects reported by the Minor Planet Center MPCORB database of 2019 September 29. In this figure, a general trend can be seen where the most populated H at a given distance decreases with a. We argue below that a specific functional form of this trend is expected due to the flux-limited completion limit of minor-planet observations.

Zoom In Zoom Out Reset image size

From this sample, we empirically determine the completion limit as a function of semimajor axis. We measure the empirical ${H}_{\mathrm{lim}}$ at each radial location by binning objects in both a and H and then identifying the H bin with the greatest number of objects, the center of which we then take to be the value of ${H}_{\mathrm{lim}}$ at that radial distance. (The measurement uncertainty on ${H}_{\mathrm{lim}}$ is discussed in Section 3.2.) This method is a generalization of the exponential-fitting method. However, our method requires fewer assumptions and fewer subjective decisions to be made in its application, making results more transportable and comparable between studies that implement it.

In Figure 2 we illustrate and compare our method of determining ${H}_{\mathrm{lim}}$ with the previously established method of using the exponential-fitting method. In the example, we perform both ${H}_{\mathrm{lim}}$ estimation methods at a single location (radial bin) containing 2041 observed objects with a ∈ (2.693, 2.695] (au).

Zoom In Zoom Out Reset image size

For this example, we find the peak of the distribution (the most populated H bin) to be at ${H}_{\mathrm{lim}}$ = 17.1 ± 0.29, where the uncertainty is the 1σ confidence interval. Regarding our bin widths in H, we note that some of the reported magnitudes in the MPC data have precision of 0.5 mag, the majority have precision of 0.1 mag, and some objects are reported to better precision. We used an H bin width of 0.25 mag, but found consistent results when using bin widths of 0.5–0.05. Users of our approach are encouraged to consider bin widths appropriate for their data set.

For the exponential-fitting method, one must determine what range of magnitude values are to be included in the fitting. In this example, we choose the range H ∈ [12.5, 16.75]. The adoption of a somewhat different range of values could be easily justified. We have chosen a range that, while not being an unreasonable range for the asteroid belt, shows the greatest disagreement between the two methods; this illustrates that the two methods result in similar ${H}_{\mathrm{lim}}$ for this sample. The data points used in fitting are shown in black. We perform Bayesian linear regression (Kelly 2007) on the data and the resulting best fit is shown as a blue line, with the corresponding confidence interval shown in lighter blue.

It is evident that one might reasonably choose "by eye" the location marked with a red "x" (at magnitude 16.75) as the completion limit because it is where the data deviate visibly from the black line (the linear best fit). As can be seen in the figure, this point lies within the 3σ confidence limit of the linear fit. In this example, the exponential-fitting method and our method both agree on ${H}_{\mathrm{lim}}$ ∼ 17, but the exponential-fitting method requires subjective choices.

We expect our method to result in values of ${H}_{\mathrm{lim}}$ that are systematically higher than those found by the exponential-fitting method. To get an idea of how much higher, we take 50 random semimajor axis bins (of width 0.002 au) between 2.6 and 2.8 au and measure the distance between the exponential fit and our estimated ${H}_{\mathrm{lim}}$. We find that the median distance is +0.24 mag with a standard deviation of 0.09 mag. This should be taken as a rather rough estimate because it is not clear that the exponential-fitting method represents the true completion limit value as the size distributions of minor planets are not consistent with single power laws (Cellino et al. 1991; Jedicke & Metcalfe 1998; Yoshida & Nakamura 2005) and its results are impacted by choices of which data is to be included and which excluded in the exponential fitting.

Our method assumes that the peak in the observed H distribution represents ${H}_{\mathrm{lim}}$. This assumption may not be valid if the distribution is multimodal. A local maximum in the intrinsic distribution could be confused with ${H}_{\mathrm{lim}}$ in the apparent distribution if it were significant enough to exceed the population of objects observed at the true ${H}_{\mathrm{lim}}$. While our method does not assume that minor-planet size–frequency distributions follow single power laws, and our method can be compatible with distributions that contain features or are not monotonically increasing, it is important to note that our method may return the wrong ${H}_{\mathrm{lim}}$ in the above case. This may be important to consider for some populations, e.g., the Kuiper Belt (Shankman et al. 2013; Singer et al. 2019).

Finally, it should be noted that our ${H}_{\mathrm{lim}}$ represents the magnitude at which observational completeness begins falling faster than the population of asteroids is increasing with decreasing size. In observational surveys for small bodies in the outer solar system, Bannister et al. (2016) show that detection efficiency can fall by 20% over 2 mag before falling off much more quickly.

3.2. Model Fitting

Our initial tests in fitting the measured ${H}_{\mathrm{lim}}$ as a function of a found that a linear modeling of ${H}_{\mathrm{lim}}$ could fit fairly well much of the data, but could not simultaneously fit the inner and outer edges of the main belt. This slight curve seen in Figure 1 (monotonically decreasing, concave up) motivated us to find a nonlinear function for our fitting, which we describe next.

We start with the premise that asteroids are observed at or near opposition with Earth-based telescopes detecting reflected sunlight. The flux of reflected sunlight from an asteroid of diameter D and geometric albedo p located at a heliocentric distance a is

$$\begin{eqnarray}&&F=\displaystyle \frac{{{pD}}^{2}{L}_{\odot }}{64\pi {a}^{2}{(a-1\mathrm{au})}^{2}},\end{eqnarray} \tag{ 1 }$$

where L_⊙ is solar luminosity. If F_lim is the smallest detectable flux, then in flux-limited observations, the diameter of the smallest asteroid that is detectable at heliocentric distance a is

$$\begin{eqnarray}&&{D}_{\mathrm{lim}}={\left(\displaystyle \frac{64\pi {F}_{\mathrm{lim}}}{{{pL}}_{\odot }}\right)}^{\displaystyle \frac{1}{2}}a(a-1\,\mathrm{au}).\end{eqnarray} \tag{ 2 }$$

To convert D_lim to H_lim, we use the standard definition of absolute magnitude for asteroids, H, as the apparent visual magnitude ² at distance 1 au from the Sun and 1 au from the observer (at zero phase angle), i.e.,

$$\begin{eqnarray}&&H={H}_{0}-2.5{\mathrm{log}}_{10}\displaystyle \frac{{F}_{1}}{{F}_{0}},\end{eqnarray} \tag{ 3 }$$

where H₀ is a standard reference magnitude and F₀ is a standard normalization flux; here F₁ is the flux of reflected sunlight from an asteroid imagined to be at heliocentric distance 1 au and geocentric distance 1 au (at zero phase angle from an Earth-based observer), i.e.,

$$\begin{eqnarray}&&{F}_{1}=\displaystyle \frac{{{pD}}^{2}{L}_{\odot }}{64\pi {\left(1\mathrm{au}\right)}^{4}}.\end{eqnarray} \tag{ 4 }$$

Setting D = D_lim in Equation (4) and using Equation (3), it follows that the absolute magnitude H_lim of the faintest detectable asteroid is

$$\begin{eqnarray}&&{H}_{\mathrm{lim}}=-5{\mathrm{log}}_{10}(a(a-1\,\mathrm{au}))+C,\end{eqnarray} \tag{ 5 }$$

where the constant C depends upon the unknown effective limiting flux of asteroid surveys to date,

$$\begin{eqnarray}&&C={H}_{0}-2.5{\mathrm{log}}_{10}\displaystyle \frac{{F}_{\mathrm{lim}}}{{F}_{0}}.\end{eqnarray} \tag{ 6 }$$

In deriving the model Equations (5)–(6), we have neglected phase angle effects; this can introduce some imperfection in our model for objects observed far from opposition, such as those at higher inclinations and at closer distances. Our model is necessarily approximate as it is meant to account for broad trends with heliocentric distance, from the Asteroid Belt to Jupiter Trojans; we do not consider near-Earth asteroids in our modeling or analysis.

We use Equation (5) for model fitting, with C as our free parameter. In our model fitting, we use semimajor axis as a proxy for the effective heliocentric distance of an asteroid; this is a reasonable approximation for our purpose because most asteroids have been observed at many points in their orbits so their semimajor axis is representative of their average heliocentric distance. We note that this is a strong assumption which appears to hold well for the main belt; we discuss it in more detail in Section 4.

Taking the values of ${H}_{\mathrm{lim}}$ (measured at each semimajor axis bin) as our data, we fit to it Equation (5) using a Markov chain Monte Carlo (MCMC) algorithm. We use the MCMC ensemble sampler implemented in the open-source code emcee (Goodman & Weare 2010; Foreman-Mackey et al. 2013). The posterior distribution of our free parameter (C, see Equation (6)) is sampled using 50 MCMC walkers. These walkers are initialized as a random uniform distribution ∼U(−50, 50). We find autocorrelation lengths (a measure of the number of steps needed for modeling convergence) for our parameter from 8 to 17 steps (depending on region, see Table 1). From this, and to be conservative, we throw away the first 500 steps (i.e., our "burn in" period) of each walker to end up with ∼10,000 samples of our posterior distribution. We use the posterior distribution generated by our MCMC fitting to find the 68% confidence interval around C.

Table 1. Model-fitting Results

Region	Inner	Outer	n	C	${H}_{\mathrm{lim}}$ a
	(au)	(au)		(mag)	(mag)
Hungarias	1.780	2.000	23527	19.57 ± 0.14	18.38
Main belt	2.120	3.250	783022	20.28 ± 0.032	17.01
Inner belt	2.120	2.500	246192	20.23 ± 0.058	17.76
Middle belt	2.500	2.960	328598	20.28 ± 0.044	17.01
Outer belt	2.960	3.250	208232	20.35 ± 0.047	16.28
Hildas	3.920	4.004	3656	21.04 ± 0.24	15.69
Trojans	5.095	5.319	6784	20.58 ± 0.13	13.88

Notes. Uncertainties reported for C are 68% confidence intervals.

^a Median ${H}_{\mathrm{lim}}$ for the given region.

Download table as:  ASCIITypeset image

To estimate the uncertainty of the computed ${H}_{\mathrm{lim}}$, consider a hypothetical bin with two objects; for such a bin, either H_min or H_max could be taken to be ${H}_{\mathrm{lim}}$. Whichever one is chosen, the uncertainty should take account of the other value because the probability is the same for each. We then adopt for the uncertainty the distance between the two magnitudes. This value respects and includes both H_min- and H_max within its bounds, and is connected to the bin population's distribution. Going beyond a sample of two objects, because we expect that the uncertainty would scale with sample size as in Poisson statistics, we divide by the square root of the number of asteroids in the bin. This approach estimates the measurement error as $({H}_{\max }-{H}_{\min })/\sqrt{n}$, where H_max and H_min are the largest and smallest observed H values within the considered region being modeled and n is the sample size in that region.

We have checked the effect of semimajor axis bin size on the measured value of ${H}_{\mathrm{lim}}$ and on the model-fitting results. Semimajor axis bin widths of 0.01, 0.05, and 0.002 au were tested, and we find that the resulting ${H}_{\mathrm{lim}}$ remained unchanged within the uncertainties. We adopted the bin width of 0.002 au for the model fitting, but use the larger bin width of 0.01 au in Figures 3 and 4 for legibility reasons.

Zoom In Zoom Out Reset image size

3.3. Results

Figure 3 shows the results of fitting our model to the entire main asteroid belt (2.12–3.25 au). We also fit several subregions independently, and we show the individual completion limit model fits for these regions in Figure 4. Table 1 lists these regions, their inner and outer edges as we have defined them, and the most likely value for the parameter C found from our fitting with the associated 68% credible region.

Zoom In Zoom Out Reset image size

As can be seen in Figure 3, the model is a good fit in the range of the main asteroid belt where the sample is most populous, 2.12 au < a < 3.25 au, although we do observe statistically significant deviations of a few individual semimajor axis bins. Across the main belt, our model finds a completion limit of ${H}_{\mathrm{lim}}$ = 18.4 at the innermost edge at a = 2.12 au, and ${H}_{\mathrm{lim}}$ = 15.9 at the outermost edge at a = 3.25 au. Regions within the main belt (inner, middle, outer), when fit individually, result in model parameters that agree closely (within 3σ) with the entire main belt fitting.

However, when we compare the regions beyond the main belt, we find C to vary by 1.5 mag. This is evident also when, as a fiducial reference, we extend the completion limit model for the main belt to other regions, we find that ${H}_{\mathrm{lim}}$ is over-predicted for the Hungaria objects, and underpredicted for the Hilda and Trojan objects. These differences likely arise from deviations from the main belt in the size, and eccentricity/inclination distributions of these populations (see Section 4). Differences in survey coverage of these three population groups likely also contribute to these deviations.

Knowing the completion limit of the minor planets' observational data is necessary to define an observationally complete sample before performing statistical analyses or testing theoretical models of the dynamical and physical evolution of asteroids. For example, before one can search for features in their size–frequency distribution, the trend imparted on the data from this observational bias must be removed.

The use of a single ${H}_{\mathrm{lim}}$ for the entire asteroid belt (or for subregions) may be appropriate for determining an observationally complete sample if a sufficiently conservative choice is made. However, for some applications, this may lead to a smaller sample size even when a significantly larger one is available, thereby "leaving data on the table" and limiting the statistical confidence of a study. Moreover, in a situation where one hopes to characterize a subtle trend throughout the asteroid belt (or a sub-region such as the inner asteroid belt), ${H}_{\mathrm{lim}}$ would need to be known not just as one value for the region, but throughout the region as a function of a. This requires estimating ${H}_{\mathrm{lim}}$ many times at high resolution in a.

Previous methods of determining ${H}_{\mathrm{lim}}$ (i.e., by eye, or by the exponential-fitting method) present several drawbacks. Determining ${H}_{\mathrm{lim}}$ by eye has an obvious problem with reproducibility; it is also impractical to do more than a few 10s of times for any extended portion of the asteroid belt. The exponential-fitting method also presents several issues to consider. To determine where the size distribution of a sample of objects deviates from a power-law function requires assumptions that make reproducibility difficult: first, a range of H values needs to be identified for inclusion in the exponential fitting, and second, one must choose the amount of deviation from an exponential which is large enough to identify the limiting magnitude. Additionally, the data samples being fit may not be well represented by single power-law size distributions (Cellino et al. 1991; Jedicke & Metcalfe 1998; Yoshida & Nakamura 2005).

We presented here our method of modeling ${H}_{\mathrm{lim}}$ which solves several of these problems.

First, by identifying the most populated H values throughout a range of a, our approach is agnostic about the underlying size distribution. Second, there are fewer decisions that need to be made when compared with the exponential-fitting method. This gives our approach the advantage that it is more easily reproducible (although the exponential-fitting method is reproducible in principle, the specific decisions in implementing it are usually not specified in the published literature).

Third, it is easily implemented in flexible ways: the ${H}_{\mathrm{lim}}$ model can be computed as a single value for an extended region, or one can compute ${H}_{\mathrm{lim}}$(a) as a function of semimajor axis with (user-determined) high resolution in a. In addition, because our approach is relatively hands-off, it scales to large data sets. This would be especially important in the coming decade as the the Vera C. Rubin Observatory's Legacy Survey of Space and Time increases the observational sample of asteroids by an order of magnitude (Ivezić et al. 2019).

As described in Section 3, we have measured the empirical completion limit of minor planets (as a function of a) within well populated regions between 1.78 and 5.30 au. We then use the results to fit to our physically motivated completion limit model (Equation (5)). Table 1 reports the regions, the number of objects in each region (n), the most likely value of our fit parameter C, and median ${H}_{\mathrm{lim}}$ for the given region, and Figure 4 plots the best-fit models.

As noted above, our model finds a completion limit of ${H}_{\mathrm{lim}}$ = 18.4 at the innermost edge of the main belt, at a = 2.12 au, and ${H}_{\mathrm{lim}}$ = 15.9 at the outermost edge at a = 3.25 au. This lowest value found is larger than many of the completion limit values adopted in recent studies, suggesting that those studies adopted somewhat conservative completion limits. (Note, however, that ${H}_{\mathrm{lim}}$ is steadily increasing over time as surveys become wider and deeper, consequently it is not straightforward to directly compare our ${H}_{\mathrm{lim}}$ with its value in previous studies.) If we were to make the conservative choice of adopting a single value, ${H}_{\mathrm{lim}}$ = 15.9 for the entire main belt, which is the approach taken in many previous studies, then our sample size would be 201,701. However, if we instead define our sample as H < ${H}_{\mathrm{lim}}$(a), the sample size is 465,837 objects, a factor of 2.31 larger. A larger sample size enables correspondingly greater statistical confidence.

The model fit, ${H}_{\mathrm{lim}}$(a) (Equation (5)), for the main belt accounts quite well for the broad trends with heliocentric distance in the main belt (see Figure 3). When we extend this model to regions outside of the main belt we find that the measured ${H}_{\mathrm{lim}}$(a) deviates significantly from the main belt model. In Figure 4 one can see that our main belt model over-predicts ${H}_{\mathrm{lim}}$(a) for Hungaria group and underpredicts it for the Hilda and Trojan groups. These discrepancies may result from differences in the distribution of sizes, eccentricities and/or inclinations among these regions; the causes of the discrepancies may be different for different regions. For the Hungarias, their smaller heliocentric distances and higher inclinations (16°–34°) are factors that reduce their detectability at low phase angles, qualitatively consistent with their observed discrepancy relative to the main belt model prediction in Figure 4. For the Hildas and Trojans, the explanation may lie in the significant dispersion of their eccentricities. It is conceivable that observational data of these objects is obtained more often when they are near their perihelion distance, so that our estimate of ${H}_{\mathrm{lim}}$ for these groups may be more representative of their perihelion distance rather than their semimajor axis. This explanation is qualitatively consistent with their observed discrepancy relative to the main belt model prediction in Figure 4: the ${H}_{\mathrm{lim}}$ values for these groups should be associated with systematically smaller heliocentric distance than their values of semimajor axis, a. The cause of these discrepancies is worth exploring quantitatively in the future.

It is our recommendation that future works consider the completion limit a function of semimajor axis, and we offer our parameterized model (and single values, for individual regions) for that purpose. H_lim—which will increase over time as surveys get deeper and wider—can be recalculated for any set of objects using the MPCORB database and our code available at https://github.com/equant/Asteroids.

We thank Jean-Marc Petit and the anonymous referee for comments that helped to improve this paper. This material is based upon work supported by the National Aeronautics and Space Administration under Agreement No. NNX15AD94G for the program "Earths in Other Solar Systems". The results reported herein benefitted from collaborations and/or information exchange within NASA's Nexus for Exoplanet System Science (NExSS) research coordination network sponsored by NASA's Science Mission Directorate. We acknowledge funding from the Marshall Foundation of Tucson, AZ.

This research has made use of data and/or services provided by the International Astronomical Union's Minor Planet Center.

Software: astropy, emcee (Foreman-Mackey et al. 2013), linmix.