Beyond Point Masses. II. Non-Keplerian Shape Effects Are Detectable in Several TNO Binaries

About 40 trans-Neptunian binaries (TNBs) have fully determined orbits with about 10 others being solved except for breaking the mirror ambiguity. Despite decades of study, almost all TNBs have only ever been analyzed with a model that assumes perfect Keplerian motion (e.g., two point masses). In reality, all TNB systems are non-Keplerian due to nonspherical shapes, possible presence of undetected system components, and/or solar perturbations. In this work, we focus on identifying candidates for detectable non-Keplerian motion based on sample of 45 well-characterized binaries. We use MultiMoon, a non-Keplerian Bayesian inference tool, to analyze published relative astrometry allowing for nonspherical shapes of each TNB system’s primary. We first reproduce the results of previous Keplerian fitting efforts with MultiMoon, which serves as a comparison for the non-Keplerian fits and confirms that these fits are not biased by the assumption of a Keplerian orbit. We unambiguously detect non-Keplerian motion in eight TNB systems across a range of primary radii, mutual orbit separations, and system masses. As a proof of concept for non-Keplerian fitting, we perform detailed fits for (66652) Borasisi-Pabu, possibly revealing a J 2 ≈ 0.44, implying Borasisi (and/or Pabu) may be a contact binary or an unresolved compact binary. However, full confirmation of this result will require new observations. This work begins the next generation of TNB analyses that go beyond the point mass assumption to provide unique and valuable information on the physical properties of TNBs with implications for their formation and evolution.


INTRODUCTION
Since the discovery and characterization of the mutual binary orbit of the transneptunian object (TNO) 1998 WW31 (Veillet et al. 2002), transneptunian bina-Corresponding author: Benjamin Proudfoot benp175@gmail.comries (TNBs) have been recognized as sensitive tracers of the history of the solar system.Acting as a detailed laboratory that enables mass measurements, TNBs open the door for remote characterization of TNOs as a whole (e.g., Grundy et al. 2007;Grundy et al. 2009;Fraser & Brown 2010; Barr & Schwamb 2016).In addition, the mutual orbital properties of a binary system provide insight into the formation and history of that binary system, as those properties encode information about the binary's formation (e.g., Brown & Schaller 2007;Brown et al. 2010), past tidal evolution (e.g., Porter & Grundy 2012;Arakawa et al. 2021), collisional history (e.g., Ragozzine & Brown 2009;Parker & Kavelaars 2011), and encounters with other bodies (e.g., Campbell et al. 2022).The statistical ensemble of mutual orbit properties of TNBs also hold valuable clues about the conditions of the protoplanetary disk from which TNBs originally formed and have revealed groundbreaking insights into the dominant formation processes in that disk (e.g., Nesvornỳ et al. 2010;Grundy et al. 2019;Nesvornỳ et al. 2019).
While the numerous studies of TNBs have enabled unprecedented understanding of processes in the outer solar system, it has become clear that current state-ofthe-art methods are lagging behind the growing observational baselines of TNBs.These methods, for the most part, rely on simple Keplerian orbital models, with only a few exceptions.In recent years, however, many authors have found that the observed relative astrometric positions of TNBs show statistically significant deviations from plain Keplerian orbits (Salacia-Actaea 3.7σ, Orcus-Vanth 2.2σ, Grundy et al. (2019); Eris-Dysnomia 6.3σ, Holler et al. (2021)).An analysis of the distribution of χ 2 values of all TNB fits also shows a statistically significant excess of poor fits.While these deviations could be the result of unidentified systematic errors in data collections/measurement, it is also likely that the deviations are the manifestation of inaccurate models including non-Keplerian gravitational effects acting in TNB systems.
Non-Keplerian gravitational effects are any gravitational effect that forces an orbit to deviate from a pure Keplerian orbit.Generally, non-Keplerian effects result in slow precession of an orbit's orientation angles.Precession of the direction of periapse is called apsidal precession, and precession of the orbit pole is called nodal precession.While there are many possible sources of non-Keplerian effects, the most relevant sources for TNBs are the non-spherical shapes of individual TNB components, systems with more than 2 components (whether known or unknown), and the external gravitational influence of the Sun (see Ragozzine et al. (2023) for more details).
In this paper, we will focus on detecting and measuring the strength of non-Keplerian effects attributable to shape or unknown components, leaving non-Keplerian effects from solar influences to future study.In Section 1.1, we discuss the causes and consequences of non-Keplerian shape effects.Next, in Section 1.2, we explain the general process of our non-Keplerian orbit fitting.In Section 2.2 and Section 2.3, we discuss our methods for both Keplerian and non-Keplerian orbital fitting applied herein.In Section 2.1, we detail our TNB sample and gather publicly available relative astrometric data.In Section 3, we discuss the results of our Keplerian orbit fits and reproduce past fitting results.Next, in Section 4, we present the results of our search for non-Keplerian effects, discuss the implications of these results, and identify the most promising targets for future investigation.Then, in Section 5, we perform a full non-Keplerian analysis of TNB (66652) Borasisi-Pabu, one of the most promising cases of non-Keplerian motion in a TNB, as a proof-of-concept of full non-Keplerian fitting.Lastly, in Section 6, we discuss our conclusions.

Non-Keplerian Shape Effects
Keplerian orbital models implicitly assume that the individual components of a TNB system are point masses (or equivalently perfect spheres).However, it is well-documented that the shapes of TNOs can be significantly non-spherical (e.g., Sheppard & Jewitt 2004;Ortiz et al. 2017).Moving beyond point masses, nonspherical shapes must cause non-Keplerian deviations in TNB orbits, though the importance of these deviations relative to present or future observational data has not yet been previously examined in detail.
Gravitational potentials of non-spherical bodies can be modeled using a spherical harmonic expansion of the gravitational potential.Current data warrant exploring the gravitational potential only at second order; by construction, higher order corrections are less important and neglected in our analysis.Ragozzine et al. (2023) provide a detailed discussion of these and other effects and how they can be modeled which we summarize here.
The second-order gravitational potential, U , of a mass M at distance r, can be written: (1) where J 2 is the second-order zonal gravitational harmonic and C 22 is the second-order sectoral gravitational harmonic coefficient, θ is the body-fixed latitude-like angle, ϕ is the body-fixed longitude-like angle (chosen to eliminate other terms), and R is a reference radius (Yoder 1995;Scheeres et al. 2000).J 2 is a measure of the oblateness of the potential and C 22 is related to the prolateness, or the ellipticity of the equator.Generally, for extremely spherical bodies (like the Earth) J 2 ≲ 0.001, extremely oblate bodies (like Haumea) J 2 ∼ 0.1, and contact binaries have J 2 ∼ 0.3.
For most TNBs, where the TNBs mutual orbit is much slower than the rotations of the individual components, C 22 has little effect on the dynamics of the mutual orbit, as the contribution from C 22 averages out.However, near spin-orbit resonances, which may be present among some TNB systems, C 22 may not fully average out and may play a significant role in the dynamics of the system (Proudfoot & Ragozzine 2021).We consider this to be a special case not relevant to most TNBs.Since we are exploring the entire ensemble of TNBs, we focus this analyses on J 2 alone.Implications for this choice are discussed further in Section 5.
The dynamics of a system with a J 2 can be described by slow apsidal and nodal precession.For a test particle orbiting around a body with a J 2 , the apsidal and nodal secular rates can be written: where n is the mean motion, a is the semi-major axis, e is the eccentricity, and i is the inclination of the orbit relative to the body's equator.As discussed in Ragozzine & Brown (2009), and is clear from Equations 2 and 3, the non-Keplerian motion induced by the gravitational harmonics require knowledge of both the strength of the gravitational harmonics and the direction of the spin axis (to determine the appropriate θ and ϕ values in Equation 1 or to determine the appropriate inclination in equations 2 and 3).Thus, detection of non-Keplerian effects allows constraints to be placed on both TNO shapes and spin poles (assumed to be identical to the quadrupolar gravitational harmonic pole appropriate for TNOs which are too large to sustain long-term non-principal axis rotation).
From Equations 1, 2, and 3, it can be seen that the measurable quantity is J 2 R 2 (and C 22 R 2 ) and not J 2 independently.This is extremely similar to the quantity GM , where the gravitational constant and mass are always paired together.The interpretation of J 2 R 2 in terms of an object's shape requires first the choice of a shape model.For example, assuming a homogeneous oblate ellipsoid with semi-axes a = b ≥ c gives a c = b c = 5 J2R 2 c 2 + 1. Different shape models can match the observed J 2 R 2 even given a specific mass (see, e.g., Marchis et al. 2005b;Ragozzine et al. 2023) and the shape model influences the choice of R and thus J 2 .Unknown contact and compact binaries can also result in an apparent J 2 R 2 , as discussed below.Except for the largest TNOs, J 2 R 2 is dominated by the overall shape of the object and not the interior density distribution.
However, the inability to determine a precise shape or interior composition with J 2 R 2 should not detract from the fact that measuring a significant J 2 R 2 has important implications for the object.For example, different formation modes can result in different J 2 R 2 values with catastrophic collisions reaccumulating into nearly spherical objects while high angular momentum formation processes tending to form more non-spherical shapes.Of course, the strength of J 2 R 2 affects long-term spinorbit-tidal dynamics and evolution of the TNBs themselves (e.g., Correia 2018).The measurement of J 2 R 2 also provides mostly orthogonal constraints on shape properties, so that it improves the interpretation of light curves, occultations, and thermal measurements.
The rotation poles of TNOs have also proved to be difficult to constrain, although some success has been seen among the Centaurs (e.g., Tegler et al. 2005;Duffard et al. 2014).Until now, determination of shape has relied on both light curve and occultation studies.While light curves are powerful tools for understanding shape and spin, the shape and spin solutions produced are non-unique (Harris & Warner 2020).Additionally, light curve inversion requires observations of the target at a variety of aspects, which is implausible for most TNOs due to their centuries-long heliocentric orbital periods.On the other hand, occultations can directly observe the shape of a TNO in a projected plane (e.g., Elliot et al. 2010;Benedetti-Rossi et al. 2016;Ortiz et al. 2017), placing limits on the full 3-dimensional shape of a body when combined with a light curve.However, due to the random timing of the events and difficulty of observation, occultations cannot yet be used to systematically understand the shapes and spins of TNOs at a population level.
Detection and measurement of non-Keplerian effects thus opens the door to a deeper understanding of TNB systems.While Keplerian orbital fitting allows for determination of mass (and therefore density), non-Keplerian orbit fitting in theory allows for study of the three dimensional shape and spin orientation of the two system components.It is able to provide unique constraints on both the shapes of TNB components and their obliquities 1 (Ragozzine & Brown 2009;Vachier et al. 2012) which are sensitive tracers of formation and evolution (e.g., McKinnon et al. 2020).

