A Measurement of the Kuiper Belt's Mean Plane From Objects Classified By Machine Learning

Chiang & Choi (2008) posed the question:"If we could map, at fixed time, the instantaneous locations in three-dimensional space of all Kuiper Belt objects [KBOs], on what two-dimensional surface would the density of KBOs be greatest?" The authors demonstrated that this surface, also known as the Laplace surface, is given by the Laplace–Lagrange linear secular theory (Murray & Dermott 1999). This theory is based on the time-variable forcing arising from the planets' secular variations; consequently, the local normal on the Laplace surface varies only slowly with time; secular timescales for KBOs are much longer than ∼10⁴ yr. The Laplace surface for particles within the Kuiper Belt is not a flat plane because it has warps owing to secular resonances in certain localized regions of semimajor axes within the belt where the rate of orbit pole precession coincides with one of the inclination secular mode frequencies of the planets; at large semimajor axes the Laplace surface converges to the solar system's invariable plane. The invariable plane is the flat plane normal to the total orbital angular momentum of the solar system; this plane has an inclination of 158 and a longitude of ascending node of 10758 with respect to the J2000 ecliptic/equinox frame (Souami & Souchay 2012).

As a foil to the Laplace surface, previous studies have also considered the solar system's invariable plane as a candidate for the mean plane of the Kuiper Belt. The mean plane measured from observational data of KBOs is of interest for its potential to test dynamical models of the solar system and to reveal unmodeled perturbations when compared with theoretical predictions.

Previous measurements of the Kuiper Belt mean plane have produced inconsistent results. Brown & Pan (2004) reported a mean plane of (i₀, Ω₀) = (186, 816), where i₀ and Ω₀ are, respectively, the inclination and longitude of the ascending node of the plane in J2000 coordinates. With a measurement error of 037 in the pole position of the mean plane, this is consistent with the Laplace surface at semimajor axis a = 44 au (the median semimajor axis of their sample) but inconsistent (at more than 99.7% confidence) with the invariable plane of the solar system as well as the orbital planes of Neptune and Jupiter. Elliot et al. (2005) used five separate methods to measure the mean plane of the Kuiper Belt; they reported i₀ in the range 165–249 and Ω₀ in the range 974–1040, with a preferred value of (i₀, Ω₀) = (151 ± 026, 1000 ± 88). They rejected the Laplace surface at 99.7% confidence as the Kuiper Belt mean plane, but did not reject the invariable plane. Chiang & Choi (2008) measured the mean plane for two small samples of KBOs near a = 38 au and a = 43 au and concluded that they could not reject either the Laplace surface or the invariable plane for either sample. After a gap of more than a decade, the next measurement of the Kuiper Belt's mean plane was reported in Volk & Malhotra (2017; hereafter vm17) when the observed sample of KBOs with well-determined orbits (<5% semimajor axis uncertainty, observed over three or more oppositions) had grown to 931. These authors carefully identified resonant KBOs and discarded them from the mean plane measurement sample (because the assumptions of Laplace theory are violated for particles in mean motion resonances). vm17 reported an overall mean plane of (i₀, Ω₀) = (18, 77°) for the classical Kuiper Belt, i.e., the sample of nonresonant KBOs in the semimajor axis range 42–48 au. For smaller semimajor axis bins, their mean plane measurements in the semimajor axis range of 35–45 au were consistent with the Laplace surface, notably including the detection of the theoretically predicted warp near a = 40 au; save for that warp, the mean plane in this range was also consistent with the invariable plane. However, beyond a = 50 au, vm17 measured mean planes that were strongly warped away from the predicted Laplace surface and inconsistent with both the Laplace surface and the invariable plane at the 97%–99% confidence level. The most recent measurement of the Kuiper Belt's mean plane, by Van Laerhoven et al. (2019), was based on a sample of KBOs discovered in two specific observational surveys with well-characterized selection biases. These authors rejected the invariable plane and did not reject the Laplace surface for semimajor axes below 44.4 au; at larger semimajor axes both planes were accepted. The significant warp of the mean plane of the distant Kuiper Belt, beyond a ≈ 50, measured by vm17 but not detected by Van Laerhoven et al. (2019) is an open puzzle in the literature.

Since vm17 published their work, the sample of known nonresonant KBOs with well-determined orbits, semimajor axes between 35 and 150 au, and perihelia above the semimajor axis of Neptune has nearly doubled, growing from 931 to 1812. Figure 1 shows a scatter plot of the semimajor axes and ecliptic inclinations of these samples. In the present work, we revisit the measurement of the Kuiper Belt's mean plane with this larger sample with the goal to test the reproducibility of vm17's results. In Section 2, we briefly describe the theoretical background for the Laplace surface. In the next two sections, we describe the method for measuring the mean plane from observations with generally unknown biases (Section 3), and a Monte Carlo method for estimating the uncertainty of this measurement (Section 4). In Section 5, we describe the machine-learning tool used for identifying the sample of nonresonant KBOs. Our results are presented and discussed in Section 6.

Zoom In Zoom Out Reset image size

Though their measurements of the mean plane differ, researchers including vm17, Chiang & Choi (2008), and Brown & Pan (2004) all start with the premise that for KBOs with perihelia outside the orbit of Neptune, the theoretically expected mean plane is adequately described by Laplace theory, i.e., the linear secular perturbation theory as set forth in Murray & Dermott (1999). In this approximation, the planetary perturbations on a test particle are orbit-averaged and the perturbing potential (called the "disturbing function") is truncated to the first order in the planetary masses and to the second order in eccentricities and inclinations. It then follows from Lagrange's equations for the variation of the orbital elements that the test particle's semimajor axis is constant, while its orbital eccentricity and the orientation of its orbit vary quasiperiodically with time over secular timescales. Our interest here is in the orbit plane, which is described by two angular elements, the inclination, i, and the longitude of ascending node, Ω. The unit vector, $\hat{{\boldsymbol{h}}}$, normal to the orbit plane can be expressed in terms of these elements (Equation (5)). The Laplace theory expositions in the literature use the more convenient two-component vector, the so-called "inclination vector" defined as

