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

Mean plane measurements of the Kuiper Belt from observational data are of interest for their potential to test dynamical models of the solar system. Recent measurements have yielded inconsistent results. Here we report a measurement of the Kuiper Belt’s mean plane with a sample size more than twice as large as in previous measurements. The sample of interest is the nonresonant Kuiper Belt objects, which we identify by using machine learning on the observed Kuiper Belt population whose orbits are well determined. We estimate the measurement error with a Monte Carlo procedure. We find that the overall mean plane of the nonresonant Kuiper Belt (semimajor axis range of 35–150 au) and also that of the classical Kuiper Belt (semimajor axis range of 42–48 au) are both close to (within ∼0.°7) but distinguishable from the invariable plane of the solar system to greater than 99.7% confidence. When binning the sample into smaller semimajor axis bins, we find the measured mean plane is mostly consistent with both the invariable plane and the theoretically expected Laplace surface forced by the known planets. Statistically significant discrepancies are found only in the semimajor axis ranges 40.3–42 au and 45–50 au; these ranges are in proximity to the ν 8 secular resonance and Neptune’s 2:1 mean motion resonance where the theory for the Laplace surface is likely to be inaccurate. These results do not support a previously reported anomalous warp at semimajor axes above 50 au.

1. INTRODUCTION Chiang & Choi (2008) posed the question:"If we could map, at fixed time, the instantaneous locations in threedimensional 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 4 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 1.58 • and a longitude of ascending node of 107.58 • 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 0 , Ω 0 ) = (1.86 • , 81.6 • ), where i 0 and Ω 0 are respectively the inclination and longitude of ascending node of the plane in J2000 coordinates.With a measurement error of 0.37 • 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 0 in the range 1.65 • − 2.49 • and Ω 0 in the range 97.4 • − 104.0 • , with a preferred value of (i 0 , Ω 0 ) = (1.51• ± 0.26 • , 100.0 • ± 8.8 • ).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 Kuiper Belt objects 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 0 , Ω 0 ) = (1.8• , 77 • ) for the classical Kuiper Belt, i.e., the sample of non-resonant 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 au to 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 non-resonant KBOs with well-determined orbits, semimajor axes between 35-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 for identifying the sample of non-resonant KBOs.Our results are presented and discussed in Section 6.

LAPLACE SURFACE
Though their measurements of the mean plane differ, researchers including vm17, Chiang &Choi (2008), andBrown &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 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 quasi-periodically 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, ĥ,  normal to the orbit plane can be expressed in terms of these elements (Eq.5).The Laplace theory expositions in the literature use the more convenient two-component vector, the so-called "inclination vector" defined as (q, p) = sin i (cos Ω, sin Ω). (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, (q, p) = (q 0 , p 0 ) + (q 1 , p 1 ). (2) The first part, (q 0 , p 0 ) = sin i 0 (cos Ω 0 , sin Ω 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 Ω 1 , sin Ω 1 ), called the "free inclination vector," is the remaining part of the total inclination vector; it is determined by initial conditions.The Laplace surface changes over secular times, 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 1 remains constant and the free longitude of node Ω 1 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, n0 , normal to the local Laplace surface is related to the forced inclination vector as follows, n0 = p 0 , −q 0 , 1 − p 2 0 − q 2 0 . (3) We briefly make a note about heliocentric versus barycentric orbital elements.The Laplace surface on 2023 February 20 for semi-major 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 eccentricity, inclination, longitude of node, longitude of pericenter, and semimajor axis.Only 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 au and 40.3 au.

MEAN PLANE MEASUREMENT
At first sight, it would seem reasonable to simply 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;Volk & Malhotra 2017).Observational surveys not only have brightness limits, but surveys of the outer solar system for KBOs also have limited and non-uniform 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, vt , 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 ĥ the unit vector normal to the orbit plane, and with r the unit vector along the radial direction to the KBO's heliocentric position, then vt = ĥ × r. (4) The unit vector, ĥ is found from a KBO's inclination and longitude of ascending node, i and Ω, and the unit vector r is found from its J2000 ecliptic latitude and ecliptic longitude, β and λ, follows, ĥ = sin i sin Ω, − sin i cos Ω, cos i , (5) r = cos β cos λ, cos β sin λ, sin β .(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 n defining the normal to the mean plane, so that n • vt = 0. (7) The mean plane can then be computed as the plane whose orbit normal is most normal to the most KBO sky-plane velocities, i.e., the plane that minimizes The function S may be minimized on 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 x 700 grid with q min = p min = −0.35,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 0.1 degrees.This is small compared to the confidence intervals that are computed later.The mean plane (q 0 , p 0 ) is trivially converted to an inclination and longitude of node (i 0 , Ω 0 ).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., Fig. 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.

MEAN PLANE MEASUREMENT UNCERTAINTY
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 non-resonant 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. (2011Petit et al. ( , 2017))) and the Outer Solar Systems Origins Survey (OSSOS, Bannister et al. (2018)).The data from these surveys is more effectively debiased with the use of survey simulators.It is 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 catalogued 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 σ = sin 18 • ; below 50 au, the free inclinations are drawn from a mixture of two Rayleigh distributions with scale parameters σ 1 = sin 3 • , σ 2 = sin 13 • .If the selected jth real object has inclination below 5 • and a random number drawn from the continuous uniform distribution on [0,1) is less than 0.9, we set σ = σ 1 , otherwise σ = σ 2 .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, where (q 0 , p 0 ) 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 7 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 0 and Ω 0 .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 0 and Ω 0 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 0 and Ω 0 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 Ω 0 spans the entire circle, [0, 2π), and the confidence interval for i 0 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 where a and b are 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 0 and Ω 0 to a precision of ∼ 0.1 • 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.