Unknown Components
1 Throughout this paper, obliquity will refer to the inclination of the mutual orbit with respect to the primary body's equator.
Another area where non-Keplerian analysis has significant potential for unique insight is in the detection of unknown components whether in "contact binaries" or simply closely separated objects ("compact binaries").Many channels of evidence suggest that these are common in the TNO population, including light curve studies (Rabinowitz et al. 2019;Thirouin & Sheppard 2018, 2019;Showalter et al. 2021), occultations (Leiva et al. 2020), and imaging from within the TNO population as enabled by the New Horizons spacecraft (Weaver et al. 2022).Indeed, the only small TNO ever visited by a spacecraft, Arrokoth, was a contact binary (Stern et al. 2019;McKinnon et al. 2020).Discovering and characterizing contact or tight components in binaries provides unique insight into the angular momentum and mass budget of TNBs.
There is a huge range of parameter space where additional components would be too far away to be seen in light curves or occultations (≳3 primary radii away) but too close to be resolved in direct imaging (≲20 primary radii away for most TNBs).Non-Keplerian effects, on the other hand, are highly enhanced in both contact and unresolved binaries.
Both non-spherical shapes and additional components primarily cause orbital precession.Indeed, to quadrupole order, the gravitational potential of a closein, previously undetected satellite is effectively the same as unusually large gravitational harmonics.This means that exploration of non-Keplerian effects using a J 2 shape model is a good starting approximation for detecting additional components.Detecting anomalously large values for J 2 (e.g., J 2 ≳ 1 when assuming a triaxial shape model), may indicate that there are additional undetected system components.
Using the definitions of the gravitational harmonics from Yoder (1995), we find that the 'effective' J 2 and C 22 of a close-in satellite would be where q = m s /m p , m p and m s are the masses of the masses of the primary and (unresolved) secondary, respectively, and a s is the semi-major axis of the secondary's orbit.In this equation, we leave the combinations J 2 R 2 and C 22 R 2 since these are the physically meaningful parameters (see Equation 1).It is possible to understand these contact/compact binaries as one end of the continuum of shape modeling that converts the measurement of J 2 R 2 into a physical configuration of mass.
These contact and compact binaries describe one (or both) of the components of an existing TNB, converting binaries into hierarchical triples (or quadruples).These systems, only one of which was previously known (Lempo; Benecchi et al. 2010), are common outcomes of gravitational collapse simulations (Nesvornỳ et al. 2010;Robinson et al. 2020;Nesvornỳ et al. 2021), although they may also be able to form through capture mechanisms (Brunini & López 2020).Regardless of their origin, discovering additional Lempo-like systems will aid in understanding the formation of TNBs.

Expectations for Non-Keplerian Effects
While non-Keplerian fitting does have a variety of drawbacks and degeneracies (see further discussion in Ragozzine et al. (2023)), it is another powerful tool for understanding TNOs at a deeper level.Detecting non-Keplerian effects may allow for a new suite of unique measurements leading to a significant improvement in understanding TNBs, but the fraction of systems with detectable effects (with current or potentially future data) has been unknown.Non-Keplerian effects of multiple moons has been robustly detected in the Haumea system (Ragozzine & Brown 2009), but no J 2 precession has ever been robustly detected in any TNB system, despite attempts (e.g., Ragozzine & Brown 2009;Gourgeot et al. 2016).
There is good reason, however, to expect that non-Keplerian shape effects are common and detectable.Firstly, aspherical shapes are common among TNOs, as shown by the results of occultation observations (e.g., Elliot et al. 2010;Benedetti-Rossi et al. 2016;Ortiz et al. 2017).These shapes, assuming that TNOs are not all differentiated, should produce detectable non-Keplerian effects over the current observational baselines.Secondly, detection of J 2 is common in both main-belt asteroids and near-Earth asteroids (e.g., Marchis et al. 2005a,b;Fang et al. 2011Fang et al. , 2012;;Vachier et al. 2012;Beauvalet & Marchis 2014;Marchis et al. 2014).For example, Marchis et al. (2005b) unambiguously detected the J 2 of the asteroid (121) Hermione based on the orbit of its small companion, with observations spanning less than a single year.While TNBs generally have longer period orbits than asteroid binaries, several TNBs have relative astrometric measurements now reaching decade long timescales.The observational baselines and high quality data strongly suggest that non-Keplerian effects would be detectable in at least a few TNBs.
To test the detectability of non-Keplerian effects in TNBs, we created synthetic relative astrometry for all TNBs, assuming the primary and secondary were interacting triaxial bodies with realistic shapes.The syn-thetic astrometry had realistic uncertainties added and was made to simulate the quantity and quality of real TNB relative astrometry by using existing datasets as templates.We then fitted the astrometry with Keplerian orbits to find if a Keplerian model provided adequate fits to the data.In many cases, we found that our synthetic astrometry was significantly inconsistent with a simple Keplerian model.We found the inconsistencies were strongest among TNBs with large primaries, although many near-equal sized binaries also exhibited inconsistencies.These results imply that non-Keplerian shape effects should be detectable and relatively common among TNBs.
If it is indeed true that non-Keplerian effects are common in TNBs, it also raises questions about the validity of past results based on Keplerian orbital analyses.In most scenarios, Keplerian fitting should yield a similar orbit to a non-Keplerian fit, especially in the parameters not describing the orbit's orientation (e.g., total system mass, semi-major axis, and eccentricity).Some results, however, depend on the orientation of the orbits (e.g., Grundy et al. 2019), which slowly precess over time in a non-Keplerian framework.The systematic error introduced by this model misspecification needs to be evaluated to determine whether past research needs to be reevaluated.

Non-Keplerian Orbit Fitting
For all of its benefits, orbit fitting of non-Keplerian orbits can be a difficult task.While the potential of triaxial body, given in Equation 1, is relatively simple and could be integrated without much trouble2 , the actual dynamics of TNBs can be much more complicated.Not only can both a TNB's primary and secondary both have aspherical shapes with arbitrary rotation poles, but the gravitational interactions between the bodies cause a torque on those bodies, forcing precession of those rotation poles.The torques, which we call "back torques," can make the dynamics of a TNB system much more complicated, especially when the binary is near-equal mass.Many past studies -both observational and theoretical -of non-Keplerian motion in asteroid binaries have neglected the effects of back torques.While this assumption is valid for very small secondaries (like many asteroid binaries), TNBs are often near-equal mass.The complexity involved in TNB spin-orbit dynamics requires the use of a coupled spin-orbit integrator, which can simultaneously integrate the translational and rotational equations of motion.
For this reason, we developed a new spin-orbit integrator, named SPINNY, which is able to self-consistently model coupled spin-orbit motion of an arbitrary number of gravitationally interacting triaxial ellipsoids, including the effect of back torques on the spinning bodies.SPINNY makes up the core of MultiMoon, a package designed for non-Keplerian orbit fitting using state-of-the-art Bayesian techniques.Both SPINNY and MultiMoon are described in full detail in Ragozzine et al. (2023) and are publicly available at https:// github.com/dragozzine/multimoon,though we discuss the main ideas here.MultiMoon's advanced statistical techniques and orbital model, combined with about a decade worth of relative astrometry for TNBs and extensive parallel computations, now make it possible to conduct a systematic survey to search for non-Keplerian shape effects in TNBs.
Exploration of possible non-Keplerian parameters is generally a difficult task, requiring at least an order of magnitude more computational power than a simple Keplerian analysis.As such, the most interesting targets need to be identified to focus future work.Here, we perform a broad search for the best targets for more detailed non-Keplerian analyses.We focus our efforts on identifying TNBs with the most detectable and statistically significant non-Keplerian effects.

METHODS
In this work, for all orbit fitting procedures, we use MultiMoon.MultiMoon is a Python-based Bayesian orbit fitting package specifically designed to complete Keplerian and non-Keplerian orbit fitting of TNBs and other small solar system binaries (Ragozzine et al. 2023).For Keplerian orbit fitting, it uses the SpiceyPy (Annex et al. 2020), a Python implementation of SPICE (Acton Jr 1996), to analytically solve the two-body problem.For non-Keplerian orbits, MultiMoon uses SPINNY which numerically integrates the coupled spin and orbit equations of motion for an arbitrary number of interacting triaxial ellipsoids.These processes produce the position of the secondary relative to the primary.These are projected into the plane of the sky using ephemerides of the TNB system (relative to the Earth) from JPL Horizons, as queried by astroquery (Ginsburg et al. 2019).These ephemerides are corrected for light-time variations and other astrometric aberrations due to Earth's orbital velocity.Light-travel time and conversion to ecliptic coordinates is done as part of the data preparation.For more details, see Ragozzine et al. (2023) or the MultiMoon code itself at https://github.com/dragozzine/multimoon. The analysis herein uses a version of MultiMoon equivalent to Version 1.0.Since our fitting was completed, MultiMoon has been updated to more accurately model the spin dynamics of small secondaries.We have confirmed that this small change has no impact on any of the results presented in this work.
The Bayesian parameter inference at the heart of MultiMoon is performed by emcee (Foreman-Mackey et al. 2013, 2019), a Python-based Markov Chain Monte Carlo (MCMC) ensemble sampler.This allows us to robustly fit orbits in a Bayesian manner.MultiMoon uses the same weighted least squares orbit fitting likelihood function that is ubiquitous in TNB orbit fitting (e.g.Grundy et al. 2008Grundy et al. , 2009;;Ragozzine & Brown 2009;Kiss et al. 2019;Holler et al. 2021, etc.).MultiMoon differs from past orbit fitters, however, in that it does not use a downhill optimization technique, but instead uses ensemble MCMC techniques that are robust to local minima.This avoids problems that have occurred when using problematic simple minimization techniques (Beauvalet & Marchis 2014).Data weights in our likelihood functions are equal to the inverse square of the measurement uncertainty.
Briefly, the outputs from MultiMoon are a posterior probability distribution of inferred parameters given the observational data and priors (which are generally flat and uninformative).Based on the output posteriors, MultiMoon produces a variety of different output plots, including both diagnostic plots (e.g., walker trace plots, likelihood plots, etc.) and publication-ready plots (e.g., corner plots, Foreman-Mackey et al. 2016).For more discussion, see Ragozzine et al., submitted.