$$\begin{eqnarray}&&(q,p)=\sin i\,(\cos {\rm{\Omega }},\sin {\rm{\Omega }}).\end{eqnarray} \tag{ 1 }$$

Laplace theory yields the solution of the linear secular equations for the test particle in which (q, p) can be written as a sum of two parts

$$\begin{eqnarray}&&(q,p)=({q}_{0},{p}_{0})+({q}_{1},{p}_{1}).\end{eqnarray} \tag{ 2 }$$

The first part, $({q}_{0},{p}_{0})=\sin {i}_{0}(\cos {{\rm{\Omega }}}_{0},\sin {{\rm{\Omega }}}_{0})$, is called the "forced inclination vector" and defines the Laplace surface; it is determined by the masses and semimajor axes of the planets and their instantaneous orbital planes, as well as the semimajor axis of the test particle. The second part, $({q}_{1},{p}_{1})=\sin {i}_{1}(\cos {{\rm{\Omega }}}_{1},\sin {{\rm{\Omega }}}_{1})$, called the "free inclination vector," is the remaining part of the total inclination vector; it is determined by the initial conditions. The Laplace surface changes over secular time, and its time variation is given by a superposition of the secular modes of the inclination vectors of the planets. The free inclination vector precesses around the Laplace surface at a constant angular rate: the free inclination i₁ remains constant and the free longitude of node Ω₁ circulates. As the orbit plane of a KBO precesses around its local Laplace surface, it follows that the mean plane of a large population of test particles in a small semimajor axis range (with dispersed orbit planes) is described by the local Laplace surface, as shown in simulations by Murray & Dermott (1999) and Chiang & Choi (2008).

For later reference, we note that the unit vector, ${\hat{{\boldsymbol{n}}}}_{0}$, normal to the local Laplace surface is related to the forced inclination vector as

$$\begin{eqnarray}&&{\hat{{\boldsymbol{n}}}}_{0}=\left({p}_{0},-{q}_{0},\sqrt{1-{p}_{0}^{2}-{q}_{0}^{2}}\right).\end{eqnarray} \tag{ 3 }$$

We briefly make a note about heliocentric versus barycentric orbital elements. The Laplace surface on 2023 February 20 for semimajor axes outside Neptune's orbit, shown in Figure 4, is computed according to the theory described in Murray & Dermott (1999), and using data for the J2000 ecliptic/equinox barycentric planetary elements retrieved from JPL Horizons and planetary masses from Standish (1995). Laplace theory as set forth in Murray & Dermott (1999) and other sources describes the Laplace surface in heliocentric coordinates, rather than barycentric coordinates. The portion of the disturbing function used to develop the theory depends on the eccentricity, inclination, longitude of node, longitude of pericenter, and semimajor axis. Only the semimajor axis is affected by the heliocentric-to-barycentric conversion. The heliocentric semimajor axis has short-period oscillations caused by the motions of planets interior to the Kuiper Belt; these average to zero over secular timescales. Barycentric elements, which lack those semimajor axis oscillations, are therefore more convenient and no less accurate than heliocentric coordinates for high-a Laplace surface calculations. As a practical matter, the Laplace surface for orbits outside Neptune's is nearly identical in heliocentric and barycentric coordinates, except in close proximity to the secular resonance locations near (barycentric) semimajor axis values of 35 and 40.3 au.

At first sight, it would seem reasonable simply to average the total inclination vectors of a group of KBOs to find their mean plane. Stated more carefully, one could average the unit vectors normal to the orbit planes of the KBOs to find the unit vector normal to their mean plane. However, without a complete or fair sample, the result of such a calculation will reflect the biases of the observational surveys that produced the sample (Brown & Pan 2004; Elliot et al. 2005; vm17). Observational surveys not only have brightness limits, but surveys of the outer solar system for KBOs also have limited and nonuniform coverage in ecliptic latitude and ecliptic longitude. The resulting KBO samples will yield average planes biased toward the observed regions, whatever the true population mean plane may be (vm17).

To mitigate the observational biases in catalogs of KBOs, Brown & Pan (2004) proposed that the unbiased mean plane of the Kuiper Belt could be identified with the plane of symmetry of the vectors, v _t , of KBOs' space velocities projected on the heliocentric celestial sphere; v _t can be called the sky-plane velocity vector. The instantaneous unit vector, ${\hat{{\boldsymbol{v}}}}_{t}$, along the direction of a KBO's sky-plane velocity vector can be computed from the KBO's observational data pertaining to its orbital plane and its position in the sky at some epoch, as follows. Denoting with $\hat{{\boldsymbol{h}}}$ the unit vector normal to the orbit plane, and with $\hat{{\boldsymbol{r}}}$ the unit vector along the radial direction to the KBO's heliocentric position, then

$$\begin{eqnarray}&&{\hat{{\boldsymbol{v}}}}_{t}=\hat{{\boldsymbol{h}}}\times \hat{{\boldsymbol{r}}}.\end{eqnarray} \tag{ 4 }$$