SAMPLE SELECTION
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 the suitable sample of KBOs for which to 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 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 fractional semimajor axis uncertainty σ a /a > 5%, as well as the eleven non-cometary objects where the semimajor axis uncertainty was unstated.Next, we downloaded the MPCORB.DAT database from the 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 at 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 non-Resonant 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 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 orbit 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 non-resonant objects.Uncertainty 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 re-classifying 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 re-examined by eye to check for errors, if desired.Because new generations of telescopes are expected to dramatically increase 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 Minor Planet Center (MPC) database 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 classification of KBOs as either Resonant, Classical, Detached, or Scattering, training 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 6x6 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 , Ω, ω, i, i.e. 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 centred 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 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 Kuiper belt object 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 non-Resonant 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 non-Resonant KBOs is provided online.A small sample is shown in Table 1.Non-Resonant KBO counts are given by 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.
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.

RESULTS AND DISCUSSION
With the non-Resonant 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.
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 0 is nearly unchanged, but Ω 0 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 ν 17 and the ν 18 secular resonance, respectively.As a increases from the lower boundary of each semi-major 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 ≈ 1.8 • , and exits to the lower right; in panel (b), the trace begins at the upper right and approaches the invariable plane.
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).
Is the measured mean plane consistent with the Laplace surface or with the invariable plane?For the classical Kuiper Belt from 42-48 au, we find the mean plane of 1242 objects as For the entire belt from 35-150 au, we find the mean plane of 1812 objects as When we exclude the classical region, we find the mean plane of 570 objects as In all three cases, this is close to the invariable plane (respectively within 0.7 • , 0.6 • , and 1.2 • ), 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 ν 8 secular resonance.In the case of the 45-48 au and 45-50 au bins, we note that both these bins are in proximity The best-fit mean plane is the dark green +.The mean planes of 40,000 Monte-Carlo samples are in light green, and the 68.2%, 95.4%, and 99.7% covariance ellipses for them are in dark green.The J2000 ecliptic/equinox pole is at the origin.The invariable plane is the black x.The theoretical prediction (from linear secular theory) for the Laplace surface as a function of semimajor axis is plotted in blue.More details are given in the main text.

0.65
Note-The best-fit and 68.2% confidence intervals of the measured mean plane's ecliptic inclination, i0, and longitude of ascending node, Ω0, and the uncertainty of the pole position, σi.
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.
• In the 35-40.3au bin (Figure 3a), almost our entire 99.7% confidence ellipse falls within their 68.2%confidence ellipse.While our 68.2%confidence ellipse is small enough to uniquely identify Ω 0 at that level of confidence, our 99.7%ellipse surrounds the origin so we cannot report a unique Ω 0 at the higher confidence level.• In the 40.3-42 au bin (Figure 3b), our 68.2%confidence ellipse falls mostly outside theirs and is closer to the origin, resulting in a similar 68.2% confidence interval for Ω 0 but a distinctly different (and lower) estimate and interval for i 0 .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 ν 18 secular resonance warps the Laplace surface so strongly with semimajor axis that a second-order theory is necessary for accurate predictions.
• In the 45-48 and 45-50 au bins (Figure 3f, g), 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.
• Most strikingly, we do not detect the strong warp vm17 reported for the distant Kuiper Belt near a > 50 au (Figure 3h, i).In the 50-80 and 50-150 au bins, the Laplace surface and invariable plane are both within 2.9 • 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 better understand this region 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.
• 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.43.8 au and 43.8-44.4au 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.8au range, 284 cold objects in the 43.8-44.4au 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.8au 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.4au 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.
• 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.
• 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 non-Resonant 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  The mean planes and confidence intervals from this work are in green, and those from vm17 are in magenta.For the sake of readability, the vertical magenta and green lines have been slightly offset from each other in semimajor axis when they would otherwise overlap.The theoretical prediction for the Laplace surface as a function of semimajor axis is plotted in blue.The invariable plane is indicated by the horizontal black line.