Data
For all orbital fits completed in this paper, we use publicly available relative astrometric data and orbital fits catalogued in Will Grundy's online TNO database3 , providing the relative astrometry and orbital solutions of 45 TNBs with solved or mirror ambiguous orbits.The list of all the TNBs in our analysis is given in Table 1; note that the meaning of some columns is explained in more detail below.With both relative astrometry and orbital solutions in hand, both are converted into ecliptic coordinates and used in our analysis as described above.All inputs and outputs from our fits are available online on Zenodo ( doi: 10.5281/zenodo.10620251).
Although the majority of our sample consists of the TNBs whose mutual orbits have been previously stud-ied in the literature, we exclude two systems, (136199) Eris-Dysnomia and (469514) 2003 QA 91 .The Eris-Dysnomia system, whose Keplerian orbit has been extensively studied in the past (Brown & Schaller 2007;Holler et al. 2021), is known to be > 6σ inconsistent with a single Keplerian orbit (Holler et al. 2021).Non-Keplerian orbit fitting (completed with MultiMoon) is discussed in Spencer et al., in prep. For 2003 QA 91 , the publicly available relative astrometry (which are compiled in Grundy et al. 2019) are clearly inconsistent with the given orbit solution.We believe that a typographical error is present in the publicly available astrometry.For these reasons, we exclude these systems from our analysis.
Most of the relative astrometry we use is derived from precise images from HST, although significant contributions are made from a variety of other ground-based observatories.Typical error bars on the astrometry from HST are less than 10 milliarcseconds, although they can even be sub-milliarcsecond when observations were taken with HST's High Resolution Camera.We also use many observations from Keck using its adaptive optics system, which typically give measurements with similar precision to HST.For some of the wider TNBs, observations can be obtained with smaller ground-based telescopes like Gemini, Magellan, and the Canada France Hawaii telescope.Given the lower resolution of these telescopes, typical uncertainties can be 100s of milliarcseconds.To get a sense of the typical precision of the observations used for our orbit fitting, we list the median uncertainty of the astrometric measurements in Table 1.
For all astrometric data, we assume that uncertainties are Gaussian and measurements are independent.Unfortunately, when combining various data sources (i.e., different telescopes and instruments), systematic errors can arise from cross-calibration issues, uncertainties on time varying plate solutions, and other issues.Past analyses have shown that these issues are not significant for the vast majority of TNOs (e.g., Grundy et al. 2008Grundy et al. , 2014Grundy et al. , 2019, etc.), etc.), most of which are not resolved (or barely resolved) in even the largest of telescopes.Further, a posteriori, we can use heuristics like reduced χ 2 to check whether our noise model (and orbit model) adequately describes the data.We find that systematics in our data are unlikely to affect our modeling significantly and add noise at (or below) the measurement uncertainties.Note-The list of 45 Trans-Neptunian Binaries which we study with both Keplerian and non-Keplerian models.This includes all TNBs (except two unusual cases) that had known orbits or mirror-ambiguous orbits.Rotation periods and primary radii used in our non-Keplerian fits (and their references) are also listed.The epoch of the fit is also listed; this is the time at which the orbital element parameters are set in non-Keplerian fits (which have time-varying orbital elements).Uncertainty refers to the typical uncertainty in the relative astrometry used to produce orbit fits.

Keplerian Orbital Fits
As a preliminary step to completing non-Keplerian fits, we first completed a round of Keplerian fits to validate our orbit fitting techniques and provide full Bayesian posteriors for the Keplerian orbits of TNBs.This allows us to analyze the quality of the Keplerian fits and provide a baseline for comparison when completing non-Keplerian fits.
Our Keplerian orbital model has seven parameters, including the usual six orbital elements (a, e, i, ω, Ω, and M) and total system mass (M sys ).All angles are referenced to the J2000 ecliptic plane (unlike many publications of these TNB orbits which use angles in the equatorial reference frame).We also consider a model with constant photocenter-barycenter offests which is discussed in Appendix A.
Our Keplerian fits were run with 100 walkers in the MCMC ensemble, with a 5000 step burn in, pruning of walkers significantly far away from the best parameter space (similar to Proudfoot & Ragozzine 2019;Hou et al. 2012), a 1000 step post-pruning burn in, and a 5000 step sample.We initialized walker positions for all parameters by drawing random, normally-distributed samples from the orbital solutions listed in the Grundy database, although we inflated uncertainties to allow for a broader search of parameter space.For systems with mirror ambiguous orbits, we ran two different fits to explore both orbital solutions.
We set uninformative priors for these fits, using uniform distributions for each parameter (see Ragozzine et al. for more details).Similar to other orbit fitting procedures, we compared our model to the relative astrometry using the χ 2 statistic (c.f.Grundy et al. 2019).After completion of the runs, we checked for convergence of the MCMC chains based on trace plots and the smoothness of the marginal and joint posterior distributions of each parameter.In a few cases where we could not clearly confirm the MCMC chains were converged, we reran the fits with more burn in and sampling steps, until we were confident the chains were converged.

Non-Keplerian Orbital Fits
Using MultiMoon's non-Keplerian orbit fitter, we completed non-Keplerian fits to all 45 TNB systems.Although MultiMoon is capable of fitting the gravitational harmonics of all bodies in TNBs system, for this analysis, we only fit the J 2 gravitational harmonic of the system's primary with an assumed size and rotation rate.This greatly reduces the complexity of the model as only 3 parameters are added (J 2 and two spin pole direction angles) rather than 14 (ellipsoid polar axis length c, J 2 , C 22 , two spin pole direction angles, a spin longitude angle, and a rotation rate, all for both objects).
Since non-Keplerian effects can break the usual mass degeneracy in Keplerian fits, we also allow both masses to float, adding one more parameter.In many TNB systems, the dominant source of non-Keplerian effects should come from the J 2 of the primary, but even if this is not the case, using a single J 2 is a good approximation to the modeling the system's "total" J 2 , as measured by the orbital precession.It also provides a way to model non-Keplerian orbital precession whatever its source, as discussed above and in Ragozzine et al. (2023).While neglecting C 22 can provide worse fits, this is only the case when a TNB is near a low-order spin orbit resonance.Given that most TNOs rotate in ∼10s of hours and TNB orbit periods are ≳10 days, low-order spinorbit resonances among our sample will be extremely rare.The only known TNB in our sample which is at spin-orbit resonance is Sila-Nunam.We discuss the consequence of neglecting C 22 for the Sila-Nunam system in Section 4. As our goal is identifying objects with the most statistically significant non-Keplerian effects, even if our assumptions do not hold, the improvement in orbital fit from the wrong model would still indicate the need for higher fidelity and more advanced fits.
Given these assumptions, our non-Keplerian orbital model has 11 free parameters, the six Keplerian orbital elements at a given epoch, the masses of the individual system components (M 1 and M 2 ), two angles describing the rotation pole of the primary at epoch (i sp and Ω sp , which are Euler angle representations of the spin pole), and ln J 2 R 2 .We opt to use the ln of J 2 R 2 to enable easier exploration over many orders of magnitude.Additionally, we use the combination J 2 R 2 instead of just J 2 so results can be interpreted with a variety of shape models.
In our MCMC ensembles, we used 100 walkers running for 26000 total steps, split between 15000 burn in steps, pruning of poorly performing walkers, 1000 postpruning burn in steps, and 10000 sampling steps.We set integration tolerance value at 10 −10 for all runs.This was chosen to balance integration quality/accuracy and computational expense.In all our testing, this tolerance level was found to be sufficient for our needs.
We initialized walker positions for mass and orbital elements identically to the Keplerian fits.In addition, we initialized walker positions for log J 2 R 2 by drawing random samples from a normal distribution centered at log (10R), where R is the object's estimated radius, a rough approximation of typical J 2 values due to nonspherical shapes over a large size range.Lastly, we initialized walker positions for the spin pole direction to be generally aligned with the orbit, within ∼10-30 • of perfectly aligned.For other required MultiMoon inputs (approximate radius and rotation period), we used values previously published in the literature (see Section 2.1 for more details).For objects with no known rotation period, we used a default value of 10 hours (similar to those found in Thirouin et al. 2014).For the few TNBs that did not have estimated radii published in the literature, we assumed a value based on an assumed albedo and the system's absolute magnitude.Information on all TNBs considered in our analysis, along with the input for approximate radius and rotation period, is located in Table 1.Our testing indicates that changing these assumed values does not substantially affect the results of our fits.Like the Keplerian orbit fits, mirror ambiguous TNBs had two orbit fits run to explore both orbital solutions.
Our priors on each parameter were uniform with the goal of being uninformative, except in a few cases.Firstly, we enforced a(1 − e) > q min , where q min was chosen to be equal to a factor of a few times the primary radii.This helped to reduce the chances of unphysical close encounters that were occasionally explored in extreme non-Keplerian integrations.We also force M 2 < M 1 .This practice reduces the effects of degeneracies between mass and J 2 .
Finally, ln J 2 R 2 < 15 was enforced to limit exploration in J 2 space.At large values of J 2 R 2 , some TNB systems can become unstable, placing this constraint prevents exploration of these unstable models.While this may prevent exploration of relevant parameter space for some systems, we found that this prior appropriately balances stability and our goal of broad exploration.In a small number cases, the prior had to be further reduced after confirming that large values of J 2 caused unphysical models4 .A uniform prior in ln J 2 R 2 is not meant to signify our actual prior knowledge of the J 2 distribution of TNBs, but rather to encourage exploration of this new parameter.We adopt as meaningful only those fits where the likelihood (calculated based on χ 2 as with Keplerian modeling) of a non-Keplerian fit strongly prefers a particular value of J 2 and thus is relatively independent of the prior.
As our goal of this project was to identify systems with the most statistically significant non-Keplerian effects, we did not aim for full convergence of our MCMC chains.Rather, we focused on achieving the best non-Keplerian fit possible in a fixed number of model evaluations.While this results in unconverged fits -so that the full posterior probability distributions are not expected to be completely accurate -the unconverged chains can still be informative, especially in light of our goal of identifying targets for future investigations.In essence, our runs are set up to determine which TNBs have the most easily detectable non-Keplerian effects.We refer to these potentially unconverged fits as "exploratory" and consider them valuable at the full catalog level of this analysis.
These fits and their implications were studied individually in detail by various co-authors as part of Brigham Young University's Physics 227 ("Solar System Astronomy") Class Project.

Non-Keplerian Orbital Fits for Borasisi-Pabu
To show that unconverged fits described above are able to identify targets worthy for future investigation, we completed a full non-Keplerian orbit analysis of the Borasisi-Pabu system, one of the more promising targets identified in our exploratory fits, as a proof-of-concept.In this fit, we relaxed several of the assumptions made in the non-Keplerian exploratory fits, most notably our assumptions regarding the C 22 of the primary.Our non-Keplerian orbital model has 13 parameters, six Keplerian orbital elements at epoch, the masses of both system components, the direction of the primary's rotation pole at epoch, the J 2 R 2 and C 22 R 2 of the primary, and the longitude of the primary's prolate axis at epoch (ω sp ).For Borasisi's rotation period we used 6.4 hours (Kern 2006) and a radius of 63 km (Vilenius et al. 2014).After our analysis was completed, Kecskeméthy et al. (2023) found a different light curve period for the combined system, we discuss the minor implications of this in Section 5.
For this full orbit fit, we ran MultiMoon with 980 walkers for 23500 steps, split between 15000 burn in steps, 1000 post-pruning burn in steps, and 7500 sampling steps.Integration tolerance was set to 10 −11 to produce the best quality of fits possible, with less regard for optimizing computational expense than the previous round of fits.Initial walker positions were drawn based on preliminary orbital fits and the results of our exploratory fits, but otherwise the fitting process was effectively identical to that above.We confirmed convergence of the MCMC chain by inspection of trace plots and marginal and joint posterior distributions of each parameter.
Our quadrupolar approximation cannot distinguish between models rotated by 180 degrees (e.g., it cannot distinguish the North pole from the South pole), resulting in a two-fold degeneracy.Since the degeneracy is fairly well-understood, we decide to fit only the prograde rotation solution, where the rotation axis is required to be inclined < 90 • relative to the binary's mutual orbit plane.