The unit vector, $\hat{{\boldsymbol{h}}}$ is found from a KBO's inclination and longitude of ascending node, i and Ω, and the unit vector $\hat{{\boldsymbol{r}}}$ is found from its J2000 ecliptic latitude and ecliptic longitude, β and λ, as

$$\begin{eqnarray}&&\hat{{\boldsymbol{h}}}=\left(\sin i\sin {\rm{\Omega }},-\sin i\cos {\rm{\Omega }},\cos i\right),\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\hat{{\boldsymbol{r}}}=\left(\cos \beta \cos \lambda ,\cos \beta \sin \lambda ,\sin \beta \right).\end{eqnarray} \tag{ 6 }$$

vm17 described a simpler implementation of the method of Brown & Pan (2004) for computing the mean plane of a KBO sample from the sky-plane velocity vectors. The sky-plane velocity of a KBO orbiting exactly on the mean plane will always be normal to the unit vector $\hat{{\boldsymbol{n}}}$ defining the normal to the mean plane, so that

$$\begin{eqnarray}&&\hat{{\boldsymbol{n}}}\cdot {\hat{{\boldsymbol{v}}}}_{t}=0.\end{eqnarray} \tag{ 7 }$$

The mean plane can then be computed as the plane whose orbit normal is most normal to most of the KBO sky-plane velocities, i.e., the plane that minimizes

$$\begin{eqnarray}&&S=\displaystyle \sum _{i=1}^{N}| \hat{{\boldsymbol{n}}}\cdot {\hat{{\boldsymbol{v}}}}_{t,i}| .\end{eqnarray} \tag{ 8 }$$

The function S may be minimized by a grid search in the (q, p) parameter space for a computationally simple, though grid-dependent, method of finding the mean plane. That is, we search on a 700 × 700 grid with ${q}_{\min }={p}_{\min }=-0.35$ and ${q}_{\max }={p}_{\max }=+0.35$. The precision of the mean plane identification is limited by the grid spacing. With our choice of grid, the precision of the pole position of the computed mean plane is ≲01. This is small compared to the confidence intervals that are computed later. The mean plane (q₀, p₀) is trivially converted to an inclination and longitude of node (i₀, Ω₀).

Brown & Pan (2004) and vm17 showed, with simulations of biased synthetic samples, that this method of estimating the midplane by using sky-plane velocity vectors substantially mitigates the systematic error and recovers the true mean pole more reliably than simply averaging the orbit pole unit vectors. For example, when this method is applied to a synthetic sample confined to only a small patch of the sky, it recovers the true mean pole whereas the average of the unit pole vectors is systematically offset from the true mean pole (see, e.g., Figure 8 in vm17). A present disadvantage of this method is that the confidence limits of the measured mean plane are not straightforward to compute (discussed in the next section). The performance of this method and its limitations and reliability with different survey designs merit further exploration.

The method of measuring the mean plane as the plane of symmetry of the sky-plane motion vectors does not render itself to a straightforward method for computing the measurement uncertainty. We follow the method of Brown & Pan (2004) and vm17, who adopted a heuristic approach which assumes a statistical model of the intrinsic inclination distribution rather than fitting for it, and constructs Monte Carlo simulations to generate synthetic samples which approximately account for the biases in the observed sample. In brief, this approach is as follows. For a given observational sample of N_K nonresonant KBOs, we generate a new list of N_K simulated objects. The simulated objects are assigned orbital elements drawn randomly from a statistical model of their intrinsic distribution. Each random draw is accepted or rejected by subjecting it to a comparison with the properties of the observed sample; this approximately accounts for the selection biases in the observational data. With a simulated sample in hand, we compute its mean plane, and repeat N_R times with additional simulated samples. In this way, we generate the statistics of the mean plane. We use the distribution of these N_R Monte-Carlo-sampled mean planes to estimate the 68.2%, 95.4%, and 99.7% confidence error bars on the measured mean plane. A more detailed description of these Monte Carlo simulations is given below.

In each of the N_R simulated samples, the ith simulated object is found by first randomly selecting the jth real object from the observational KBO sample, where j is an integer randomly selected from the discrete uniform distribution on [1, N_K ]. (How the objects are ordered, 1 − N_K , does not matter.) We then choose its semimajor axis a_i randomly from the continuous uniform distribution on [0.99a_j , 1.01a_j ), and we choose its eccentricity e_i randomly from the continuous uniform distribution on [0.95e_j , 1.05e_j ), so that the model population has approximately the same semimajor axis and eccentricity distribution as the observed population.

The most important assumption is a model for the intrinsic distribution of free inclinations. Our current best understanding of the intrinsic inclination distribution in the Kuiper Belt is found in estimates obtained from a few well-characterized observational campaigns, such as the Canada-France Ecliptic Plane Survey (CFEPS; Petit et al. 2011, 2017) and the Outer solar systems Origins Survey (OSSOS; Bannister et al. 2018). The data from these surveys are more effectively debiased with the use of survey simulators. They are then fit to physically motivated models for the intrinsic inclination distribution function. These models use the Rayleigh distribution function, or a combination of two Rayleigh distributions, one with a narrow width and one with a wider width, to accommodate models for different dynamical classes of KBOs (e.g., Brown 2001; Petit et al. 2011). vm17 adopted these models for the intrinsic distribution of free inclinations; we refer the reader to this paper for a detailed description. Van Laerhoven et al. (2019) obtained debiased estimates with slightly updated data and reported very similar results for the intrinsic inclination distribution models.