Figure 2 .
Figure 2. Kuiper Belt mean plane and its confidence ellipses in the (q, p) plane, by semimajor axis bin.Semimajor axis bins are (a) the classical belt from 42-48 au, (b) the entire belt from 35-150 au, and (c) the entire belt minus the classical region.The best-fit mean plane is the dark green +.The mean planes of 40,000 Monte-Carlo samples are in light green, and the 68.2%, 95.4%, and 99.7% covariance ellipses for them are in dark green.The J2000 ecliptic/equinox pole is at the origin.The invariable plane is the black x.The theoretical prediction (from linear secular theory) for the Laplace surface as a function of semimajor axis is plotted in blue.More details are given in the main text.

Figure 3 .
Figure 3. Kuiper Belt mean plane and its confidence ellipses in the (q, p) plane, by semimajor axis bin.Semimajor axis bins are (a) 35-40.3au, (b) 40.3-42 au, (c) 42-43 au, (d) 43-44 au, (e) 44-45 au, (f) 45-48 au, (g) 45-50 au, (h) 50-80 au, and (i) 50-150 au.The best-fit mean plane is the dark green +.The mean planes of 40,000 Monte-Carlo samples are in light green, and the 68.2%, 95.4%, and 99.7% covariance ellipses for them are in dark green.The J2000 ecliptic/equinox pole is at the origin.The invariable plane is the black x.The best-fit mean plane from vm17 is the magenta +.The theoretical prediction (from linear secular theory) for the Laplace surface as a function of semimajor axis is plotted in blue.More details are given in the main text.The 68.2%, 95.4%, and 99.7% covariance ellipses from vm17 are plotted in magenta.Each subplot has the same scale, showing that increased sample size does reduce mean plane uncertainty: in each semimajor axis bin, we have more objects and smaller covariance ellipses than vm17, and the more heavily populated semimajor axis bins on the second row have smaller covariance ellipses than the less heavily populated bins on the first row.The larger covariance ellipses on the third row, above 50 au, reflect the wider model inclination distribution used for those bins.

Figure 4 .
Figure4.Kuiper Belt mean plane and 68.2% confidence intervals by semimajor axis bin.The mean planes and confidence intervals from this work are in green, and those from vm17 are in magenta.For the sake of readability, the vertical magenta and green lines have been slightly offset from each other in semimajor axis when they would otherwise overlap.The theoretical prediction for the Laplace surface as a function of semimajor axis is plotted in blue.The invariable plane is indicated by the horizontal black line.

Table 1 .
Sample Table of Non-Resonant KBOs., and r are provided by JPL Horizons, and β and λ are calculated from Eq. 6.The complete table is available online.

Table 3 .
Statistical significance of the measured mean plane compared to the Invariable Plane (IP) and Laplace surface (LS) Note-Acceptance/rejection decisions for invariable plane and Laplace surface compared to the measured mean plane in each semimajor axis bin.Checks are accepted, ellipses are rejected.A check indicates that the 68.2% (95.4%, 99.7%) confidence ellipse of the measured mean plane pole for the semimajor axis bin contains the pole of the invariable plane, or it contains the pole of the Laplace surface for any semimajor axis within the boundaries of the bin.An ellipsis indicates the contrary.

•
Van Laerhoven et al. (2019)computed the plane of the cold classical Kuiper Belt, with uncertainties, for 107 objects in the 42.4-43.8au bin, 82 objects in the 43.8-44.4au 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.8au bin).They computed the Kuiper belt plane and its uncertainties for an unspecified number of hot classical objects between 42.4-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-