Evaluating Fit Quality
To ensure the quality of all of our fits (both Keplerian and non-Keplerian), we have thoroughly evaluated each fit using a variety of statistical techniques.Most simply, and easily output from MultiMoon is the calculation of reduced χ 2 , which is equivalent to the best fit parameter set's χ 2 per degree of freedom.For almost all of our fits, reduced χ 2 was ≲ 1, indicating good quality fits.For 8 of our Keplerian fits, however, we found statistically significant cases of elevated reduced χ 2 .Several of these systems have been previously identified as inconsistent with Keplerian orbits (e.g.Grundy et al. 2019).These systems all had significant, or nearly significant, non-Keplerian effects detected (see Section 4).
Another way we can analyze the quality of our orbit fits is to calculate the root-mean-square (RMS) residual of the fit.For our fits, we find that most of our fits have RMS residuals ≲ 15 milliarcseonds, indicating an excellent fit to the data.Compared with the typical uncertainties in the astrometry we fit to (see Table 1), we find that our fits are very robust.A few systems have much larger RMS residuals, but this is caused by large uncertainties in the data.As discussed in Section 2.1, the ultra-wide TNBs tend to have many observations taken with relatively imprecise ground-based observatories, which can result in large RMS residuals.Combined RMS residuals (with each dimension added in quadrature) are reported in all tables below.
Like RMS residuals, we also looked at the weighted RMS residuals.These residuals are essentially in units of observational error bars (as shown in Table 1).We find that all of our fits have weighted RMS residuals < 1.4, with most being < 1.This confirms the quality of our fits.
We also closely examined the residual plots for each fit (see Figures 2 and 7 for examples).These plots show the normalized residuals (residual divided by measurement uncertainty) and were examined to ensure no systematic trends were visible in the residuals (e.g., uncentered residuals, significantly larger residuals in one dimension, correlations among subsets of the data taken with different facilities, etc.).We found no serious issues with any of our fits in individual fits, and also found no problems at the ensemble level.These normalized residuals are available for both our Keplerian and non-Keplerian fits publicly available on Zenodo ( DOI: 10.5281/zenodo.10620251).
Another sign of the quality of our fits is their close match to orbit fits published in the literature.In Section 3, we show that our orbit fits closely reproduce fits in the literature.This close agreement is strong evidence that our orbit fits are of high quality.

KEPLERIAN FITS
The results of our Keplerian fits are contained in Table A1.It contains information regarding the posterior distributions of the mass and orbital elements of each TNB (median values with 1σ confidence intervals).In addition to values from the marginal (single variable) posterior distributions (as seen in the table), the full posterior distribution can be displayed as a corner plot.
As an example, we show the corner plot of the Keplerian posterior for (66652) Borasisi-Pabu in Figure 1 and a residual plot in Fgiure 2. The results for all 45 TNBs (including MCMC chains, diagnostic plots, statistical information, etc.) are publicly available5 .
For all mirror ambiguous TNBs, two orbit solutions are shown in Table A1, with the exception of 2002 WK 183 and 2000 OJ 67 .For both of these objects, despite attempting two orbit fits for the separate solutions, the mirror ambiguous solutions are so close together that the posteriors of each solution significantly overlap.Separating the two distributions in a statistically rigorous way is difficult due to the blended nature of the posterior.Note that the emcee algorithm at the heart of MultiMoon is not optimized for rigorous exploration of multi-modal posterior distributions.As such, we report the statistics of the blended distribution and urge caution when using the orbit fits for these objects.
As a test of MultiMoon's performance in completing Keplerian fits, the posterior distributions of the ensemble of fits can be compared to the best fits from the Grundy database.We calculated the z -score of every Grundy best fit for every parameter and every object, using our posterior distribution as the reference distribution.The distributions of these scores for each parameter is shown in Figure 3 using kernel density estimates.The z -score distributions clearly show that the Grundy best fits are consistent with the MultiMoon posteriors.Technically, since we are analyzing the same data with the same model, the agreement should be much better than drawing parameters from a random distribution.While not a perfect match, the fact that the vast majority of parameters are within 1-σ means that any systematic errors are minimally important compared to the statistical errors.We note that a similar analysis of an ensemble of Keplerian TNB orbit fits also found consistency with the Grundy results (Emelyanov & Drozdov 2020), further supporting that these results are robust to the analysis method.
Likewise, when comparing the Grundy database parameter uncertainties (which are determined via Monte Carlo techniques) to the MultiMoon Keplerian posterior distributions, we find excellent agreement.For all parameters (except argument of periapsis ω and mean anomaly M), the error bars are similarly sized for all the TNBs.For ω and M, our model is differently parameterized, where the Grundy database uses longitude of periapse (ϖ) and longitude at epoch (ϵ).This slightly different parameterization obscures a like-to-like comparison.
For those interested in precise understanding of the uncertainties in Keplerian fit parameters, our reporting of full posteriors provides an improvement to the state-of-the-art which approximates parameters as having a mean value and uncertainty.Posterior samples are also particularly useful for scheduling of follow-up observations.Each sample in the chain provides a predicted position of the binary components at any given time.Taken together, the ensemble of samples (or a sufficiently large subset of the ensemble) then provides a probability distribution for the predicted positions of the binary components, as further discussed below.
The results for the Keplerian model with constant photocenter-barycenter offsets are discussed in Appendix A. We find that the offsets presented are attributable to overfitting.We also find that they do not significantly affect the Keplerian results, except for possibly adjusting the eccentricity.
The agreement of our Keplerian outputs with previously published fits confirms that MultiMoon can effectively fit TNB orbits, despite the different parameterization and methods used.This validates MultiMoon's fitting procedures and techniques, lending confidence in our non-Keplerian fits presented in Sections 4 and 5.
Although our Keplerian fits only reproduce past work, this work allows us to identify possible systematic issues with TNB orbit fitting at the catalogue level.Furthermore, it allows us to understand if the uncertainties reported in astrometric data analysis are accurate.Given our high-quality fits, we find that Keplerian orbit fitting and current generation astrometric data analysis adequately account for systematic effects.We measure this using fit quality heuristics like reduced χ 2 and p-values.Aside from a few systems that are expected to have non-Keplerian orbits (e.g., Salacia, Orcus, etc.), our fits give reduced χ 2 ≲ 1, which implies that the reported astrometric uncertainties and orbit model describe the data well.Although this does not preclude all systematic effects, it significantly reduces the number of possible sources of error.Future investigations into sources of systematic error should use flexible noise models that can parameterize the sources of error present in the astrometric data.

NON-KEPLERIAN FITS
Aligned with our goal of identifying candidates for future non-Keplerian analysis, we present the non-Keplerian best fit model for each TNB in Table A2.The best fits in the table do not have uncertainties attached to them because our MCMC chains are not converged.In an MCMC framework, uncertainties are drawn from the 16th and 84th percentiles of the output chains.Without convergence, the resulting 16th and 84th percentiles are extremely unreliable estimations of the true uncertainties.Despite this, we still believe that the best fits are still useful for showing that non-Keplerian effects are detectable.To evaluate the detectability of non-Keplerian effects, we use a likelihood ratio test comparing the Keplerian and non-Keplerian orbit fits.The likelihood ratio test compares the goodness-of-fit of two nested models and quantifies the improvement in fit quality between the models.In our case, the Keplerian model is a subset of a non-Keplerian model where J 2 = 0.Under the likelihood ratio test if the null hypothesis (no currently detectable non-Keplerian effects) is supported by the data, the likelihood ratio, L K L N K ≈ 1, where L K and L N K are the Keplerian and non-Keplerian likelihoods, respectively.If the null hypothesis can be rejected (non-Keplerian effects are currently detectable), L K L N K ≪ 1. Traditionally, some threshold is chosen at which the improvement in fit is deemed great enough and the null hypothesis is formally rejected.However, the fits presented here are exploratory in nature and may not have reached the global maximum likelihood, as evidenced by our results in Section 5. Further analysis would only increase the improvement of non-Keplerian models, as these results act as a lower limit on the improvement of the fit.In keeping with our goal to identify targets for future analysis, we choose to give special attention to cases where L K L N K < 0.1 which are bolded in Table A2.Inspection of the posterior distributions (though not necessarily converged) also show that in these systems, J 2 R 2 = 0 is disfavored.Further work (i.e.converged fits) is needed to truly determine the true statistical significance of these improvements, although we confirm the improvement in the likelihood is more than expected from the addition of new parameters.
We do find that change in RMS residuals between our Keplerian and non-Keplerian fits are not strongly correlated with improvement in fit as measured by L K L N K .
Some exploration of possible reasons as to why this is the case showed that the result is due to heterogeneous uncertainties in the data.When instead using normalized/weighted residuals (as in Figure 2), we see the expected correlation between change in fit quality and change in weighted RMS residuals.The 8 TNBs satisfying our improvement threshold have non-Keplerian effects that are currently detectable.In this section, we will first discuss population-level trends identified in our fits, and then individually discuss the 8 TNBs we identify as prime targets for non-Keplerian analysis.
In addition to Table A2, we display the ensemble of non-Keplerian fits in Figure 4 where we plot the best fit J 2 value against the likelihood-ratio statistic.In this figure, we assume a triaxial shape model for easy comparison between our fits and to the literature.We emphasize that only the best fit orbit (or J 2 ) is shown, as our fits are not completely converged.As such, no uncertainties can be associated with these measurements.Likewise, we do not claim that these best fits are the global best fits.The results should not be treated as precise measurements of physical parameters.While they obviously hint at the TNBs true physical properties, these can only be used as starting points for future analyses.
Using our fits we can identify several broad trends.First, and most strikingly, TNBs with large primaries are more likely to have smaller best fit J 2 values than small primaries, in line with theoretical predictions and observations of small bodies across the solar system.Another striking trend is that larger objects are more likely to have improvements using a non-Keplerian model, with the three greatest improvements (Salacia, Orcus, and Gonggong) being among the largest TNBs in our sample (which excludes Pluto, Eris, Haumea, and Makemake).This is likely because TNBs with large primaries tend to have shorter orbital periods (and thus stronger and faster non-Keplerian effects) in addition to having more high quality data.
Another interesting discovery is a group of objects with orbit fit improvements at relatively high J 2 (Borasisi, Altjira, 1999 RT 214 ).The high J 2 indicated by our exploratory fits potentially hint that one of the system components may be a contact binary.While this is an unusual configuration, among asteroid binaries several primaries have been found with extreme values for J 2 , perhaps most notably Kleopatra, with J 2 = 0.765 (Brož et al. 2021).Their prevalence in our findings, however, may be attributable to selection bias.Since these systems may have large J 2 values, their non-Keplerian effects would be far easier to detect, when compared to objects with similar data quality, quantity, and observa- The best fit J2 presented here is found by assuming a triaxial shape and dividing our measured J2R 2 by the resulting body's volumetric radius squared.The ratio L K L N K compares the Keplerian models presented in Section 3 to the best exploratory fit presented in this section.If the likelihood ratio is ≪ 1, the null hypothesis (no detectable non-Keplerian effects) could be rejected, and non-Keplerian effects are indeed detected.In general, throughout this paper, we adopt a threshold of L K L N K < 0.1 to indicate systems which have orbit fit improvements indicative of non-Keplerian motion.We have labeled all of our systems with improvements below our threshold, as well as several systems with marginal improvements.Asterisks in labels indicate either a prograde (*) or retrograde (**) orbital solution for mirror ambiguous TNBs.Color corresponds to the radius of the TNB's primary.We find that many TNBs are far better explained by a non-Keplerian model, with many systems showing moderate improvements.We also find a wide range of acceptable J2 values, ranging from oblate spheroids to unresolved binaries.Future work and observations should focus on confirming the significance of our detections and extending the observational baselines of systems with moderate improvements in fit quality.We urge caution in using the values of J2 in this plot as there may be significant uncertainties that our exploratory fits are unable to fully constrain/measure.tional baseline.As such, these systems may not be as common as would be suggested by our findings.
Similar to a possible discovery of contact binaries, our results show three TNBs with unresolved binary-like J 2 values (2005 EO 304 , Teharonhiawako ("QT297"), and 2006 BR 284 ), although none of the fits for these bodies cross our improvement threshold.As discussed above, hierarchical triple systems may be sensitive tracers of planetesimal formation in the early solar system, making a discovery of additional hierarchical triple systems an important goal.In fact, Nesvornỳ et al. (2021) conclude that Lempo-like triple systems "should be found in the Kuiper Belt when observations reach the threshold sensitivity."More observations of these bodies, and subsequent orbit reanalysis, are required to find if these systems can be confirmed as having detectable non-Keplerian effects.Preliminary investigations suggest that these hierarchical triples are not resolvable with imaging.
Another possible explanation for the large J 2 values found for these three objects is that these TNBs are strongly affected by the Sun's gravitational influence, with this influence manifesting as a large measured J 2 .All three are ultra-wide binaries, a class of TNBs with extremely large separation.Ultra-wides are most affected by the Solar tide due to their extremely long orbital periods.Future modeling efforts should aim to include the effects of the Solar tide to more fully model all gravitational dynamics at play.
Our ensemble of exploratory fits is also able to identify systematic errors which may stem from the use of Keplerian fits in a variety of past analyses.To do this, it is most useful to compare our non-Keplerian best fits with the full Keplerian posteriors discussed in Section 3. Using the same methods as the comparisons to the Grundy best fits, we compare our fits in Figure 5, using our Keplerian posteriors as the reference distribution.As can be seen, our non-Keplerian best fits are somewhat different to our Keplerian posteriors, especially for ω and Ω.This is expected for non-Keplerian analyses since the addition of J 2 shape effects allows for precession of these orbital angles.The absolute difference in these angles (at epoch), which can be several standard deviations from the Keplerian fit, are usually only a few degrees.
Notably, the masses, semi-major axes, and eccentricities found in our non-Keplerian best fits are consistent with the Keplerian posteriors.This shows that systematic errors in orbital orientations can occur when assuming Keplerian orbits, but that masses, semi-major axes, periods, eccentricities, and inclinations are not significantly affected.Thankfully, the large number of analyses that rely on Keplerian fits are unaffected by systematic errors induced by the exclusion of non-Keplerian effects.