For the purposes of our goal to test and update the results of vm17 with the larger current observational sample of cataloged KBOs, we follow the choices made in that work for the functional form and parameters of the model distribution of free inclinations. Thus, as in vm17, for semimajor axis bins above 50 au we draw free inclinations from a Rayleigh distribution with scale parameter $\sigma =\sin 18^\circ ;$ below 50 au, the free inclinations are drawn from a mixture of two Rayleigh distributions with scale parameters ${\sigma }_{1}=\sin 3^\circ$ and ${\sigma }_{2}=\sin 13^\circ$. If the selected jth real object has an inclination below 5° and a random number drawn from the continuous uniform distribution on [0,1) is less than 0.9, we set σ = σ₁; otherwise σ = σ₂. vm17 used this rule to generate simulated objects approximately evenly split between the low-inclination and high-inclination Rayleigh distributions. The longitude of ascending node of the free inclination vector is drawn randomly from the continuous uniform distribution on [0, 2π). The total inclination vector of the simulated object is then computed as the sum of this randomly generated free inclination vector and the inclination vector of the mean plane of the N_K objects in the observational sample, as computed in Section 3; that is

$$\begin{eqnarray}&&({q}_{i},{p}_{i})=({q}_{0},{p}_{0})+({q}_{1,i},{p}_{1,i}),\end{eqnarray} \tag{ 9 }$$

where (q₀, p₀) is the inclination vector of the mean plane of the N_K observed objects and (q_1,i , p_1,i ) is the randomly generated free inclination vector. The remaining angular orbital elements—the mean anomaly and the argument of pericenter—are each drawn randomly and independently from the continuous uniform distribution on [0, 2π). With all the orbital elements of the simulated object in hand, we then compute its ecliptic latitude β_i and ecliptic longitude λ_i . Simulated objects are generated until there is exactly one simulated (β_i , λ_i ) pair for each (β, λ) pair in the observational sample, where, as in vm17, an acceptable match is within one degree in β and within five degrees in λ. Unfortunately, some objects (typically one or two per semimajor axis bin) are difficult to match in this manner. If an object has not been matched after 10⁷ random draws, we accept the following draw and move on to the next object. As explained in vm17 and Brown & Pan (2004), this procedure of generating the synthetic samples by matching their sky locations to those of the real observational sample approximately accounts for the selection biases present in the real observational sample, subject to the key assumption that each patch of sky in ecliptic coordinates (where a real object is found) has been thoroughly observed. It follows from this assumption that a true population will have the same number of objects in each patch of sky as the observed catalog, and we roughly control for the brightness and magnitude biases of the many surveys that comprise the entire KBO catalog by requiring the semimajor axis and eccentricity of each synthetic object to be close to that of a real object.

Having generated N_K simulated objects, we compute the mean plane of the simulated sample as described in Section 3. We repeat this procedure N_R times, producing a set of simulated mean planes (q_0,r , p_0,r ), where r is an integer from 1 to N_R . We use the R command dataellipse to compute the 68.2%, 95.4%, and 99.7% covariance ellipses of the set of simulated mean planes. We are now ready to produce uncertainty intervals for the inclination and longitude of ascending node of the mean plane of the observational sample, i₀ and Ω₀. If these were independent quantities defined on infinite domains, we would simply report low and high percentiles from the inclinations and longitudes of ascending node of the N_R synthetic mean planes, but both i₀ and Ω₀ have finite angular domains. Instead, we use the 68.2% (95.4%, 99.7%) covariance ellipse from the simulated mean planes (q_0,r , p_0,r ) and take the maximum and minimum values of i₀ and Ω₀ on the ellipse as the 68.2% (95.4%, 99.7%) confidence intervals for the inclination and longitude of node of the mean plane of the observational sample. If the covariance ellipse surrounds the origin in the (q, p) plane, the confidence interval for Ω₀ spans the entire circle, [0, 2π), and the confidence interval for i₀ has its lower endpoint at zero inclination. We also report a single number, σ_i , to describe the uncertainty of the pole position, in degrees, on the celestial sphere, computed as

$$\begin{eqnarray}&&{\sigma }_{i}=\arcsin \displaystyle \frac{1}{2}(a+b),\end{eqnarray} \tag{ 10 }$$

where a and b are, respectively, the semimajor and semiminor axes of the 68.2% covariance ellipse of the simulated mean planes.

For each semimajor axis bin, N_R = 40,000 synthetic mean planes are computed, which is sufficient to get convergence of the confidence intervals. Note that our grid search computes i₀ and Ω₀ to a precision of ∼01 and ∼1°, respectively, but our Monte Carlo procedure produces 68.2% confidence intervals of these parameters that are an order of magnitude larger; the latter is the 1σ uncertainty of the measured mean plane.

Laplace theory has been developed with the assumption that no two planets in the system are in mean motion resonances with each other, and no test particle is in a mean motion resonance with any planet. Before the mean plane of any observational population of KBOs can be computed for comparison with the Laplace surface, all resonant KBOs must be identified and removed.

To identify a suitable sample of KBOs with which we could calculate the mean plane, we first used the JPL Solar System Dynamics Group's Small Body Database query to retrieve all asteroid-type (as opposed to comet-type) objects with e < 1 and a heliocentric semimajor axis constrained to 30 < a < 200 au. On 2023 February 20, this returned heliocentric elements for 4149 objects. We eliminated all objects with a fractional semimajor axis uncertainty σ_a /a > 5%, as well as 11 noncometary objects where the semimajor axis uncertainty was unstated. Next, we downloaded the MPCORB.DAT database from the Minor Planet Center (MPC) on 2023 February 20 and cross-referenced it against the Small Body Database to eliminate all objects that have been observed for fewer than the three oppositions recommended by Gladman et al. (2008) and eliminated all objects with cometary designations. For each of the remaining objects, we retrieved barycentric elements on 2023 February 20 from JPL Horizons using the Python package Astroquery (Ginsburg et al. 2019). We then eliminated any remaining comets according to the criteria in Gladman et al. (2008), i.e., objects with a Tisserand parameter of T_J < 3.05 and a perihelion of q < 7.35 au. Next, we enforced barycentric semimajor axis limits of a < 150 au. We introduced a perihelion cutoff of q > 30.34 au (which is equal to the barycentric aphelion distance of Neptune on 2023 February 20) to eliminate any remaining planet-crossing objects. We then eliminated objects with barycentric a > 150 au. After the initial stages of sample selection, 2810 objects remained. We classified the remaining objects as resonant or nonresonant as described below, and eliminated the resonant objects from the sample.