Identified Targets
Here, we discuss in turn each of the eight TNBs identified as having detectable non-Keplerian effects.

(120347) Salacia-Actaea
Salacia-Actaea is one of the largest and most massive TNBs in our sample.Our fits used 14 individual observations over a 10 year time span.Detection of non-Keplerian effects is expected as previous analyses showed that observations were 3.7σ inconsistent with a Keplerian orbit (Grundy et al. 2019).
Our analysis had a best fit J 2 = 0.0195, which corresponds to a slightly oblate spheroid, but fits with similar likelihoods (although slightly worse) were different by a factor of a few.This large range of possible values for J 2 is expected for systems where only apsidal or nodal precession is detected.In this case, Salacia's J 2 is degenerate with its obliquity with respect to the mutual orbit.Since Salacia-Actaea's mutual orbit is nearly circular (e = 0.0062 +0.0031 −0.0027 , Table A1) the non-Keplerian effects present in the system are likely to be nodal precession, since apsidal precession is difficult to detect at low e.Detection of nodal precession implies that Salacia's rotation pole is misaligned with the binary's mutual In this figure, we compare the non-Keplerian best fits found in our exploratory fits to the Keplerian posteriors presented in Section 3. We find that significant deviations are present in the argument of periapsis (ω) and longitude of ascending node (Ω) when assuming a non-Keplerian model.This is expected since this model allows for precession of these two angles.The distributions are still centered at zero presumably because we are comparing the non-Keplerian angles at an epoch typically near the middle of the observations and because many systems show insignificant non-Keplerian effects.Notably, however, the system mass, semi-major axis, eccentricity, and inclination are still consistent with Keplerian fits, implying that systematic errors due to non-Keplerian effects are small in most previous analyses.
orbit, but the degree of misalignment is not possible to determine without further analysis.Future observations of the Salacia-Actaea system may be able to detect apsidal precession, allowing for a more full determination of the system's properties.Additionally, constraints provided by Salacia's low amplitude light curve (Thirouin et al. 2014) and future occultations may aid in this effort.
Our analysis had a best fit J 2 = 0.016, corresponding to a slightly oblate spheroid.Similar to Salacia-Actaea, due to Orcus-Vanth's nearly circular orbit only nodal precession was detected, introducing large uncertainties on our J 2 and obliquity measurements.
In the future, observations of Orcus-Vanth may help understand the system better, however, our prelimi-nary results show that detection of apsidal precession may be an order of magnitude more difficult than for Salacia-Actaea.Understanding the shape of Orcus (or Vanth) will probably require additional multi-chord observations of stellar occultations, similar to the 2017 occultation of Vanth (Sickafoose et al. 2019).

(225088) Gonggong-Xiangliu
Gonggong-Xiangliu, the most massive binary in our sample, is the most surprising target with detection of non-Keplerian shape effects.Our fits only used 6 individual observations over about 8 years.The dearth of data, especially when compared to other TNBs, is due to the discovery of the moon in 2016 (Kiss et al. 2017), although subsequent reanalysis found the moon in 2009 and 2010 HST images.Our detection is only robust when considering the prograde orbit solution.
Our analysis had a best fit of J 2 = 0.084, corresponding to a fairly oblate spheroid.The high J 2 we find is quite unexpected, and may possibly contradict light curve measurements (Pál et al. 2016).In our analysis, we believe we have detected both apsidal and nodal precession, allowing us to place constraints on Gonggong's obliquity.Our best fit obliquity is 116.4 • , however, our fits cannot distinguish between north and south pole.By folding the obliquities of our best fit solutions (the best fit and those with similar likelihoods) we find obliquities between 55 • and 65 • are favored.The large obliquity is surprising, but the binary's mutual orbit is also quite eccentric (e = 0.2852 +0.0086 −0.0079 ), possibly implying a complicated history of unusual tidal pumping (Kiss et al. 2019;Arakawa et al. 2021).We note that Gonggong's density has significant uncertainty that can be resolved when the obliquity of the satellite is clearly detected (Kiss et al. 2019), making Gonggong-Xiangliu an extremely promising observational target.
Based on the small dataset, another possible explanation for this detection is overfitting by our model.With small amounts of data, model fitting can become more susceptible to overfitting.In general, 5 or more observations are required to fit a Keplerian orbit (Grundy et al. 2008), so it seems suspicious that we find such high confidence detections of non-Keplerian effects, especially with no previously detected deviation from a Keplerian orbits (like Salacia and Orcus).As such, we do not place much confidence in this detection and consider the detection potentially spurious.
To better understand Gonggong-Xiangliu's unexpected eccentricity and obliquity, resolve the mirror ambiguity, and/or understand issues with overfitting, additional observations need to be taken of this interesting system.The large difference in brightness between Gonggong and Xiangliu make these observations difficult without using the most advanced telescopes (e.g., HST, JWST, or possibly Keck).

(66652) Borasisi-Pabu
Borasisi-Pabu, a small near-equal size binary in the Cold Classical Kuiper belt, was fit using 9 observations taken over about 8 years.No previous indications of poorly fitting Keplerian orbits have been reported in the literature, but our Keplerian analysis, which was consistent with past analyses, showed a Keplerian fit with a reduced χ 2 of 1.9, indicating a somewhat poor fit to the data.
We find a best fit J 2 = 0.279, consistent with an extremely elongated object, possibly similar to a contact binary.We also found an obliquity of ∼50 • for Borasisi, after detecting both apsidal and nodal precession.For Borasisi, the detectability of apsidal precession was enhanced by the mutual orbit's large eccentricity (e = 0.4696 +0.0017 −0.0017 ).We discuss Borasisi-Pabu in detail in Section 5 and present fully converged, high quality non-Keplerian fits.

(148780) Altjira
Altjira (and its unnamed secondary), a small equalsize binary just outside the Cold Classical belt, was fit with 8 observations taken over about 4 years.No past analysis indicated poor quality Keplerian fits, but our Keplerian analysis indicated a marginally bad fit (reduced χ 2 of 1.4).
Our analysis resulted in a best fit J 2 = 0.479, roughly consistent with a contact binary.The measured J 2 may also be caused by an unresolved system component, possibly making this a hierarchical triple system.We are able to constrain Altjira's obliquity to be between ∼ 30 • and 60 • by detecting both apsidal and nodal precession.
With V = 23 mag, little more is known about this system, with no known light curve or occultation measurements.Any additional constraints on this system will contribute to understanding the configuration of this system.