The defining property of a resonant object is the libration of its critical resonance angle. The generally accepted method for identifying resonant KBOs in a sample population is to integrate the orbit of each KBO in an accurate n-body integrator, including perturbations from all the planets, recording all the resonant angles of interest, and then examine plots of the time series of the resonant angles by eye to look for persistent librations. The orbits must be integrated long enough to detect the libration of the longest-period resonant angles. Standard resonance identification methods, most notably that of Gladman et al. (2008), integrate for 10 Myr. For the longest-period angles, for angles that librate over large fractions of the circle, or for objects that alternate between libration and circulation, it may be difficult to distinguish resonant objects from nonresonant objects. Uncertainties in KBO orbital elements can also complicate resonant classification.

The current count of KBOs with well-determined orbits to classify (2810 objects) is rather large for efficient classification by eye. When assigning resonant classification by eye, it is important to maintain mental consistency so that objects at the beginning, middle, and end of the list receive the same scrutiny. Consistency of resonant evaluations across the sample may be improved by repeatedly shuffling the sample and reclassifying the objects, then accepting the most common classification for each object, or by deliberately evaluating each object according to a checklist of features that can be seen by eye. If a checklist is to be used, however, it is faster, more reliable, and more reproducible to automate the process by encoding those features in a classification algorithm (one example is Yu et al. 2018). The machine-classified objects may then be reexamined by eye to check for errors, if desired. Because new generations of telescopes are expected to increase dramatically the number of known KBOs in the upcoming decade, it is desirable to have an automated method to classify KBOs.

Smullen & Volk (2020) used the criteria of Gladman et al. (2008) to classify 2305 KBOs from the MPC as of 2016 October 20 after fitting new orbits to each object and integrating them for 10 Myr. They then used the Python machine-learning package Scikit-Learn (Pedregosa et al. 2011) to develop a gradient boosting classifier for the classification of KBOs as either resonant, classical, detached, or scattering, training the said classifier upon the orbits newly classified using the Gladman et al. (2008) criteria. Their code integrates a KBO orbit from initial barycentric elements in the n-body integrator rebound (Rein & Liu 2012) for 100 kyr and records 55 features of the orbit for use by the machine-learning algorithm. Full details of the settings used for the machine-learning algorithm, and a full explanation and justification of the 55 recorded features, are given in their paper. This machine-learning classifier was demonstrated to reproduce with 97% accuracy the classifications of their testing sample of 542 objects while using orbit integrations of only 1% the length of the 10 Myr standard. Smullen & Volk (2020) posted user-friendly python sample code and training data to GitHub to allow other researchers to use their gradient boosting classifier to classify KBOs by simply providing the barycentric orbital elements as inputs.

We downloaded the python sample code and training data (KBO_features.csv) from the Smullen & Volk (2020) GitHub repository. We used their gradient boosting classifier without modification, trained it on the same training set they provided, exactly as in the sample code, and used it to classify our sample of 2810 KBOs remaining after the down-selection procedure described above.

To account for orbital uncertainties, we used the JPL Small Body Database API (JPL Solar System Dynamics Group 2022) to download a JSON file for each object to be classified. The JSON file contained a nominal heliocentric orbital state, a 6 × 6 covariance matrix for the heliocentric orbit, and an epoch for the nominal orbit and the covariance matrix. The heliocentric orbital elements and their covariance were given as e, q, t_p , Ω, ω, and i, i.e., the eccentricity, perihelion distance in au, time of perihelion passage (Julian date), longitude of the ascending node, argument of perihelion, and inclination, where all angles are in degrees and referenced to the J2000 plane, and the epoch is a Julian date. JSON covariances for two of the 2810 objects were unavailable, so we could only classify 2808 objects.

We generated 301 heliocentric orbital element sets for each object: the nominal orbit and 300 clones from a Gaussian distribution centered at the nominal orbit, from the given covariance. The mean anomaly for each orbital element set was computed as the mean motion for the semimajor axis, times the elapsed time between the time of perihelion passage and the epoch. We used Horizons to download heliocentric orbital elements for the giant planets at each epoch. For each of the 2808 objects, we then built 301 rebound simulations consisting of the outer planets at the appropriate epoch and the orbital element set of the nominal orbit or the clone at the same epoch, treating the planets as massive particles and the KBO or the clone as a massless test particle.

Each rebound simulation was then classified as classical, scattering, detached, or resonant using the unmodified Smullen & Volk (2020) gradient boosting classifier. Of the 2808 cloned objects, 1812 had zero clones classified as resonant, 304 had 1–300 clones classified as resonant, and 692 had all 301 clones classified as resonant. To be absolutely sure no resonant objects contaminated our sample, we only accepted the 1812 objects for which no clones were resonant. Had we selected a different cutoff (50% resonant clones), our sample size would have been 150 objects larger. A further classification of the nonresonant objects as classical, detached, or scattering was not needed for our mean plane computations, so we did not further examine the classifications of the clones. Mean plane calculations then proceeded using the nominal orbits of the remaining objects. The complete set of 1812 nonresonant KBOs is provided online. A small sample is shown in Table 1. Nonresonant KBO counts are given by the semimajor axis bin in Table 2. This table also contains KBO counts by bin for vm17, for comparison. The online supplementary material for this paper includes a ZIP archive with the JSON files that contain the JPL covariance matrix for each object.

Table 2. Measured Kuiper Belt Mean Plane

Semimajor Axis Bin (au)	Count	Count (vm17)	i0 (deg)	Ω0 (deg)	σi (deg)
35–40.3	92	43	${3.0}_{-2.3}^{+1.8}$	${121}_{-37}^{+25}$	1.59
40.3–42	164	82	${4.1}_{-0.9}^{+2.3}$	${272}_{-12}^{+17}$	1.40
42–43	207	100	${1.6}_{-0.7}^{+0.8}$	${64}_{-46}^{+24}$	0.81
43–44	360	186	${2.1}_{-0.6}^{+0.6}$	${96}_{-15}^{+11}$	0.54
44–45	273	141	${1.8}_{-0.7}^{+0.6}$	${85}_{-22}^{+15}$	0.59
45–48	402	194	${1.3}_{-0.4}^{+0.7}$	${3}_{-31}^{+25}$	0.60
45–50	440	217	${1.4}_{-0.4}^{+0.8}$	${355}_{-26}^{+22}$	0.59
50–80	221	125	${2.0}_{-2.0}^{+3.2}$	${114}_{-114}^{+246}$	4.02
50–150	276	162	${4.1}_{-4.1}^{+3.1}$	${138}_{-138}^{+222}$	3.66
Total (35–150)	1812	931	${1.6}_{-0.4}^{+0.2}$	${82}_{-15}^{+6}$	0.28
Classical (42–48)	1242	617	${1.6}_{-0.4}^{+0.3}$	${82}_{-14}^{+8}$	0.31
Total – classical	570	314	${1.7}_{-1.0}^{+0.4}$	${149}_{-15}^{+39}$	0.65

Note. The best fit and 68.2% confidence intervals of the measured mean plane's ecliptic inclination, i₀, longitude of ascending node, Ω₀, and the uncertainty of the pole position, σ_i .

Download table as:  ASCIITypeset image

Note that Table 1 is given for illustration only. Neither the low precision of the orbital elements in the printed table nor the high precision of the elements in the electronic table represents the true statistical precision of the JPL Horizons orbits. This is because the uncertainties reported in the JSON covariance matrix are not straightforwardly related to the uncertainties of the orbital elements reported in Table 1. The JPL covariance matrix is given in a different set of orbital elements that must be transformed to the standard set (q to a, t_p to M). Table 1 also contains elements (r, λ, β) that are not linearly related to the standard set. Also, the JPL covariance matrix is given at a different epoch for each object, and Table 1 reports the elements of each object at a common epoch. It is well known that Gaussian uncertainty regions around a nominal orbit do not stay Gaussian as the nominal orbit is integrated forward or backwards in time, although the deviation from Gaussianity is small for small elapsed times.

With the nonresonant objects identified as in Section 5, and following the mean plane calculation method described in Section 3, and its confidence interval calculation method described in Section 4, the mean planes and confidence intervals we found are shown in graphical form in Figures 2–4, and tabulated in Table 2.

Zoom In Zoom Out Reset image size

In Figure 2, the left panel shows the results for the classical Kuiper Belt (42–48 au, 1242 objects), the center panel shows the results for the entire Kuiper Belt (35–150 au, 1812 objects), and the right panel shows the results for the entire Kuiper Belt with the classical region excluded (570 objects). The best-fit mean plane is shown as a dark green +, and the 68.2%, 95.4%, and 99.7% confidence ellipses are also indicated in dark green. For comparison, the semimajor-axis-dependent Laplace surface (from linear secular theory as developed in Murray & Dermott 1999, recomputed with updated planetary parameters for the epoch 2023 February 20) is indicated in blue. For context, we also indicate the location of the invariable plane with a black x. The J2000 ecliptic/equinox pole is located at the origin. We observe that the classical region dominates the results: when it is removed, i₀ is nearly unchanged, but Ω₀ shifts by +67°, and the dispersion σ_i of the synthetic mean planes more than doubles.

In Figure 3, we plot the results for the mean planes computed at higher resolution in semimajor axis: each panel shows the results for a single semimajor axis bin from Table 2. As in Figure 2, we plot the measured mean plane and its uncertainty distribution, and indicate the invariable plane and the Laplace surface. In these panels, we also indicate the best-fit mean plane from vm17 as a magenta +, with the 68.2% confidence ellipses from vm17 in magenta and with estimated 95.4% and 99.7% confidence ellipses scaled from the 68.2% ellipses. Note that the Laplace surface appears as a very concentrated set of dots in most panels, but in panels (a) and (b) it is an extended quasi-linear set of dots owing to the warps caused by the ν₁₇ and the ν₁₈ secular resonances, respectively. As a increases from the lower boundary of each semimajor axis bin to its upper boundary, the Laplace surface traces these paths: in panel (a), the trace begins from the left, at high i and Ω ≈ 120°, passes near the origin with i ≈ 18, and exits to the lower right; in panel (b), the trace begins at the upper right and approaches the invariable plane.

Zoom In Zoom Out Reset image size