(174567) Varda-Ilmarë
Varda-Ilmarë is a somewhat high mass TNB; however, unlike the other high mass TNBs with detected non-Keplerian motion, the primary-to-secondary brightness ratio is far more equal, putting it in a class that is somewhat intermediate between the low-mass (with near-equal masses) and high-mass (large primary-tosecondary mass ratio) TNBs.Our analysis used 12 observations taken over about 4 years.Previous Keplerian analyses did not indicate a significant deviation from Keplerian motion, but our Keplerian analysis showed that Keplerian orbits were somewhat poor quality fits for both the prograde and retrograde orbital solutions at the ∼1σ level.Like with Gonggong-Xiangliu, our detection of non-Keplerian effects is significant for only the prograde orbital solution, although the retrograde solution almost reaches our improvement threshold.
Our fits find a best fit J 2 = 0.07, consistent with a significantly flattened body.As with other large TNBs with near-circular orbits, we have only robustly detected nodal precession, leaving a degeneracy between J 2 and obliquity.. Varda-Ilmarë is somewhat unique in having a mirror ambiguous orbit given the large quantity of high-quality data.This is primarily caused by the system's orbit which is viewed nearly face-on.Understanding Varda's role as a possible transitional system, may enable better understanding of the internal structure and composition of TNOs. 4.1.7. 1999 RT214 1999 RT 214 , a small near-equal size TNB in the Cold Classical belt, was fit using 6 observations taken over about 10 years.Past analysis did not show a significant deviation from a Keplerian orbit, but our Keplerian analysis indicated an extremely poor fit (reduced χ 2 of 3.2), although with so few data points, such a result is possible ∼1% of the time.With this amount of data, a detection of non-Keplerian motion is surprising, but the long observational baseline allows a fairly robust detection.
Our best fit analysis for this system gave J 2 = 0.436, corresponding to a contact binary-like shape.Another possible explanation is that a system component is an unresolved binary, possibly implying that 1999 RT 214 is a hierarchical triple.A wide range of permissible obliquities are allowed by our exploratory fits, with a range of ∼ 30−55 • being most likely.Our fits imply that both nodal and apsidal precession are visible.
Additional observations are needed to confirm our detection and to further characterize this system.These will have to be completed by HST, JWST, or large ground-based observatories since V = 24 mag for the system.Combined with new observations, high quality non-Keplerian fits may be able to enable a better understanding of this system.

(79360) Sila-Nunam
Sila-Nunam, the largest near-equal mass binary in the Cold Classical belt, was fit using 20 individual observations taken over 11 years.Uniquely, six of those observations are mutual events, where one system component passes in front of the other, from our view on Earth, producing a measurable drop in light.No previous modeling, including our Keplerian fits, hinted that non-Keplerian effects were detectable in this system, although previous measurements of Sila-Nunam's combined light curve showed that one, or both, components were flattened by ∼12% (Rabinowitz et al. 2014).
Our fits find a best fit J 2 = 0.199, consistent with a considerably flattened body.In this case, since both components are possibly flattened, it is possible that both bodies contribute significantly to the overall system's J 2 .Roughly speaking, a system where both bodies have a J 2 ∼ 0.1 would match our results.Our exploratory fits cannot confidently determine whether nodal or apsidal precession (or both) are detected with high confidence, especially since the geometry of the system introduces several degeneracies in the model parameters.For example, ω and M had solutions that were good fits 180 • away from our best fit.These issues resulted in quality issues in our fits.Future non-Keplerian fits may be helped by introducing priors or reparameterizations which reduce the effects of degeneracy.
In a tidally evolved system where the primary and secondary are locked in a spin-orbit resonance, like Sila-Nunam is thought to be, C 22 can play a prominant role in the system's dynamics (Proudfoot & Ragozzine 2021).This tidally evolved state is at odds, however, with the small, but potentially non-zero, eccentricity measured in a variety of analyses (Grundy et al. 2012;Benecchi et al. 2014), including both our Keplerian and non-Keplerian fits (e = 0.0158 +0.0167 −0.0098 ).This may be the result of relatively recent excitation which has yet to be damped out, or may be caused by unusual tidal dynamics.
Another explanation for the improvement in fit is overfitting.Sila-Nunam is unique in its use of mutual event data.Using mutual event data to aid in orbit fitting is unusual, but has been used in previous studies of this system without issue (Benecchi et al. 2014).We do, however, acknowledge that this data may introduce unforeseen systematic errors into our analysis, which may somehow result in better fits with a non-Keplerian model.Additionally, we found that a Keplerian model was a reasonably adequate fit to the data (reduced χ 2 ≈ 1), calling into question the need for a more complicated model.
The Sila-Nunam system is a prime target for future non-Keplerian analysis, especially with more complicated models that include the effects of C 22 .This future analysis should also rigorously test if any systematic effects are introduced by the use of mutual event data and test whether overfitting may effect any resulting non-Keplerian fits.Given the mutual events and large set of photometric data, Sila-Nunam an ideal target for a thorough orbital analysis.

Non-Keplerian Fitting of Other TNBs
While the goal of this work is to identify the targets with the most detectable non-Keplerian shape effects, (and considerable work needs to be done to complete converged, high quality fits for these systems) other TNBs should not be ignored.Although the current data is not able to robustly detect non-Keplerian shape effects in systems not meeting our improvement threshold, full non-Keplerian fits can still be used to place upper limits on the presence of those effects.These upper limits may, in some cases, provide valuable constraints on the shapes of TNB components.
To illustrate this, we consider our non-Keplerian fits to the 2002 UX 25 system.While our analysis only shows a slight improvement with the addition of a J 2 and rotation pole, our fits show a rough upper limit on J 2 of ∼ 0.05.While this rough estimate of an upper limit cannot be fully verified without converged non-Keplerian fits, it requires that 2002 UX 25 must have a remarkably spherical shape.When combined with its unusually small density (0.82 ± 0.11 g cm −3 , Brown 2013), 2002 UX25 proves to be an unusual object and is possibly a touchstone for TNO formation (e.g.Bierson & Nimmo 2019).The low density suggests a significant porosity (potentially both macroporosity and microporosity) which seems unusual given its large size (>300 km radius), though central pressures are still relatively small as discussed by Grundy et al. (2019).Additional formation modeling is needed to explain the presence of a small (collisionally-formed?) eccentric satellite around a large nearly-spherical rubble-pile.

Future Observations
Based on our results, we believe that future observations should focus on those TNBs with detectable non-Keplerian effects (those described in the last section) and, possibly more importantly, TNBs that are near our improvement threshold.These TNBs show the most promise for future detections of non-Keplerian effects and provide the most efficient route for constraining the shapes of the largest number of objects.As observational baselines are extended, non-Keplerian effects become more detectable and uncertainties on the fits shrink.Observations can also be targeted at certain times when the uncertainty of the predicted position is largest, providing the most constraining observations (see Section 5 for more details).
In addition to high resolution imaging to resolve binary components, other types of observations, like light curves or occultations, still provide important orthogonal channels of information on the shapes and spins of TNBs.The constraints given by these observations are able to be meaningfully used as additional data sources in future more complicated models of TNB orbits.For example, the axial precession of a TNB primary due to back torques from the secondary, may be detectable in certain scenarios, which may enable a direct measurement of the primary-secondary mass ratio.Likewise, precession of the axes of TNO moons may provide information about the shapes and rotation poles of those moons (as suggested by Hastings et al. 2016).