Our results for the mean plane inclination and longitude of node, and their 68.2% confidence intervals for each semimajor axis bin are tabulated in Table 2. These are also plotted in Figure 4 as a function of semimajor axis. The best-fit mean plane results are in green, where the horizontal error bars indicate the width of the semimajor axis bin and the vertical error bars indicate the 68.2% confidence interval. For comparison, we plot the Laplace surface (varying with semimajor axis) in blue and the invariable plane is indicated by the black horizontal line. We also indicate the results of vm17 in magenta (with a small horizontal offset, for legibility).

Zoom In Zoom Out Reset image size

Is the measured mean plane consistent with the Laplace surface or with the invariable plane? For the classical Kuiper Belt from 42 to 48 au, we find the mean plane of 1242 objects as ${i}_{0}=(1\buildrel{\circ}\over{.} 6{)}_{-0.4}^{+0.3}$, ${{\rm{\Omega }}}_{0}=(82^\circ {)}_{-14}^{+8}$. For the entire belt from 35 to 150 au, we find the mean plane of 1812 objects as ${i}_{0}=(1\buildrel{\circ}\over{.} 6{)}_{-0.4}^{+0.2}$, ${{\rm{\Omega }}}_{0}=(82^\circ {)}_{-14}^{+8}$. When we exclude the classical region, we find the mean plane of 570 objects as ${i}_{0}=(1\buildrel{\circ}\over{.} 7{)}_{-1.0}^{+0.4}$, ${{\rm{\Omega }}}_{0}=(149^\circ {)}_{-15}^{+39}$. In all three cases, this is close to the invariable plane (respectively within 07, 06, and 12), but distinguishable from it at greater than 99.7% confidence, as evident in Figure 2.

We do not comment on any distinction between the measured mean plane and the Laplace surface for these broad ranges, as the changing location of the Laplace surface renders the question ill formed, but we can answer this question for the smaller semimajor axis bins. In Table 3, we tabulate simple acceptance/rejection decisions with respect to the invariable plane and the Laplace surface for each semimajor axis bin. At the 99.7% confidence level, we can only reject the invariable plane or the Laplace surface for the bins 40.3–42, 45–48, and 45–50 au. In each of those bins, we reject both planes. Elsewhere, both planes are consistent with the sample population.

The discrepancy between the theoretically expected Laplace surface and the measured mean plane in the 40.3–42 au bin may arise from inaccuracies of the linear theory for the theoretical estimate of the strong warp in proximity to the ν₈ secular resonance. In the case of the 45–48 au and 45–50 au bins, we note that both these bins are in proximity to Neptune's 2:1 mean motion resonance but the theoretical Laplace surface does not account for the effects of mean motion resonances.

Because our methodology closely follows that of vm17 and our main motivation was to test the reproducibility of their results with the updated sample of KBOs, we comment in some detail on the comparisons between those results and ours. In general, the results reported here have smaller measurement uncertainties of the mean pole positions than those in vm17, due to the larger sample sizes now available.

• 1.  
  In the 35–40.3 au bin (Figure 3(a)), almost our entire 99.7% confidence ellipse falls within their 68.2% confidence ellipse. While our 68.2% confidence ellipse is small enough to identify Ω₀ uniquely at that level of confidence, our 99.7% ellipse surrounds the origin so we cannot report a unique Ω₀ at the higher confidence level.
• 2.  
  In the 40.3–42 au bin (Figure 3(b)), our 68.2% confidence ellipse falls mostly outside theirs and is closer to the origin, resulting in a similar 68.2% confidence interval for Ω₀ but a distinctly different (and lower) estimate and interval for i₀. Our results are closer to the Laplace surface than theirs, but still differ from the Laplace and invariable planes by more than 99.7% confidence. With future increases in the size of KBO samples, splitting this bin into two or more narrower semimajor axis ranges could enable higher precision comparisons of the measured and predicted Laplace surface within the constraints of the existing theory. Without splitting the bin, the ν₁₈ secular resonance warps the Laplace surface so strongly with semimajor axis that a second-order theory is necessary for accurate predictions.
• 3.  
  In the 42–43, 43–44, and 44–45 au bins (Figures 3(c), (d), and (e), respectively), our measured mean planes fall within their 68.2% covariance ellipses, and we produce overlapping, though smaller, 68.2% confidence intervals.
• 4.  
  In the 45–48 and 45–50 au bins (Figures 3(f) and (g), respectively), our measured mean planes again fall within their 68.2% covariance ellipses, but our 95.4% and 99.7% covariance ellipses are small enough that, unlike vm17, we reject the Laplace surface and invariable plane at all three confidence levels. These bins are close to Neptune's 2:1 outer mean motion resonance.
• 5.  
  Most strikingly, we do not detect the strong warp vm17 reported for the distant Kuiper Belt near a > 50 au (Figures 3(h) and (i)). In the 50–80 and 50–150 au bins, the Laplace surface and invariable plane are both within 29 of our measured mean planes, and fall within our 68.2% covariance ellipses, while for vm17 they are more than 6° away from the measured mean planes and fall near the 99.7% covariance ellipses.We cannot currently rule out the invariable plane, the Laplace surface, or even the ecliptic plane as the mean plane of the Kuiper Belt in this region at even the 68.2% confidence level. We note that our sample size in this region is 272, which is 70% larger than vm17's sample. Although the best estimate of the mean plane of the distant Kuiper Belt differs by more than 2σ in the two studies, there is significant overlap of the 3σ confidence ellipses of our and their measurements. Larger sample sizes in the future would help to understand this region better by reducing the measurement uncertainties. If we assume, as a rule of thumb, that the size of the uncertainty ellipse varies with the sample size as ∼N^−1/2, then reducing the measurement uncertainty by a factor of 2–3 will require a sample size 4–9 times as large. This is likely to be achieved over the next decade as the Vera C. Rubin Observatory carries out the Legacy Survey of Space and Time (LSST; Ivezić et al. 2019).