BORASISI-PABU
As a proof-of-concept we have performed detailed non-Keplerian fits to the orbit of Borasisi-Pabu, accounting for non-Keplerian effects due to both J 2 and C 22 .The marginal posterior distributions for each parameter are displayed in Table 2 and joint posterior distributions are shown as a corner plot in Figure 6.We also show a residual plot in Figure 7.
Using the framework introduced in Section 4, we can assess the statistical significance of our non-Keplerian fit, when compared to the Keplerian fit from Section 3 (see also Table A1 and Figure 1).As before, the Keplerian model we use is a subset of our non-Keplerian model (where J 2 and C 22 are 0).Using a likelihood-ratio test, we find that L K L N K = 8.34 × 10 −3 , giving confidence that the non-Keplerian model is statistically significant.This likelihood ratio statistic is somewhat smaller than that found in our exploratory fits (compare with Table A2).Both Table 2 and Figure 6 show that a Keplerian orbit (J 2 , C 22 = 0) for the mutual orbit of Borasisi-Pabu is strongly excluded, confirming that a non-Keplerian model is preferred.A non-Keplerian orbit model is also indicated by the poor quality Keplerian fit.Using the Keplerian fit described in Section 3, we find that there is a ∼3% chance that random chance would result in a fit that is as bad (or worse) than the best fit we found.While not conclusive, this improvement merits further investigation.
To test how robust these results are to our input values (rotation period of Borasisi, c-axes of Borasisi and Pabu, integration tolerance), we also performed many additional non-Keplerian fits with these values substantially changed.These tests produced posteriors that are nearly identical to the one presented here, except when integration tolerance was set to a too high level.These tests showed that the choices we made for input values  .Instead of individually plotting the primary and secondary masses, which are completely degenerate, the system mass (Msys) is plotted to better show parameter correlations.Additionally, to allow for comparison with the literature, we show J2 and C22 rather than J2R 2 and C22R 2 , where we take the shape to be a triaxial ellipsoid with R as the volumetric radius.Of particular interest is the marginal distribution for J2 that clearly shows a high value of J2 is preferred, with J2 = 0 (equivalent to Keplerian motion) being excluded with a significance > 3σ.−1.87 Note-All fitted angles are relative to the J2000 ecliptic plane on JD 2451900.0 (2000 December 21 12:00 UT).Assumed c-axes for primary and secondary are 63 km and 54 km, respectively (Vilenius et al. 2014).We use a rotation period of 6.4 hours (Kern 2006).Our fit presented here is only to the prograde obliquity solution, another solution with obliquity of ∼ 135 • also exists.Our fitted parameters J 2 R 2 and C 22 R 2 are presented here as J 2 and C 22 to make comparison with the literature easier; in this transformation we take R to be volumetric radius.As mentioned in the text, we are not able to break the primary-secondary mass degeneracy, so we urge caution in using our fits for the masses.The best fit in the MultiMoon non-Keplerian fit corresponded to a χ 2 of ∼11.5 (compared to 21 for the Keplerian fit).We find that using our best fit Keplerian case, for this fit we find L K L N K = 8.34 × 10 −3 .RMS residuals are 2 milliarcseconds in both longitude and latitude.did not affect the fit presented here, and validate the inputs used in our exploratory fits.
In addition to validating input parameters, we also ran a series of recovery tests on MultiMoon.These tests used MultiMoon to fit synthetically generated data (with artificially added uncertainties) using the best fit solution in Table 2, allowing us to determine if MultiMoon could recover the input parameters.These tests returned similar fits to the one presented in this section, showing similar uncertainties and degeneracies, justifying the use of MultiMoon given our current data and mitigating the risk of overfitting for this system.
For most parameters, there is agreement with both our Keplerian solution and the Grundy best fit, at the ∼1σ level, except for Ω, which shows a strong disagreement with Keplerian fits.This is an expected result of non-Keplerian fitting as the orbital model now allows for precession of the node.Our posterior shows J 2 = 0.4446 +0.1192  −0.0928 implying an extremely non-spherical shape of Borasisi.The spin pole direction is nearly perpendicular to the ecliptic.
Our fits show no preference for a certain value of C 22 , as can be seen in Table 2 and the posterior distribution in Figure 6.The posterior for C 22 is purely a reflection of the prior limiting C 22 ≤ 1 2 J 2 R 2 .Since C 22 is not detected, the corresponding spin longitude parameter, ω sp , is unconstrained, with a uniform distribution between 0-360 • .This is expected since the mutual orbit is much longer than Borasisi's rotation period (∼46 days and ∼20 hours, respectively).In future analyses, exclu-sion of C 22 may be appropriate, when the orbit to spin period ratio is large.This may reduce the computational expense of the fitting process.
In addition to constraining Borasisi's J 2 , and by proxy, its shape, we can also place constraints on Borasisi's obliquity w.r.t. the mutual binary orbit.Our analysis shows an obliquity of ε sp,orb = 45.04 +6.65  −4.66 .This fit only reflects the prograde obliquity solution, where Borasisi's rotation pole and the binary's mutual orbit lie in the same hemisphere.Another obliquity solution (with ε retrograde = 180 • − ε prograde ) exists and provides a similar fit to the data.Distinguishing between these solutions may be possible with additional relative astrometry, depending on how the system's angular momentum is partitioned.In a future work, we will explore why two solutions exist and when they are distinguishable.
Although we do present the individual mass posteriors in Table 2, these are highly degenerate in the current fits and should not be used individually.Prospects for breaking the mass degeneracy will be examined in a future work.
The match between this full non-Keplerian fit and our exploratory fits supports our conclusion that the exploratory fits have identified promising targets for future non-Keplerian analysis and additional observations.5.1.Discussion on the inferred J 2 of Borasisi Based on the large value of Borasisi's J 2 from our analysis, several interpretations of the results are possible.First, and most simply, Borasisi is extremely flattened, and probably elongated, with a shape most likely resembling a contact binary.Contact binaries are common among the TNOs (Thirouin & Sheppard 2018, 2019;Showalter et al. 2021) and can be part of TNB systems (Rabinowitz et al. 2019).With current constraints on the J 2 , it is difficult to say more about Borasisi's shape, but with narrowed uncertainties, detailed shape modeling may be able to distinguish between snowman (like Arrokoth) or peanut contact binaries (like Kleopatra) (as in Marchis et al. 2005b).
Our measurement of Borasisi's obliquity (∼ 45 • ) is somewhat surprising, as most asteroid binaries have orbits fairly well-aligned with the rotation pole of the primary.TNBs, however, are formed differently and exhibit properties not seen in asteroid binaries.The closest comparison to Borasisi-Pabu is Manwë-Thorondor; Manwë is possibly a contact binary inclined to its mutual orbit by 27 • (Rabinowitz et al. 2019).Both are Cold Classicals.Combined with a similar mutual orbit eccentricity, total system mass, and separation, the Manwë-Thorondor system is strikingly similar to the Borasisi-Pabu system.Indeed, if we calculate Manwë's J 2 (based on the shape model from Rabinowitz et al. (2019)), we find a value very similar to Borasisi's.These objects may be part of a larger population of contact (or compact) binaries with oblique satellites.
The study of this population of objects should focus on explaining how these objects formed and evolved over the history of the solar system.One possible hypothesis is that these systems formed with relatively circular and coplanar orbits.Then, due to their somewhat loose separations, close encounters with other TNOs (or TNBs) may have excited their mutual orbits (similar to Campbell et al. 2022).Tidal dissipation is unlikely to be effective at damping eccentricity and inclination in these systems due to the system's large separation and J 2 (Porter & Grundy 2012), so any past excitation would remain visible in their orbits today.Alternatively, these systems may have formed with eccentric and inclined orbits.More advanced simulations of TNB formation, especially of gravitational collapse of pebble clouds, need to be completed before distinguishing between these (or other) possibilities.
One shortcoming of our model, that may affect our interpretation, is the assumption that Pabu's shape does not contribute to the system dynamics.This is likely an oversimplification given Pabu's large relative size.However, assuming that Borasisi and Pabu have the same obliquities (w.r.t. the mutual orbit) and shapes, we can calculate the required J 2 for each body to match our median J 2 R 2 , assuming the volumetric radius ratio of Borasisi and Pabu is ∼0.8 (Vilenius et al. 2014).We find that each body would have J 2 ∼ 0.3, still a shape that is consistent with a contact binary (see Figure 9 of Marchis et al. 2005b).In the case where the obliquities are not aligned, J 2 for each object would have to be even higher.This still implies that at least one system component must be a contact or compact binary.In the future, if more data becomes available for the Borasisi-Pabu system, a more complicated model of the system's dynamics including quadrupole gravity for Pabu may be warranted.
A second possible interpretation of the large J 2 value found in our analysis is that the Borasisi-Pabu system hosts another unresolved component.The possible properties of this unresolved component, assuming all three components are spheres, can be explored using Equation 4. In Figure 8, we plot the possible configurations (in terms of mass ratio and separation) of the unresolved component based on our median J 2 R 2 and its 1, 2, and 3σ uncertainties.It clearly shows that a wide range of allowed system configurations exist.While some are potentially resolvable with HST's Wide Field Camera 3 (WFC3), it is unlikely that re-solvable satellites (with separations of ≳2000 km) would be dynamically stable for long times with Pabu (periapse of 2400 km).However, unresolvable near-equalmass inner binaries with separations of a few hundred kilometers are quite possible.We also draw comparisons with other known tight binaries 2011 JY 31 , 2014 WC 510 , and Lempo-Hiisi, the inner binary in the Lempo triple system.These three comparisons have a range of masses, with Lempo being much more massive than Borasisi-Pabu and 2011 JY 31 being much less massive, but they illustrate that an unresolved system component is a plausible explanation for the detected non-Keplerian effects.
Under the unresolved binary interpretation, the obliquity measured in our analysis is a proxy for the inclination of the unresolved object's orbit w.r.t. the orbit of Pabu.An inclination of ∼ 45 • is above the critical inclination where the Kozai-Lidov mechanism can be active (i crit ≈ 39.2 • ), but if the central components are significantly aspherical, Kozai-Lidov cycling can be avoided.The short-and long-term dynamics of triple systems are extremely interesting and complex, but we leave a full investigation of them to future work.
Assessing the relative probability of each of our interpretations (contact binary and unresolved component) is impossible at the 2-quadrupole level.However, additional 3-body analysis with MultiMoon may allow us to infer the orbit of an undetected system component since it would induce somewhat different non-Keplerian motion.Additionally, other observations of light curves, thermal emission, or occultations may enable us to distinguish between these two interpretations.
After completing all the analysis discussed in this Section, a new study of TNO light curve found a different light curve period for Borasisi of 19.868 ± 0.032 hours with an amplitude of ∆m = 0.216±0.057(Kecskeméthy et al. 2023).However, our testing showed that our results were robust to different primary rotation periods, so we do not include this new rotation period in our analysis.The fairly high amplitude light curve, however, is generally compatible with our results.

Possible Systematics
Although our analysis is promising, there is still a distinct possibility that our results could be due to systematic errors.When comparing Figures 2 and 7, it seems that the improvement in orbit fit is driven by just one data point near the end of the dataset.This is expected from a true non-Keplerian orbit, as non-Keplerian deviations grow linearly in time, but could also indicate unknown systematics.The next data point seems to be consistent with a Keplerian orbit, further Borasisi radius Median J 2 1¾ limit 2¾ limit 3¾ limit Figure 8. Constraints on the mass-ratio and separation of a putative unresolved binary in the Borasisi system.The thick black line shows the mass ratio and separation given by the median value of the inferred J2 implied by equation 4. Colored contours indicate the 1, 2, and 3σ ranges from the J2 posterior.Gray dashed lines show several separations for reference.Here, we adopt a radius of 63 km for Borasisi, but this could be smaller, especially at mass ratios close to 1.The separation equivalent to a single Wide Field Camera 3 (WFC3) pixel on HST is calculated based on Borasisi-Pabu's typical distance from Earth and the instrument's 0.04" px −1 pixel scale.The separation of Pabu's orbit is shown in Table 2.The black points show other close binaries known in our solar system.Data for Lempo-Hiisi is from Ragozzine et al. (2023), 2014 WC510 is from Leiva et al. (2020) and assumes equal densities and albedos, and for 2011 JY31 is from Weaver et al. (2022) and assumes equal masses.The J2 found in our analysis is consistent with an unresolved binary interpretation, as the separation of the binary would be less than a single HST pixel across a wide range of mass ratios.
raising concern for systematics.However, upon examination, that observation was taken near periapse, when the non-Keplerian signal would be the smallest.
To check for possible systematics, we completed a fit without the second-to-last data point (from ∼2007).In this fit, we find that L K L N K ∼ 0.3, implying that the non-Keplerian fit shows only modest improvement over a Keplerian fit.Given the lack of data (just 9 observations), this is expected behavior for a true non-Keplerian orbit, but it still does not conclusively rule out possible systematics.At this time it is impossible to conclusively rule against (or for) systematic errors.We do point out however, that our non-Keplerian model is able to match the observations extremely well, which would be unlikely to occur given randomly distributed systematic errors.Additionally, our model outputs are physically realistic and similar to other known binary systems containing a contact binary (e.g., Manwë).Understanding if systematic errors can masquerade as non-Keplerian effects for the Borasisi-Pabu system will require additional observations of the system.We note that Borasisi, Sila, 1999 RT214, and Altjira are targets for HST Program 17206 (PI: Proudfoot) where two additional observations will be made to improve our ability to interpret non-Keplerian motion and systematic uncertainties.
In future investigations, sources of systematic errors should be thoroughly examined.These could come from a variety of sources, including: cosmic ray hits on images, inaccurate PSF models, uncertainties in the plate solution, under/overestimated astrometric uncertainties, center-of-light-center-of-mass offsets (see Appendix A), and others.In the future, it may be useful to use a likelihood model that includes a flexible noise model so MultiMoon can explicitly infer the noise properties of the astrometric measurements.

Potential For Future Observations
Future observations are key to confirming the results we present here, and will enable us to significantly narrow the uncertainties associated with our modeling.Any additional observations will strengthen the detection of non-Keplerian effects, and may enable us to break the degeneracies associated with our modeling.To illustrate the ability of future observations to improve our orbital models, in Figure 9 we have plotted the "information gain" as a result of additional observations at any given time during the first quarter of 2023, along with the separation of the binary.We define information gain as the size of the cloud of projected positions of the secondary, in units of typical HST observations.This can be written: where σ x,pred , σ y,pred are the standard deviations of the projected positions from a sample of posterior draws and σ x,typ , σ y,typ are typical error measurements of Pabu's position, based on past HST observations of Borasisi-Pabu.To calculate this, we have taken 1000 random samples from our MCMC chain and used them to predict the position of Pabu at a given time.We then use the ensemble to calculate the information gain.This is equivalent to calculating the average expected χ 2 of a new relative astrometric observation given the current posterior distribution (as listed in Table 2).For an observation with no additional constraining power, information gained is √ 2. Observations can be optimized to best shrink uncertainties by scheduling when information gain is at its maximum.In the case of Borasisi-Pabu, observations in 2023 can significantly shrink measurement uncertainties, which is unsurprising given that Borasisi-Pabu was last observed by HST in 2008.Targeting observations at times of high information gain allow for optimization of observing schedules.
In addition, observing when Borasisi and Pabu have the greatest separation ensures high data quality, further motivating observations at certain times.

CONCLUSIONS
Using a new non-Keplerian Bayesian orbit fitter, MultiMoon, we have completed a set of Keplerian and non-Keplerian orbit fits to 45 TNBs.Our Keplerian orbit fits are in close agreement with those previously completed in the literature, validating MultiMoon's orbit fitting procedures.Our exploratory non-Keplerian fits were run with the goal of identifying targets for full non-Keplerian analysis and possible observational campaigns.Almost 20% (8/45) of TNBs analyzed have improvements in fit quality when using a non-Keplerian orbit model, with many systems with moderate improvement in fit quality.The 8 systems we identify for future analysis and observation are Salacia-Actaea, Orcus-Vanth, Gonggong-Xiangliu, Borasisi-Pabu, Altjira, Varda-Ilmarë, 1999 RT 214 , and Sila-Nunam.
Our analysis is consistent with previous expectations that large TNOs are nearly spherical.The large TNB systems with detectable non-Keplerian effects are all consistent with an oblate spheroidal shape.We find that small TNBs host extremely aspherical components.Of the small TNBs with detectable non-Keplerian effects, almost all are consistent with contact binary shapes, though we are biased to highly non-spherical objects.We also identify a possible population of widely separated hierarchical triples, although these detections individually do not quite reach our improvement threshold.Our ensemble of exploratory non-Keplerian fits is consistent with their Keplerian counterparts for the TNB system's mass, semi-major axis, eccentricity, and inclination.This indicates that there are no systematics introduced by non-Keplerian effects that may invalidate past Keplerian orbital analyses of TNBs.
As an illustration of how our target list can be used to inform future non-Keplerian analyses, we have completed a full non-Keplerian fit to the mutual orbit of Borasisi-Pabu.Our analysis confirms the finding of our exploratory fits, finding that non-Keplerian motion is detectable in the system.We find that, assuming a triaxial shape, Borasisi (and/or Pabu) is extremely flattened and elongated, and may be a contact binary.Another interpretation of our results is that Borasisi (and/or Pabu) may themselves be a compact binary, making the system a hierarchical triple system like Lempo.Alternatively, our results may be caused by unidentified systematic errors present in a small fraction of our dataset.In any case, future observations are necessary to fully understand the system's architecture.6).Information gained clearly shows large variations; observing at times of maximum information gained allows for maximum leverage over reducing uncertainties.Additionally, avoiding times of close separation can ensure high quality data.
Future observations of potentially non-Keplerian TNBs are of high priority since these effects are a unique probe into TNB shapes and sizes.Many of the systems we have identified as having currently detectable non-Keplerian effects have not been observed in more than a decade allowing significant orbital precession to accumulate.High resolution relative astrometry combined with non-Keplerian models can further confirm our results, narrow uncertainties, and break degeneracies.
undetected satellite, but given a series of observations, the misalignment from a close-in satellite on a circular orbit will generally average out to zero.For an eccentric satellite however, the satellite will spend more time near apoapse, naturally causing an average photocenter-barycenter misalignment.Incomplete sampling of the putative undetected satellite's orbit and an eccentric orbit can both cause an average change in the primary's photocenter.
Third, objects with complex shapes and potentially heterogeneous interiors can have center-of-mass-center-of-figure offsets.Even neglecting the above effects, this would lead to a center-of-light offset.
A full model of potential center-of-light/center-of-mass/center-of-figure variations for non-spherical shapes under varying illumination and observation conditions would be highly complex.Most of this complexity would be well below the detectability threshold.However, to explore these effects we add two parameters to our Keplerian model that represent a simple constant offset to the modeled astrometry (in ecliptic longitude and latitude).These fits allow us to consider any potential astrometric offsets produced by aspherical objects (e.g., the center of reflected light from an object may be offset from the center of mass) or contamination from previously undetected components of the TNB system (e.g., small moons close to the primary that are unresolvable with current telescopes).These fits were run and analyzed identically to the normal Keplerian fits.
In our ensemble of Keplerian fits with constant offsets, results were varied.In most cases, the fits, when compared to the normal Keplerian fits, showed little to no improvement in quality.Despite this, a few objects had large improvements to the quality of the fits.However, the majority of these large improvements came at the cost of unrealistically large constant offsets (no significant improvements were obtained with offsets less than 5 mas).While this may indicate the presence of undetected components in the TNB system, many of the offsets were probably much too large to be caused by plausible undetected system components (offsets could be as large as ∼90 mas).
The improvement in fit is most likely due to overfitting.Our parametrization will always enable an improvement in fit by shifting the 2-dimensional residual cloud towards the origin.This is amplified by the small dataset used, with some TNBs only having 5 observations.This, combined with the physically unrealistic offsets found, convince us that the use of constant photocenter offsets is not warranted at this time.
Although we discard these fits, some insights can still be learned from examining their results (which are available upon request).In cases where there is improvement to the quality of the fit over the Keplerian (small or large improvements), the inclusion of constant offsets generally alters the fitted eccentricity of the TNB's orbit.In most cases, the fits with offsets show a small decrease in the fitted eccentricity, although increases can occur.When the Grundy best fits (without offsets) are compared to the posterior distributions from the constant offset fits (in a similar way to the data in Figure 3), the z -scores are similarly small for all parameters except eccentricity.This indicates that even if constant photocenter offsets are present in TNB relative astrometry, results based on the ensemble of TNB orbits (e.g., the mutual inclination distribution) are unaffected.
The non-detection of any believable constant photocenter offset, is not overly surprising.In reality, photocenter offsets (from any source) are time-varying.For offsets due to viewing geometry, most of the variation is caused by the Earth's orbit around the Sun, with the TNB's heliocentric orbit making a small contribution.Offsets from albedo variations and aspherical shapes will vary with the rotational period of the object in question.Future investigations of photocenter offsets should use a more physically realistic model, which accounts for the phase angle and rotation phase/rate of the modeled objects.These models may be able to constrain the shapes and albedo variations of TNB components once sufficiently high resolution data is available.Note-All fitted angles are relative to the J2000 ecliptic plane.Residual indicates the combined RMS residual for the best fit.In our fits, we use ln J 2 R 2 as our fitted parameter.Here, we present only the J 2 to enable comparison with the literature.For R, we use the volumetric radius of the modeled body (see Table 1).Bold-face data indicates TNBs that have detections of non-Keplerian motion.

Figure 1 .Figure 2 .Figure 3 .
Figure 1.Corner plot for the Borasisi-Pabu Keplerian orbit fit.Along the tops of the columns are the marginal posterior distribution for each parameter; in the case of a Gaussian distribution, this can be used to determine the best-fit and uncertainty.The contour plots show the joint posterior distribution between every pair of parameters.For example, the system-mass semimajor axis degeneracy is seen in the top left which is a result of measurement of only the orbital period and Newton's Version of Kepler's Third Law.Contours show the 1, 2, and 3σ levels of each joint posterior distribution.The horizontal and vertical lines show the orbit solution publicly available in the Grundy database.Of particular interest is the strong agreement between the Grundy solution (horizontal and vertical lines) and the solution derived by MultiMoon.The best fit in the MultiMoon fit corresponded to a χ 2 of ∼21.RMS residuals are 2 milliarcseconds in both longitude and latitude.All angles are relative to the J2000 ecliptic plane on JD 2451900.0 (2000 December 21 12:00 UT).

Figure 4 .
Figure 4. Improvements in fit quality by adopting a non-Keplerian model.The best fit J2 presented here is found by assuming a triaxial shape and dividing our measured J2R 2 by the resulting body's volumetric radius squared.The ratio L KL N K compares the Keplerian models presented in Section 3 to the best exploratory fit presented in this section.If the likelihood ratio is ≪ 1, the null hypothesis (no detectable non-Keplerian effects) could be rejected, and non-Keplerian effects are indeed detected.In general, throughout this paper, we adopt a threshold of L K L N K < 0.1 to indicate systems which have orbit fit improvements indicative of non-Keplerian motion.We have labeled all of our systems with improvements below our threshold, as well as several systems with marginal improvements.Asterisks in labels indicate either a prograde (*) or retrograde (**) orbital solution for mirror ambiguous TNBs.Color corresponds to the radius of the TNB's primary.We find that many TNBs are far better explained by a non-Keplerian model, with many systems showing moderate improvements.We also find a wide range of acceptable J2 values, ranging from oblate spheroids to unresolved binaries.Future work and observations should focus on confirming the significance of our detections and extending the observational baselines of systems with moderate improvements in fit quality.We urge caution in using the values of J2 in this plot as there may be significant uncertainties that our exploratory fits are unable to fully constrain/measure.

Figure 5 .
Figure 5. Kernel density estimates in the style of Figure 3.In this figure, we compare the non-Keplerian best fits found in our exploratory fits to the Keplerian posteriors presented in Section 3. We find that significant deviations are present in the argument of periapsis (ω) and longitude of ascending node (Ω) when assuming a non-Keplerian model.This is expected since this model allows for precession of these two angles.The distributions are still centered at zero presumably because we are comparing the non-Keplerian angles at an epoch typically near the middle of the observations and because many systems show insignificant non-Keplerian effects.Notably, however, the system mass, semi-major axis, eccentricity, and inclination are still consistent with Keplerian fits, implying that systematic errors due to non-Keplerian effects are small in most previous analyses.

Figure 6 .
Figure 6.Corner plot for the Borasisi-Pabu non-Keplerian orbit fit in the style of Figure 1.Horizontal and vertical lines correspond to the best fit solution in the Grundy database.The best fit in the MultiMoon non-Keplerian fit corresponded to a χ 2 of ∼11.5 (compared to 21 for the Keplerian fit).Instead of individually plotting the primary and secondary masses, which are completely degenerate, the system mass (Msys) is plotted to better show parameter correlations.Additionally, to allow for comparison with the literature, we show J2 and C22 rather than J2R 2 and C22R 2 , where we take the shape to be a triaxial ellipsoid with R as the volumetric radius.Of particular interest is the marginal distribution for J2 that clearly shows a high value of J2 is preferred, with J2 = 0 (equivalent to Keplerian motion) being excluded with a significance > 3σ.
Figure 7. Residual plot for the Borasisi-Pabu non-Keplerian orbit fit in the style of Figure 2. The best fit in the MultiMoon non-Keplerian fit corresponded to a χ 2 of ∼11.5.RMS residuals are 2 milliarcseconds in both longitude and latitude.Error contours here are the same size as in Figure 2.

Figure 9 .
Figure9.A plot showing the change in "information gained" by a single HST observation (in gray) and separation of Borasisi and Pabu (in red) in the first quarter of 2023.The information gained by single observations is a statistic which compares the uncertainty in the projected position of a system's secondary to the uncertainty of a typical observation (Equation6).Information gained clearly shows large variations; observing at times of maximum information gained allows for maximum leverage over reducing uncertainties.Additionally, avoiding times of close separation can ensure high quality data.

Table 1 .
System Properties of TNBs