We briefly comment on comparisons of our results with those of other previous studies.

• 1.  
  Chiang & Choi (2008) computed the Kuiper Belt mean plane and its uncertainty for 10 objects with 38.09 < a < 39.10 au, e < 0.1, and ecliptic inclination i < 10°, and for 80 objects with 42.49 < a < 43.50 au and the same eccentricity and inclination cutoffs. They did not rule out either the Laplace surface or the invariable plane as the true Kuiper Belt plane for those two specific semimajor axis ranges with greater than 3σ (99.7%) confidence.When we impose the same eccentricity and inclination restrictions on our sample, we respectively obtain 47 and 347 objects in these semimajor axis ranges. When we compute the mean plane and its uncertainty for these samples, we do not rule out the Laplace surface at any confidence level for either sample, and we do not rule out the invariable plane for either sample at a confidence level above 68.2%. Removing the eccentricity and semimajor axis cutoffs adds no objects to either semimajor axis range, and does not change the outcome.
• 2.  
  Van Laerhoven et al. (2019) computed the plane of the cold classical Kuiper Belt, with uncertainties, for 107 objects in the 42.4–43.8 au bin, 82 objects in the 43.8–44.4 au bin, and 67 objects in the 44.4–47 au bin, where a cold object was defined as having an inclination of 4° or less relative to the Laplace surface at the center of the semimajor axis bin (i.e., 43.1 au for the 42.4–43.8 au bin). They computed the Kuiper Belt plane and its uncertainties for an unspecified number of hot classical objects between 42.4 and 47 au, where a hot object has an inclination of 9° or more relative to the Laplace surface at 44.7 au. They repeated the computations, without inclination restrictions, for 57 objects in the 50–80 au bin and for 83 objects in the 48–150 au bin. In the 42.4–43.8 au and 43.8–44.4 au bins of cold objects, they rejected the invariable plane but not the Laplace surface at 99% confidence. In the 44.4–47 au bin of cold objects, the 42.4–47 au bin of hot objects and the higher-a bins (a > 47 au), they rejected neither the invariable plane nor the Laplace surface at 99% confidence.When we impose the same semimajor axis and inclination restrictions on our sample, we find 450 cold objects in the 42.4–43.8 au range, 284 cold objects in the 43.8–44.4 au range, and 444 cold objects in the 44.4–47 au range. We find 1178 hot objects in the 42.4–47 au range, and 314 total objects in the 48–150 au range. In each case, our sample sizes are larger than Van Laerhoven et al. (2019)'s by factors of three or more. In the 42.4–43.8 au bin of cold objects, we reject the invariable plane at 95.4% confidence but do not reject the Laplace surface at any confidence level. In the 43.8–44.4 au bin of cold objects, we reject the invariable plane at 68.2% confidence but do not reject the Laplace surface at any confidence level. These results are similar to those shown for nearby semimajor axis bins in Figure 3. In the 44.4–47 au bin of cold objects, we reject neither the invariable plane nor the Laplace surface at any confidence level. In the 42.4–47 au bin of hot objects, we reject both the invariable plane and the Laplace surface (for 44.7 au, the center of the bin) at well above 99.7% confidence. As previously noted, for the higher-a bins we reject neither the invariable plane nor the Laplace surface at 95.4% confidence.Our results for the cold population agree with those of Van Laerhoven et al. (2019), but we reject at high confidence both the invariable plane and the Laplace surface for the hot objects in the 42.4–47 au range where they rejected neither. Our simulated mean planes in that bin are tightly clustered around our observed mean plane, while theirs were tightly clustered around the invariable plane. We suspect that this discrepancy is explainable by our much larger sample size; perhaps the mean plane of their much smaller sample was located much closer to the invariable plane.
• 3.  
  Brown & Pan (2004) rejected the invariable plane at the 3σ (99.7%) confidence level for the Kuiper Belt as a whole. We also reject it for the entire belt between 35–150 au.
• 4.  
  Elliot et al. (2005) rejected the Laplace surface at the 1σ (68.2%) confidence level for 85 "classical" objects in the range 37.9 < a < 47.0 au. In their definition, a classical object is a nonresonant object with e ≤ 0.2 and a Tisserand parameter relative to Neptune of T_N > 3.When we impose the same eccentricity and Tisserand restrictions on this semimajor axis range, we obtain 1507 objects. When we compute the mean plane and its uncertainty for these 1451 objects, we reject the invariable plane at 99.7% confidence but do not reject the Laplace surface at any confidence level. Removing the eccentricity and Tisserand cutoffs adds no objects to the semimajor axis range, and does not change the outcome.

We close with the observation that future increases in the sample sizes of outer solar system minor planets will enable higher confidence results in measurements of the Kuiper Belt mean plane, and should be examined for their potential to detect the effects of unseen distant planets or other unmodeled perturbations.

We thank Aaron Rosengren, James Keane, and especially Kathryn Volk for correspondence, comments, and sample code. We acknowledge funding from the NSF (grant AST-1824869) and from the JPL SURP and SIP student support programs. We also thank the anonymous reviewer and the data editor and statistics editor for their comments, questions, and advice.

Software: astroquery (Ginsburg et al. 2019), rebound (Rein & Liu 2012), scikit-learn (Pedregosa et al. 2011), shapely (Gillies et al. 2022), KBO classifier (Smullen & Volk 2020).