Beyond Point Masses. III. Detecting Haumea’s Nonspherical Gravitational Field

The dwarf planet Haumea is one of the most compelling trans-Neptunian objects to study, hosting two small, dynamically interacting satellites, a family of nearby spectrally unique objects, and a ring system. Haumea itself is extremely oblate due to its 3.9 hr rotation period. Understanding the orbits of Haumea’s satellites, named Hi’iaka and Namaka, requires detailed modeling of both satellite–satellite gravitational interactions and satellite interactions with Haumea’s nonspherical gravitational field (parameterized here as J 2). Understanding both of these effects allows for a detailed probe of the satellites’ masses and Haumea’s J 2 and spin pole. Measuring Haumea’s J 2 provides information about Haumea’s interior, possibly determining the extent of past differentation. In an effort to understand the Haumea system, we have performed detailed non-Keplerian orbit fitting of Haumea’s satellites using a decade of new, ultra-precise observations. Our fits detect Haumea’s J 2 and spin pole at ≳2.5σ confidence. Degeneracies present in the dynamics prevent us from precisely measuring Haumea’s J 2 with the current data, but future observations should enable a precise measurement. Our dynamically determined spin pole shows excellent agreement with past results, illustrating the strength of non-Keplerian orbit fitting. We also explore the spin–orbit dynamics of Haumea and its satellites, showing that axial precession of Hi’iaka may be detectable over decadal timescales. Finally, we present an ephemeris of the Haumea system over the coming decade, enabling high-quality observations of Haumea and its satellites for years to come.


INTRODUCTION
Almost all of the largest transneptunian objects (TNOs) are known to host satellites (e.g.Christy & Harrington 1978;Brown et al. 2005Brown et al. , 2006;;Brown & Suer 2007;Noll et al. 2007;Parker et al. 2016;Kiss et al. 2017).These satellites are generally small relative to the system primary, and are thought to have formed by collisions (Barr & Schwamb 2016;Arakawa et al. 2019).The current density of the transneptunian region is far too low to have formed so many satellites by collision (Campo Bagatin & Benavidez 2012;Abedin et al. 2022), implying that these systems must not have formed insitu.The emerging consensus is that large TNOs formed in a relatively massive primordial disk exterior to the giant planets, which was subsequently scattered by Nep-Corresponding author: Benjamin Proudfoot benp175@gmail.comtune's outwards migration (Nesvorný 2018;Gladman & Volk 2021).The large TNOs we see today, which are on excited orbits, are the remnants of this primordial disk.By understanding how the satellites of large TNOs formed and evolve, we can probe the conditions of the early primordial disk where these systems formed.
(136108) Haumea, the third most massive TNO known, is host to two satellites: Hi'iaka on a ∼50 day orbit and Namaka on a ∼20 day orbit (Brown et al. 2005(Brown et al. , 2006).Haumea's shape, determined by both light curve observations (Rabinowitz et al. 2006) and stellar occultations (Ortiz et al. 2017), is significantly nonspherical due to its 3.9 hour rotation period.Haumea and its satellites may have formed during a collision, which simultaneously spun up Haumea, created the satellites, and also formed Haumea's unique, icy collisional family (Leinhardt et al. 2010;Proudfoot & Ragozzine 2022), but there remains some disagreement on these circumstances (e.g., Ortiz et al. 2012;Campo Bagatin et al. 2016;Noviello et al. 2022).Connecting a formation model to all of the system's unique characteristics has been difficult despite many proposals (e.g., Proudfoot & Ragozzine 2019).
Our understanding of the complex nature of the Haumea system can be advanced by detailed study of the satellites' orbits.Study of the orbits potentially allows for a measurement of the masses of each component, as the two satellites strongly interact with each other (Ragozzine & Brown 2009, hereafter RB09).In addition to satellite-satellite interactions, perturbations from Haumea's nonspherical gravitational potential are also present.Haumea's gravitational potential is determined both by its shape and internal density distribution.Since Haumea's shape is fairly well known (due to observations of a stellar occultation, see Ortiz et al. 2017), measuring the gravitational harmonics of Haumea may allow us to constrain its internal density distribution.The internal density distributions of TNOs are almost completely unconstrained, although large TNOs are expected to be differentiated (McKinnon et al. 2008;Dunham et al. 2019).
Haumea, in particular, has evidence for differentiation in the form of the collisional family.Haumea's family members are spectroscopically unique showing strong water ice features and few other major constituents (Souza Feliciano et al., in prep.).In combination with very high albedos (Elliot et al. 2010;Vilenius et al. 2014), it seems like the family members could be primarily composed of nearly pure water ice.This may suggest that the proto-Haumea was a differentiated "ocean world" and that Haumea family members are pieces of the water ice mantle.Thus studying presentday Haumea could give us unique insight into the interiors of ocean worlds.
Haumea's gravitational potential can be described by a spherical harmonic expansion.To quadrupole order, the gravitational field of Haumea, U , with mass, M , at a distance r can be written as: (1) where J 2 and C 22 are the second-order gravitational harmonic coefficients, θ is the body-fixed latitude angle, ϕ is the body-fixed longitude angle, and R is an "reference" radius (Yoder 1995;Scheeres et al. 2000).In this work, we assume that R is equivalent to volumetric radius.The coefficients J 2 and C 22 describe Haumea's shape and internal density structure.Taking the shape found by Ortiz et al. (2017) and assuming a homogeneous density, Haumea is expected to have J 2 = 0.24 and C 22 = 0.05.However, when taking differentiation into account and using the Dunham et al. (2019) model of Haumea's interior, these harmonics would be J 2 = 0.16 and C 22 = 0.03.While both of these models are simplified, they serve as a useful guide.
In the original work that determined the orbits of Haumea's satellites (RB09), Haumea's nonspherical gravitational field was not clearly detected, although they were able to robustly detect satellite-satellite interactions.Subsequent follow up studies have also been unsuccessful in detecting the nonspherical field (Gourgeot et al. 2016).However, with new ultra-precise HST observations from the past decade, another analysis of the satellites' orbits is in order.By leveraging these observations, as well as new computational techniques, we present a new, updated set of orbital fits to the Haumea system.We are able to detect the nonspherical gravitational potential of Haumea, constrain the masses of Haumea's satellites, and study the spin-orbit evolution of the system.
The paper will be formatted as follows.In Section 2, we describe the observations used in our analysis.Then in Section 3, we describe our non-Keplerian orbital model and fitting procedure.Results of the fitting are presented in Section 4, and discussed in Section 5. We then conclude in Section 6 and discuss future work.

OBSERVATIONS AND DATA ANALYSIS
The data we use in our orbit fitting comes from a variety of sources, but can be broadly broken into three separate datasets.The first dataset comes directly from RB09, which extracted satellite positions from Keck and HST observations.The second dataset consists of HST observations from HST Programs 12243 and 13873.The last dataset is made up of Keck observations from 2020-2022.For our orbit fitting, we combined the relative astrometry from each data set and simultaneously fit all the data.Our compiled data is presented in Table 1.
The published astrometry from RB09 was found to have a sign error in their listed RA offsets (in their Table 1).This error can be seen in their Figure 1 as RA decreases towards the east, opposite to convention.This mistake affected their orbital modeling, preventing them from correctly determining the orbital plane of the system, although the rest of their analysis is relatively unaffected.In our analysis, we use the RB09 data, although we correct the error and use the mirrored RA values.
HST Programs 12243 and 13873 used HST's Wide Field Camera 3 (WFC3) to observe the Haumea system with a combined 13 orbits of coverage.Program 12243 imaged the system, over 10 consecutive HST orbits, in an attempt to observe a Haumea-Namaka mutual event.Program 13873 used 3 single-orbit visits to measure satellite relative astrometry to better constrain orbit models.Both of these programs took ∼30 individual exposures per orbit, with exposure times of 49 and 39 seconds for the programs respectively, using the F350LP filter to maximize signal-to-noise.The images from these programs were analyzed using the same method used in RB09, although changes were made to fit the WFC3 data, replacing the older PSF models used in previous studies.
Keck observations used the laser guide star adaptive optics system (LGS AO) (Wizinowich et al. 2006) with the narrow camera of NIRC2 1 .In the 2020 and 2021 observations, nearby field stars were used for tip-tilt correction, since they were brighter than Haumea.All Keck observations were done in the infrared H filter, covering wavelengths from ∼1.48 to 1.77 µm, with a series of dithered exposures for sky subtraction and to minimize the effect of bad pixels.120 second exposures were taken during the 2020 and 2022 observing runs, while 60 second exposures were used in 2021.Unfortunately, due to the short exposure time, Namaka was not visible in the 2021 images.During the 2022 observing run, unfavorable orbital phase also prevented detection of Namaka.Pairs of dithered images were flat-fielded and pairwise subtracted to remove sky background using practices common in TNO binary observations (e.g., Grundy et al. 2011).
To extract relative astrometry of the satellites from the processed Keck data, we simultaneously fit 2dimensional Gaussian PSFs to each visible object in individual processed images.While a Gaussian is a relatively poor approximation for the NIRC2 PSF, it is still able to measure the center of each PSF quite accurately.Relative detector positions were then converted to relative right ascension and declination assuming a mean plate scale of 9.952 mas/pixel and an orientation offset of 0.252 • (Konopacky et al. 2010;Yelda et al. 2010;Service et al. 2016).The median and standard deviation offsets of individual measurements are used for the astrometric offsets and error for each night, although we implemented a conservative noise floor of 10 milliarcseconds to account for unknown systematics.This method has been extensively used to extract relative astrometry from Keck NIRC2 images (e.g.Fraser & Brown 2010;Grundy et al. 2015).
1 https://www2.keck.hawaii.edu/inst/nirc2 As can be seen in Table 1, both satellites are not always detected at each epoch.In principle, nondetections can be used to help constrain the satellites' orbits, but in practice, given the already well-known orbits, they barely constrain the fits.Hence, during our orbit fitting process, we do not use non-detections in any way.

METHODS
For our orbit fitting, we use MultiMoon, a state-of-theart orbit fitter designed for use with TNOs (Ragozzine et al., submitted).MultiMoon is built around an nquadrupole integrator that can simulate the gravitational and rotational dynamics of an arbitrary number of tri-axial ellipsoids to quarupole order.Internally, it uses emcee (Foreman-Mackey et al. 2013, 2019), a popular Markov chain Monte Carlo (MCMC) ensemble sampler, allowing us to treat the orbit fit as a Bayesian inference problem.In its simplest form, MultiMoon uses a least squares method for evaluating the goodness-of-fit of a given orbital model.It can also accommodate a more complicated goodness-of-fit metric which we describe later in this paper.
In our fits, we only consider the J 2 since it is by far the most dominant harmonic.The other second order harmonic C 22 -which is related to the prolateness of Haumea-is relevant for understanding dynamics of orbits around Haumea, but only within a few times the corotation radius (Proudfoot & Ragozzine 2021).Even within this range, C 22 averages out over many orbits except when close to a spin-orbit resonance (SOR).Haumea also certainly has substantial higher-order harmonics (most notably J 4 ), but their effect is small due to their r −5 distance dependence.As further justification of this assumption, we can analytically estimate the precession induced by Haumea's J 4 , and find it is only ∼0.1% the strength of J 2 precession for Namaka, and even smaller for Hi'iaka.Thus we believe that our simple model of Haumea's gravitational potential is sufficient to describe the dynamics taking place in the Haumea system.
For the orbit fits presented here, we only model the gravitational harmonics of Haumea ignoring the (presumably) nonspherical shapes of Hi'iaka and Namaka.We revisit this assumption later in the paper.Our baseline orbit model has 18 free parameters including the mass, J 2 , and 2 spin pole direction angles of Haumea, in addition to the masses and 6 orbital elements of each satellite.Our model also requires the input of Haumea's rotational period to correctly model any axial preces-    Note-The relative right ascension and declination positions of Haumea's satellites, Hi'iaka (H) and Namaka (N).At some epochs, Hi'iaka or Namaka were not visible in the images, for a variety of reasons.For these entries, no data is listed and our orbit fits were not constrained by their non-detection.Data from before 2010 are taken from Ragozzine & Brown (2009), although we correct their sign error in the ∆x columns.
sion which the satellites may cause.Although this value could, in principle, be a free parameter in the model, it is known with high precision and has very little influence on the orbital dynamics of the system.Hence we opt to use a fixed value of 3.915 hours (Rabinowitz et al. 2006).
To account for possible systematics arising from the use of a variety of data sets, we have implemented a sophisticated likelihood function within MultiMoon.This likelihood function is adapted from the outlier pruning methods presented in Hogg et al. (2010).Since we, a priori, do not describe the systematic errors that may arise in the fitting process, we use an extremely flexible framework.Our likelihood model is a mixture model that combines two least-squares terms.The first is a common least-squares likelihood model, the standard technique for orbit fitting.This term is combined with another least-squares model with an additional error term.The error term, which we call σ sys , is combined with the measured uncertainties of our observations in quadrature.Also included is a normalization factor (f sys ) describing the fraction of data displaying systematic errors, which also acts as a penalty for exclusion of data.The entire likelihood function can be written: where y i and σ i are the N observations and uncertainties, and y i,m is the model.Technically, y i,m and σ i are vectors where each have two dimensions (∆α cos δ, ∆δ) and there is an implied summation over both of these dimensions.For brevity, however, we exclude this summation, although it is implemented in its full form internally.If there are significant outliers in the data, this prescription downweights them relative to the typical least-squares assumptions and thus qualifies as a "robust" (to outliers) statistical method.In this sense, it operates similar to an automated sigma-clipping technique.It also allows for the expansion of systematic uncertainties when the quoted statistical uncertainties are too small to explain the scatter in residuals relative to the model.The factors (1 − f sys ) and f sys are critical for normalizing the two likelihood models and provide an implicit prior that penalizes overestimation of systematic effects.However, when σ sys ≪ σ i (i.e.there are no systematic errors present in the data), f sys is not well defined as both likelihood functions asymptotically approach one another.To prevent this degeneracy from becoming problematic, we implement a prior forcing σ sys ≥ 1 milliarcsecond.This likelihood model adds an additional two free parameters to our model (σ sys , f sys ).However, instead of fitting σ sys , we opt to fit log 10 σ sys , allowing the MCMC algorithm to more easily explore a greater range of values.This "robust" likelihood model has now been implemented into MultiMoon and is publicly available on GitHub2 .We have extensively validated this likelihood model using synthetically produced data sets that have large systematics applied.We find that when using this model, MultiMoon can recover the original orbital parameters even when systematic uncertainties of 10s of milliarcseconds are applied to ∼50% of the data set.
During our data fitting process, we found that large systematics were present when combining both the Keck and HST data that necessitated the use of this robust likelihood model.Unfortunately, our model could not resolve these issues and unusual systematics remained unaccounted for.To remain as conservative as possible, we elected to complete an orbit fit using the HST data only with a standard least-squares likelihood model.We discuss the drawbacks of the HST+Keck fit further in Section 4 and 5.
As part of the Bayesian framework MultiMoon uses, we set priors for all parameters to be uninformative (except for σ sys as discussed above), allowing the data to constrain the posterior distribution.However, in our HST only fit, we set uniform priors on the spin pole direction of Haumea to prevent walkers from getting stuck in a lower dimensional subspace.The priors were chosen to bracket the best region of likelihood space within ∼ 10 • , as identified in preliminary fits.After the fit was completed, we confirmed that this prior did not significantly prevent walkers from exploring favorable parts of likelihood space.
We drew initial walker positions from Gaussian distributions centered near the location of highest likelihood that was identified in preliminary runs.These preliminary runs were conducted to broadly search parameter space and used very broad initial guesses, allowing for a rigorous search of the 18-dimensional parameter space (20-dimensional for the HST+Keck fits).Our preliminary fits showed no signs of other likelihood maxima with acceptable fit quality.Our baseline orbit fits used 960 walkers in the MCMC ensemble, which were run for 5000 burn-in steps.We then pruned underperforming walkers, replacing them with random linear combinations of highly performing walkers, after which the Figure 1.A corner plot for the HST only orbit fit to the Haumea system.We include all 18 fitted parameters along with 2 derived parameters.To facilitate easy interpretation, we list J2 rather than J2R 2 by taking the volumetric radius R from the occultation derived shape model (Ortiz et al. 2017).We also show the inclinations of each satellite with respect to Haumea's equator in the last two columns.Along the top of each column is the marginal posterior distribution for each parameter in our fit.Below the marginal distributions are the 2-dimensional joint posterior distributions for every pair of parameters.Contours correspond to 1, 2, and 3σ levels.Small black points mark individual samples from our MCMC chains.The best fit parameter set in our MCMC chains corresponds to a χ 2 of 99.1 with 90 degrees of freedom.Of particular note is the strong exclusion of J2 = 0 in the marginal posterior for Haumea's J2, alongside strong correlations between J2 and a variety of other parameters.
ensemble was run for 1000 more burn-in steps.The ensemble was then run for 20000 steps to sample the posterior distribution.We confirmed that the resulting chains were converged by visual inspection of walker trace plots and marginal parameter-likelihood plots.

RESULTS
When comparing our different orbit fits, we find that there is strong disagreement between the two datasets (HST+Keck and HST only).When fitting to the combined dataset, we find that the most recent Keck observations of the system are at odds with the 2014-2015 HST observations.Using our robust likelihood model and the combined HST+Keck dataset, we find that our best fit is ∼10σ inconsistent with the 2014 HST observation of Namaka.This inconsistency was attributable to the data rather than the model, as shown by fits both with and without our robust likelihood model.Our robust likelihood model parameters indicated that approximately 10% of the data had uncertainties underestimated by ∼10 milliarcseconds.Fits without the robust likelihood model were extremely similar to those with it, except the fit quality was much worse.When closely examined, no obvious problems were seen in the Keck or HST images, and no difficulties arose during our analysis of the images.To examine whether our image analysis techniques were to blame, we attempted to extract astrometry from the images with a variety of techniques (e.g., Gaussian PSF fitting, WFC3 model PSF fitting, etc.), all of which yielded similar results.
In addition to the internal inconsistency, the HST-Keck combined fit also produced a measurement for Haumea's J 2 that was too high to be compatible with other observations of the Haumea system (see Section 5 for more discussion).In comparison, the HST only fit showed no such issues.When the orbit fits are compared, very little changes between the models with the exception of Hi'iaka's mass, Haumea's J 2 , and Haumea's spin pole direction.As unknown systematic errors are affecting our orbit fit, we choose to proceed by eliminating possible sources of these systematic errors.Since HST's PSF is extremely stable and has been extensively cross-calibrated across instruments, we adopt the HST only orbit fits for the purpose of this work.This choice results in larger uncertainties within the model, but allows us to be more confident that our results are not affected by systematics.Although we adopt the HST only fit, we still discuss the implications of our combined orbit fit in Section 5, as well as reporting its results in Table 2. x and y correspond to ecliptic longitude and latitude, which is the primary coordinate system used in MultiMoon.Hi'iaka residuals are shown as circles and Namaka with triangles.The color of each point corresponds to the observation date.The circles correspond to 1, 2, and 3σ error contours.As reported in the text, this fit corresponded to a reduced χ 2 ∼ 1.10.We find that although the fit quality is worse than would be desired, the p-value of associated with the χ 2 statistic is 0.23, indicating an acceptable fit.
The results presented here are our most refined orbital fits.Including our preliminary analysis, exploratory fits, and fits using different likelihood models, our nominal orbit model is the result of well over 10 9 individual orbit integrations.We show the HST only orbit model in it entirety as a corner plot (Foreman-Mackey et al. 2016) in Figure 1.Each column in the corner plot displays the marginal posterior distribution for each parameter as a histogram (at the top) and 2-dimensional joint posterior distribution as a contour plot.In addition, we display the marginal posteriors for both fits in Table 2.Both Figure 1 and Table 2 also contain several derived parameters, parameters that are functions of our fitted parameters.To display the fit quality, we show the residuals of the best fit parameter set in Figure 2, alongside 1, 2, and 3 σ error contours.This best fit parameter set is only one realization of our posterior distribution, but it illustrates the quality of our fit.
One of the outstanding features of our orbit fit is our detection of Haumea's J 2 .When assuming the volumetric radius derived from stellar occultation measurements (798 km, Ortiz et al. 2017), we find J 2 = 0.262.However, our orbit fit shows that Haumea's J 2 and Hi'iaka's mass are highly degenerate with one another.In Figure 3, we show, in detail, the degeneracy between these parameters as a function of reduction in fit qual-  Note-Reported values represent the median value taken from the posterior distribution, while the stated uncertainties represent the 16th and 84th percentiles.All fitted angles are relative to the J2000 ecliptic plane on Haumea-centric JD 2454615.0 (2008 May 28 12:00 UT), chosen to match the epoch used in Ragozzine & Brown (2009).Assumed c-axis for Haumea is 537 km (Dunham et al. 2019) and spin period is 3.915 hours (Rabinowitz et al. 2006), however, altering these values produces no meaningful change to the fit.
To transform to J 2 from only the more physically meaningful J 2 R 2 , we use a volumetric radius of 798 km (Ortiz et al. 2017).
ity.It is clear that a large range of values for these two parameters are acceptable, with nearly no reduction in fit quality.In our HST+Keck orbit fit, we find that Haumea's J 2 has much lower uncertainties, but is unexpectedly high, J 2 = 0.431.Although probably attributable to unidentified systematic errors in our dataset, we will discuss possible causes/interpretations of this unusual measurement in Section 5. Our detection of Haumea's J 2 is significant in both orbit fits.The HST only fit detects J 2 at ∼2.5σ confidence, while the HST+Keck fit detects it at > 5σ confidence.Alongside our detection and measurement of Haumea's J 2 , we also provide a measurement of Haumea's rotation pole.We find that Haumea's pole (or more precisely, the pole of Haumea's gravitational quadrupole) points toward (α p , δ p ) = (282.9+0.6 −0.7 • , −9.7 +0.6 −1.0 • ), very close to the occultation derived rotation pole of (α p , δ p ) = (285.1 ± 0.5 • , −10.6 ± 1.2 • ) (Ortiz et al. 2017).
In our orbit fit, we are able to significantly detect the masses of both satellites at significance > 3σ.RB09 previously detected Namaka and Hi'iaka's masses, but only with 1.2σ confidence for Namaka.While our fit strongly detects both, the uncertainty on Hi'iaka's mass is substantial due to its degeneracy with Haumea's J 2 .Along- side mass measurements, we are also able to constrain the satellites' inclinations with respect to Haumea's equator.We find inclinations of 12.8 +0.8 −0.6 • and 1.0 +0.6 −0.5 • , for Namaka and Hi'iaka, respectively.We also measure the satellites' mutual inclination of 13.2 +0.2 −0.2 • .Our orbit fits are significantly different from past orbit fits (Ragozzine & Brown 2009;Gourgeot et al. 2016).While this difference is expected since we include more dynamical effects (e.g.including J 2 ), some important differences are still present.Most notable is the change in orbit angles, which stems from RB09's incorrectly tabulated astrometry, allowing for close agreement with the orbit planes found in Gourgeot et al. (2016).We find a lower eccentricity for Namaka (0.2179 +0.0032 −0.0033 ) compared to RB09 (0.249 ± 0.015 using the same epoch), also presumably due to their incorrect astrometry.Another notable difference is the change in Hi'iaka's mass (12.13 +3.22  −3.11 × 10 18 kg) when compared to RB09 (17.9 ± 1.1 × 10 18 kg), due to our inclusion of Haumea's J 2 .Our preliminary fits showed that our orbit model, when eval-uated with a small J 2 approximately reproduces RB09's measurement of Hi'iaka's mass.
When compared to the orbit model presented in Gourgeot et al. (2016), we find quite large differences in orbital parameters especially in the fit for Namaka's orbit.This is unsurprising since their orbit model was a pure Keplerian orbit fit, neglecting both Haumea's J 2 and satellite-satellite interactions.Their analysis claimed that there was no signature of non-Keplerian effects caused by Haumea's J 2 in the system, although they use a much shorter span of data than our analysis.We find that non-Keplerian effects from both satellitesatellite interactions and Haumea's J 2 are strongly detected, however it remains uncertain how strong each effect is.Assuming a homogeneous density structure and using the equations found in Yoder (1995), the occultation derived shape implies a J 2 of 0.24 (Ortiz et al. 2017).Allowing for differentiation decreases J 2 significantly.The model for a two-layer differentiated Haumea presented in Dunham et al. (2019), gives an overall J 2 of ∼0.16.Our model fitting to all available data (HST+Keck) is 3σ inconsistent with both of these models.However, the fit with only HST data is consistent with both, encouraging us to explore possible reasons Haumea's J 2 may be higher.
One possible reason could be the gravitational contributions from Haumea's ring.Assuming a circular ring, the following expression can be derived for the J 2 contribution of a ring: where M r and M P are the masses of the ring and Haumea, respectively, and r is the radius of the ring.
For the ring to contribute ∼1% of the measured J 2 R 2 of our HST only fit, the ring would need to have a mass of ∼ 10 18 kg, about the mass of Namaka, given the known ring radius of 2287 km.For it to be the cause of Haumea's unexpectedly high J 2 in the HST+Keck fit, the ring would need to be two orders of magnitude more massive, equivalent to tens of Hi'iaka masses.While no mass constraints on the ring are found in the literature, this value is absurdly high.For comparison, Saturn's rings are ∼1 Hi'iaka mass (Iess et al. 2019).Hence, the ring is unlikely to contribute significantly to our measured J 2 .There is a distinct possibility that more rings may be detected in future occultations (e.g.Pereira et al. 2023), but even when combined, a ring system is unlikely to contain enough mass to substantially contribute to Haumea's J 2 .Another potential source is an undetected satellite within Namaka's orbit.Averaged over an orbit, the putative inner satellite would act similar to a solid ring of material.Hence, using Equation 3 above to calculate the J 2 of a putative inner satelite orbiting at 10000 km, we find that an inner satellite with a mass near Namaka's mass could significantly contribute to Haumea's measured J 2 .Unfortunately, Burkhart et al. (2016) significantly ruled out satellites more than 60 km in diameter closer than 10000 km by using a nonlinear shift-andstack image analysis technique.This diameter implies a mass approximately one-tenth of Namaka's mass, which would scarcely contribute to the overall J 2 .Given this, we believe it is unlikely that our results could be caused by an unknown inner satellite.Likewise, undiscovered satellites external to Hi'iaka's orbit would not produce the observed signature.
An alternative reason could be an extreme mass anomaly on Haumea's surface, either positive or negative, which would cause Haumea's J 2 to be larger than expected.Unfortunately, the mass surplus (or deficit) would have to be substantial, of order ∼10% of Haumea's total mass.Maintaining a mass anomaly of that magnitude would require Haumea to have implausibly high material strength.Using the method developed by Johnson & McGetchin (1973), we optimistically estimate that Haumea may support a maximum topographic feature of ∼10 km, amounting to far less than the ∼100 km required to produce an unusually high J 2 .Hence, it seems unlikely that a surface feature could plausibly explain the HST+Keck measurement.
Likewise, Haumea could have an unusual interior density distribution.If Haumea formed in the aftermath of a large collision, it may have an unusually shaped core left over as a remnant of this impact.Alternatively, its core could be offset relative to its external figure, potentially explaining Haumea's "Dark Red Spot" (Lacerda 2009).Assuming Haumea's core is triaxial and adopting the external figure of Haumea as measured by stellar occultation, we can calculate the J 2 of Haumea with an arbitrary triaxial core shape with an offset.We find that for any realistic core shape or offset, Haumea's J 2 can not increase by more than 50%, still well below the constraint provided by the HST+Keck fit.In any case, the extreme versions of these hypotheses are unlikely to be geophysically viable as they ignore fluid-like relaxation of Haumea's core and mantle.We believe that an unusual interior is unlikely the cause of our results.
Due to the implausibility of all of these solutions, we conclude that our result is due to factors dependent on our modeling techniques or data.One possible explanation is higher order non-Keplerian dynamics taking place within the system.Since our model only includes Haumea's J 2 , other gravitational harmonics may be needed to fully model the system.To investigate the effect of C 22 , we ran a MultiMoon fit that included Haumea's C 22 potential.This fit gave nearly identical results and found no constraint on C 22 , indicating that the orbital dynamics of the system are not strongly coupled to C 22 .Indeed, previous work exploring the effects of C 22 found that it is only relevant when P orb ∼ P spin , or near a low order SOR (Proudfoot & Ragozzine 2021).Beyond quadrupole dynamics, fourth-order, or hexadecapole, dynamics could contribute to orbital precession, but as previously discussed their r −5 dependency ensures that their contributions are small.Taking the J 4 harmonic as an example, we find that the ratio of nodal precession of Namaka's orbit caused by the J 4 and J 2 harmonics is ΩJ4 / ΩJ2 ≈ 0.003 J4 J2 .Apsidal precession has a similarly small strength.Typically, the J 4 harmonic is much smaller than J 2 , implying the nodal precession induced by J 4 is a very small effect.
Odd harmonics (e.g., J 3 , C 31 , etc.) could, in theory, also play a role in the system's orbital dynamics.(Note that the "dipole" term is 0 due to using the center-ofmass as the coordinate system; a center-of-mass-centerof-figure offset can contribute to J 2 which was included in the calculations with the offset core above).Again, finding the ratio of nodal precession for Namaka, we find ΩJ3 / ΩJ2 ≈ 0.006 J3 J2 when ΩJ3 is at its maximum.Odd harmonics are typically extremely small in bodies at (or near) hydrostatic equilibrium, producing ∼no detectable precession.While it seems unlikely that J 3 or other odd harmonics cause significant changes in the dynamics on the timescale of our observational data, future investigations should explicitly test whether odd harmonics are necessary for accurate modeling of the Haumea system.
Another possible effect that we do not account for is the satellites' putatively nonspherical gravitational fields.Our model assumes that Hi'iaka and Namaka are point masses, although they are likely to be substantially nonspherical.In some cases, however, the gravitational harmonics of a system's secondary can play a major role in the overall orbital dynamics (e.g.Ragozzine & Wolf 2009).MultiMoon is well poised to explicitly test this assumption.Rather than adding six parameters to our overall model, which would be computationally expensive, we can simply add the satellites' harmonics as fixed values.We can then compare a model without their harmonics to one with them, allowing us to see the change in system dynamics.In this comparison, we use the characteristics of all our observations (HST+Keck) to explicitly connect the dynamical change to actual observability.In Figure 4, we show the change in orbit fit quality as a function of Hi'iaka and Namaka's J 2 , when both satellites have a moderately high obliquity.We define change in fit quality as ∆χ 2 = i (x i,J2 −x i,J2=0 )/σ i,obs , where x i,J2 and x i,J2=0 are synthetic astrometric measurements from models including satellite J 2 and those without and σ i,obs are the measurement uncertainties for the system as tabulated in Table 1.Using those measurement uncertainties allows us to connect the system's dynamics to the actual observations.Overall, we find the fit quality is barely decreased when reasonable values for J 2 are tested.Hi'iaka would need J 2 > 2 for a detectable change, while Namaka would need J 2 > 10.For comparison, Arrokoth, a contact binary, has a J 2 ∼ 0.3.
While not an exhaustive search of parameter space, this is strong evidence that Hi'iaka and Namaka's nonspherical shapes do not significantly contribute to the system's dynamics with the present data.
In our view, the only remaining option is the presence of unknown systematic errors plaguing our dataset.Despite our use of novel statistical techniques, our model cannot account for all systematic errors arising from combining our dataset.For example, time-varying distortions in the NIRC2 field cannot be appropriately accounted for by our model.Likewise, wavelengthdependent offsets between Haumea's center of light and center of mass, potentially caused by the known wavelength-dependent rotational variability known as the "Dark Red Spot" (Lacerda 2009) may introduce unwanted systematics when combining the data sets.Indeed, when combining the datasets, we find that the Keck data from the 2020s is incompatible with the HST visit from 2014.Our combined dataset produces large residuals for the 2014 HST visit, while the HST only fit shows no such effect.While disconcerting, this conclusion is not extremely surprising given a similar result in RB09, from which we draw much of our data.Those authors similarly found that the Keck data was inconsistent with the HST-only fits.We thus argue that unknown systematic errors are the source of our unusually high measurement of J 2 .HST instruments are extremely well-studied and have been rigorously crossvalidated, so we view the HST only fit as more trustworthy.To remain as conservative as possible, we adopt the HST only fit as our nominal model for the rest of the analysis in this work.

Masses and Densities of Hi'iaka and Namaka
)/σi,obs, where xi,J 2 and xi,J 2 =0 are synthetic astrometric measurements from models including satellite J2 and those without.σ i,obs are the true measurement uncertainties for the system as tabulated in Table 1.This formulation allows us to directly determine whether the satellites' J2values produce detectable changes in the system relative astrometry, given our current observations (HST+Keck).The parameters for the models used were taken from the best fit of our orbit fits.The rotational poles of the satellites were chosen to have high obliquities (w.r.t.their Haumea-centric orbits) to enhance the effect of J2.The three separate lines show each satellite's individual contribution, as well as a comparison where both satellites have (the same) J2.For reasonable values, J2 ≲ 0.5, very little change in fit quality is found, although for extremely nonspherical shapes, Hi'iaka's J2 could begin to alter the fit.
While our model does not fully constrain the masses of the satellites, especially Hi'iaka, it is still scientifically valuable given some assumptions.If we assume Haumea has J 2 of a differentiated body (J 2 = 0.16 using the Dunham et al. ( 2019) model), we are able to estimate the mass of Hi'iaka to be roughly 1.6 × 10 19 , with uncertainties of ∼10%.This is similar to the mass determination found in RB09, albeit 10% smaller.Given an estimated radius of ∼160 km (RB09), this gives a density of ∼900 kg m −3 , similar to other TNOs in the same size range (Grundy et al. 2019).Even given the undifferentiated model, Hi'iaka is not substantially less massive, altering the density marginally.
While the mass measurement for Hi'iaka must be improved to place better constraints on Hi'iaka's density, this needs to be matched with improvement in the radius determination.As density is proportional to r −3 (as opposed to linear in mass), uncertainties in radius have a stronger influence on the overall density uncertainty.Thermal observations and modeling of the satellites in-dicate high albedos and small sizes (Müller et al. 2019), but uncertainties are large.At this time, given the massive uncertainties in both radius and mass, robust interpretation of the density of Hi'iaka is impossible.
Fortuitously, however, a recent occultation campaign has captured a multi-chord stellar occultation of Hi'iaka (Fernandez-Valenzeula et al., in prep).The results of this campaign are forthcoming, but early indications point towards Hi'iaka being larger than expected (E. Fernandez-Valenzeula, pers. comm.).When using the results from those observations, future work should be able to tightly constrain the density of Hi'iaka if the Hi'iaka mass-J 2 degeneracy is broken.
While our uncertainty on the mass of Namaka is much smaller than for Hi'iaka, Namaka still suffers from a poorly constrained size.However, given size estimates from thermal measurements (Müller et al. 2019) and past photometric analysis (RB09), Namaka could have a radius of ∼75 km, giving a density of ∼650 kg m −3 .While quite low, it is typical of small TNOs (Grundy et al. 2019).Such a low density is indicative of a porous interior, consistent with relatively gentle accretion in a satellite forming disk around Haumea in the aftermath of an impact.To improve confidence in this measurement, Namaka should be a high-priority target for an occultation campaign.

Haumea's Pole
Among TNOs, very few spin poles have been constrained.When disregarding non-Keplerian fitting, the only techniques currently able to characterize the spin poles of TNOs are long-term light curve monitoring and occultations.Light curve inversion techniques require observations of a TNO over a significant portion of their orbit, which is difficult due to TNOs' long heliocentric orbital periods.Occultations are extremely powerful for inferring spin poles, but observations of multiple multichord events are required, which are only available for a few TNOs.Non-Keplerian orbit fitting now adds an additional tool which can be used to understand the spin poles of TNOs 3 .Normally, non-Keplerian fitting cannot find a unique pole solution, but is able to determine the angle between the primary's spin pole and the secondary's orbit normal.However, since Haumea has two satellites, the spin pole can be found unambiguously.
Our measurement in this work represents the first dynamically determined spin pole among TNOs.We are able to place tight constraints on the spin pole direction of Haumea, finding (α p , δ p ) = (282.9+0.6 −0.7 • , −9.7 +0.6 −1.0 • ).Interestingly, Haumea's pole is almost perpendicular to its orbit.Based on Haumea's heliocentric orbit, we find Haumea's obliquity to be 87.1 • .Our dynamically determined spin pole measurement lies 2.3 • from the ring pole determined by stellar occultation (α p , δ p ) = (285.1±0.5 • , −10.6±1.2 • ) (Ortiz et al. 2017).Given that our methods are completely different from those used in analyzing the stellar occultation and we use no constraints from that work, the close match is encouraging.Formally, however, these two spin pole measurements are ∼2σ apart.The difference could be a real effect (i.e., Haumea's ring is inclined with respect to Haumea's nonspherical gravitational field), but the nodal precession of ring particles induced in this case is likely to erode the rings excessively (Marzari 2020).It seems much more likely that the ring is coplanar with Haumea's gravitational field, and the disagreement is due to model dependent factors or random chance.For example, Ortiz et al. ( 2017) modeled Haumea's ring as a flat, circular annulus, but the true ring may have substantial eccentricity and/or thickness.Alternatively, since the spin pole in our fit is highly correlated with J 2 , our large uncertainties may cause/contribute to the disagreement.When examining our HST+Keck fit, we find the spin pole measurements are even further apart, lending further credibility to the HST only model.

Rotational Dynamics
The rotational dynamics of the Haumea system are extremely interesting to study.For Haumea itself, the torque from both satellites on Haumea's nonspherical body cause a small amount of axial precession.We show a 1000 year integration of Haumea's spin dynamics in Figure 5.The integration is performed by SPINNY, the spin-orbit integrator at the heart of MultiMoon (Ragozzine et al., submitted).We find that Haumea's axis precesses by < 1 • on a timescale of 100s of years.Visible is a complex precession cycle in Haumea's spin pole direction with one 'fast' frequency and one 'slow' frequency.These two frequencies are caused by the torque from each satellite with periods corresponding to each satellite's nodal precession period.Namaka's nodal precession period is strongly coupled with both Hi'iaka-Namaka interactions and Namaka-J 2 interactions.Since Namaka is coupled to Hi'iaka's nodal precession, Namaka's nodal precession then weakly couples to Hi'iaka-J 2 nodal precession, although this is a much smaller effect.Hi'iaka's nodal precession has a fast, low-amplitude component caused by Hi'iaka-Namaka interactions as well as a slow, high-amplitude component from the .The spin precession of Haumea over a 1000 year integration.In the top two panels, we show the spin precession of Haumea in terms of pole ecliptic longitude and latitude.In the bottom two panels we show the precession of the orbit normal of Hi'iaka and Namaka, again in terms of ecliptic longitude and latitude.The integration parameters have been chosen to be representative of the posterior found in Table 2. Similar integrations with different values for Haumea's J2 change the long-term precession period, but are qualitatively similar.Haumea's spin precession is coupled with the nodal precession of the satellites.High frequency, low amplitude variations in Haumea's pole direction are caused by Namaka's rapid precession, while low frequency, high amplitude variations are coupled with Hi'iaka's J2 precession.As can be seen, the precession of all components are strongly coupled, both through Hi'iaka-Namaka gravitational interactions and interactions with Haumea's J2.
Hi'iaka-J 2 interaction.The low-amplitude, high frequency precession of Haumea's pole caused by Namaka would produce little detectable change in Haumea's light curve or occultation shadow.The Hi'iaka coupled precession has a much higher amplitude, but has a period of 100s of years, severely hampering detectability.Given Haumea's prolate shape, the satellites' torques can also alter Haumea's rotation period, however, this effect is tiny due to Haumea's large angular momentum.Using SPINNY simulations, we estimate the period variations are ∼ 10 −8 hours, approximately two orders of magnitude smaller than the uncertainty in the measured rotation period (Rabinowitz et al. 2006).More amenable to detectability is possible precession of Hi'iaka's rotational axis.In Hastings et al. (2016), the light curve of Hi'iaka was studied using resolved photometry from HST images.They found that Hi'iaka is rapidly rotating (∼9.8 hour double-peaked period) with an unusual sawtooth shaped light curve of amplitude ∆m ≈ 0.23.Using a simplified model of axial precession, they found that Hi'iaka's axial precession would be detectable on decade timescales if there was significant obliquity (w.r.t.its orbit).The detectability is significantly enhanced by Hi'iaka's high amplitude light curve.Detection of any precession would require longterm monitoring of Hi'iaka's light curve which would slowly change amplitude and/or shape across the precession cycle.SPINNY provides an ideal framework for validating this possible method.Using the best fit from our nominal orbit model, we have explicitly modeled Hi'iaka's axial precession.In Figure 6, we show the evolution of Hi'iaka's rotation axis assuming differing starting obliquities.Then in Figure 7, we illustrate how that precession translates to variation in Hi'iaka's light curve amplitude assuming triaxial shapes as in Hastings et al. (2016).Interestingly, even in the case where Hi'iaka's pole is initially aligned with its orbit, precession still occurs.While initially surprising, the precession is due to Hi'iaka's nodal precession in its orbit around Haumea which will always misalign Hi'iaka's spin pole.Encouragingly, even for a relatively small obliquity of 10 • , the precession is substantially different from the no precession case.This allows us to confirm previous results (e.g.Hastings et al. 2016) and show that small perturbations (e.g., Hi'iaka's eccentric orbit, torques from Namaka, etc.) seem to make little difference to the overall evolution.In the future, SPINNY and/or MultiMoon could be modified to explicitly model changes in light curve amplitudes.This method would provide a detailed model with which to understand the spin dynamics of Hi'iaka; we defer this to future work.
Even though Namaka has been solidly detected in several epochs of HST observations, no periodic brightness variations have been found, although photometry from HST programs is suggestive of a long rotation period (>1 day).Purely based on theoretical dynamical arguments, Namaka is expected to be significantly despun, except if its initial rotation period was extremely short (Hastings et al. 2016).Given Namaka's eccentric orbit, overlap between SORs inevitably causes spin-orbit chaos, very similar to Saturn's satellite Hyperion (Wisdom et al. 1984).As an illustration of chaotic rotation, the resonance overlap criterion, which predicts chaotic spin-orbit evolution if satisfied, near the 1:1 and 3:2 SORs is: where A, B, and C are Namaka's principal moments of inertia and e is Namaka's orbital eccentricity.In this case, to avoid spin-orbit chaos Namaka's shape would have to satisfy (B−A) C < 0.025.Satellites that are similar in size to Namaka generally have (B−A) C ≳ 0.1.For example, with a mass between Namaka and Hi'iaka (∼5×10 18 kg), Hyperion has (B−A) C ≈ 0.26.While the above calculation is a simplistic example comparing just two resonances, in general, resonance overlap and spin-orbit chaos are expected for a slowly rotating Namaka.In Figure 8, we show the initiation of chaotic tumbling where Namaka's attitude (i.e.pole direction) and rotation period evolve chaotically.Given the difficulty of acquiring resolved photometric observations of Namaka and the long timespan needed to detect a slow (and possibly chaotic) rotation, confirming it may be extremely difficult.Similarly, searching for Namaka's light curve in unresolved photometry of the entire system is extremely difficult as Namaka's brightness only contributes ∼1% of the total system flux.

Ring-Satellite Interactions
Haumea's ring, first discovered during a stellar occultation, lies close to Haumea's equatorial plane, as shown in the previous section.This matches theoretical expectations, which show that ring particles should collisionally damp to the equatorial plane in the presence of Haumea's J 2 .When accounting for interactions between ring particles and Haumea's satellites, however, ring particles can preferentially settle in the satellite orbit plane (Marzari 2020).However, the satellites are too far away and/or not massive enough to cause this to occur for the ring.The satellites' main dynamical roles are to act as small perturbers that increase velocity dispersion.Marzari (2020)  .The precession of Hi'iaka's spin axis based on its initial obliquity.Similar to Figure 5, we show the precession of Hi'iaka's spin axis in terms of ecliptic longitude and latitude.Initial integration parameters have been chosen to be representative of the posterior found in Table 2.The moments of inertia of Hi'iaka were chosen to be similar to objects of similar size, although their values only change the frequency of the precession.For different initial obliquities, we find different precession periods, matching analytical theory and results in the literature (e.g.Hastings et al. 2016).Interestingly, there are small variations when Hi'iaka's spin is initially aligned with its orbit.Although initially there would be no net torque and no precession when aligned, since Hi'iaka's orbit precesses, the alignement is broken and precession begins.
collision velocity of ring particles under the influence of Haumea's satellites and found that collision velocities still remained low with typical velocities < 1 m s −1 .
It may be possible that if other rings exist external to the currently known ring, satellite-ring dynamics may be more important.As the strength of satellite-ring interactions increases, the collision velocities between ring particles may become large enough to completely disperse the ring.Rings external to the Roche limit may be possible at the temperatures of the Kuiper belt (Morgado et al. 2023;Pereira et al. 2023), but in Haumea's case they are likely to be in a regime where perturbations make them long-term unstable.

Future and past observations
To enable future observations of the Haumea system, we have created an ephemeris of the system containing the predicted ∆α cos δ and ∆δ of each satel-lite and the uncertainties on the positions.We compute the ephemeris using 500 random draws from our posterior distribution.The ephemeris contains the predicted positions every 8 hours between 2005 and 2035.In Table 3, we show the first 10 lines of the ephemeris.The ephemeris is published in its entirety in machinereadable format.Uncertainties on our predictions for both satellites are quite small (≲50 mas) through the 2020s until 2030, at which time, the uncertainties begin to grow rapidly, especially for Namaka.By 2035, uncertainties on Namaka's position are of order ∼0.1 arcseconds.Rapid growth in uncertainty is attributable to the large degeneracy in our model (see Table 2 and Figure 3).
To ascertain whether future observations may be able to break the mass-J 2 degeneracy in our fits, we have analyzed an ephemeris (similar to that presented above) where the predicted positions are also a function of Haumea's J 2 (or equivalently Hi'iaka's mass).We find that the future positions of Hi'iaka and Namaka are After just a few decades, Namaka becomes attitude unstable and begins to tumble.Chaos is a natural consequence of Namaka's eccentric orbit and is inevitable if Namaka has been substantially despun.In this sense, Namaka is very similar to Saturn's satellite Hyperion, which chaotically rotates due to its eccentric orbit around Saturn (Wisdom et al. 1984;Klavetter 1989).
strong functions of Haumea's J 2 (or Hi'iaka's mass) indicating that the degeneracy will be broken with additional HST observations.Based on our fits, we can roughly estimate that ∼2-4 new epochs of observations in the mid-2020s (or beyond) are necessary to constrain Haumea's J 2 with ∼10% precision.Past work has shown that timing observations at certain, well-selected times can substantially improve the quality of the resulting orbit fits (Proudfoot et al., submitted).These times occur when the uncertainties in ∆α cos δ and ∆δ in Table 3 are at (local) maximum.We recommend that continued observations be taken by HST to prevent any systematic errors from arising, as we found in our HST+Keck fits.
As these observations will probe Haumea's interior and place strong constraints on its differentiation, we regard them as high priority.
As well as its use for future observations, our ephemeris can evaluate past predictions of mutual events in light of our new orbit solution.We find that the mutual event predictions made in RB09 are generally accurate, even given their sign error in the system astrometry.We predict events that are similar in depth and duration, but are somewhat offset in time.The timing differences are only a few hours in 2009-2011, but steadily grow to tens of hours by the end of the mutual event season.Since the mutual event timings are quite sensitive to Namaka's eccentricity, and we find an eccentricity ∼15% lower than previously found, it is unsurprising that we find timing differences.The next mutual event season will occur in approximately half a heliocentric Haumea orbit, about midway through the 2200s.The exact timeframe will be dependent on Namaka and Hi'iaka's precise precession rates which future observations will be able to precisely measure.

The Origin and Evolution of Haumea's Satellites
Our new information on the properties of the Haumea system give more accurate insight into the formation and evolution of Haumea and its satellites.
Presumably, Haumea's satellites (and possibly its ring) formed in the same event that created the Haumea family.For a full review of family formation proposals, see Proudfoot & Ragozzine (2019).The most compelling proposals involve mass shedding due to an excess of angular momentum, explaining Haumea's rapid rotation and the family members' low ejection velocity.There is still much debate in the literature on the original reason that proto-Haumea had too much angular momentum: Proudfoot & Ragozzine (2022) propose the merger of a near-equal-mass binary similar to Pluto-Charon perhaps triggered by Kozai cycles initiated due to Haumea's placement on a higher inclination orbit.Noviello et al. (2022) point out that internal evolution could change Haumea's moment of inertia, spinning it up past the point of break-up.Ortiz et al. (2012) and others suggest that it may be the cumulative effect of many smaller impacts, though starting with a rapid rotation would significant increase the probability that a random walk would lead to such a rapid rotation.
Whatever the origin, the general agreement is that Haumea goes through a phase of excess angular momentum where water ice chunks are ejected at low velocities (relative to catastrophic collisions) from the tips of Haumea.Thus, a plausible starting point for the formation of the satellites is a disk of satellitesimals ejected from the fast-spinning Haumea, mostly composed of pure water ice with sizes similar to known satellites and family members.Given that the distribution of ejection (or sub-ejection) velocities should be somewhat smooth, it is likely that satellitesimals are both ejected and remain in orbit.
This initial disk of satellitesimals should rapidly evolve through ejections and collisions.In this scenario, the unejected material experiences a collisional cascade that leads to a final configuration of a ring of near-circular, near-coplanar disk of collisional debris.This debris would eventually reaccrete into Hi'iaka and Namaka.Haumea's ring could also derive from the parts of this initial ring that did not form into satellites.Combining smoothed particle hydrodynamic (SPH) modeling with long-term dynamical evolution is necessary to understand this process.We strongly encourage work on this problem, which may provide understanding applicable to the formation of other TNO satellites.
After their initial formation, a variety of physical processes can influence the evolution of the satellites until they reach their current configuration.The most important effects are expected to be tidal evolution, Hi'iaka-Namaka interactions, excitation from passing TNOs and binaries, and possible interactions from previous satellites ( Ćuk et al. 2013;Hastings et al. 2016).These are discussed in detail in Ćuk et al. (2013) and we focus here only things that are updated in our new fit.Ćuk et al. (2013) found that long-term orbital stability would be significantly improved if the satellites were ∼50% of their nominal masses reported in Ragozzine & Brown (2009).Indeed, our new results are most consistent with satellites that are ∼60% of the initially estimated masses with the Namaka/Haumea mass ratio of 3.0 ± 0.6 × 10 −4 and the Hi'iaka/Haumea mass ratio of 3.1 ± 0.8 × 10 −3 in the HST-only fit (see Table 2), though this is strongly affected by the degeneracy discussed above.Ćuk et al. (2013) also investigate in detail the effect of the 8:3 and 3:1 mean-motion resonances between the satellites, especially in the presence of tidal evolution.This was based on the observation by Ragozzine & Brown (2009) that the satellites were possibly in or near the 8:3 resonance, suggesting that tidal evolution in resonances could explain the source of the moderately eccentric and inclined orbits.This idea was further strengthened by the observation -still true with the new orbit fit -that the excitation is similar in eccentricity and inclination and is inversely proportional to the masses of the satellites.This means that the "Angular Momentum Deficit" (Laskar 1997) is approximately evenly partitioned between the two satellites, suggesting they have been strongly dynamically coupled at some point in the past (or present).We find a period ratio of Hi'iaka and Namaka of 2.689 ± 0.004 (including corrections to Newton's Version of Kepler's Third Law from J 2 ), slightly closer to the 8:3 mean-motion resonance that reported in Ragozzine & Brown (2009).With the residual uncertainty in Hi'iaka's mass and Haumea's J 2 , we leave to future investigation whether the 8:3 resonance is currently dynamically active in the system.Confirmation of an active resonance is key to understanding the recent history of the satellites.Overall, however, the general conclusion of ( Ćuk et al. 2013) that resonance passage could explain the excited orbits remains consistent with the updated fits.The primary challenge in explaining the current orbital configuration of the satellites is their distant orbits from Haumea at ∼36 and ∼70 primary radii.Tidal evolution to such distances is challenging in standard tidal theories, requiring Haumea tides to be extremely dissipative with an implausible combination of tidal parameters.This is exacerbated by a factor of ∼2 with the lower masses for the satellites, but Quillen et al. (2016) find that the triaxial shape of Haumea increases tidal evolution by a factor of a few (though not as large as hoped for by Ragozzine & Brown (2009) and Ćuk et al. (2013)).It is possible that more realistic tidal and geophysical models are able to resolve these issues.It is interesting to note that the satellites are close to the expected outcome of tidal evolution given their mass and semi-major axis ratios.
One potential resolution to the tidal dissipation problem is to form Hi'iaka and Namaka near their current locations.Since tidal expansion is very strongly dependent on separation, even starting at ∼90% of the present distance does not relieve pressure on tidal theories.Given that we observe objects ejected from this disk, it could have been extended out to the Hill sphere.A disk of such size may allow satellite formation further out than previously thought.Along these lines, we note that although the satellites seem well-separated, dynamically speaking they are only separated about 5 mutual Hill radii, suggesting they are dynamically packed.Hence, intermediate satellites could not fit dynamically, so perhaps Namaka and Hi'iaka are the natural outcome of an extended disk near their present locations.Further modeling is needed to understand the plausibility of this scenario.
In conclusion, the formation and evolution of Namaka and Hi'iaka are plausibly connected to the same process that spun up Haumea and created the family.One possible formation hypothesis is that water ice chunks which do not escape to form the family eventually collide and grind down to a disk near the present location of the satellites.Namaka and Hi'iaka perhaps form directly from this disk and recent dynamical interactions (e.g., from the nearby 8:3 resonance) lead to the orbits seen today, as proposed in Ćuk et al. (2013).Once Hi'iaka's mass is better known, a more detailed investigation into the secular, resonant, and tidal dynamics could confirm or refute this hypothesis.However, the most important next step is more detailed modeling of the post-spinup and family ejection process, extending to the longer timescals needed to understand the subsequent epoch of satellite formation.

CONCLUSION
Using a state-of-the-art orbit fitter, MultiMoon, combined with several new epochs of observations from Keck and HST, we have refit the orbits of Haumea's satellites.The model we use can account for both satellite-satellite interactions and Haumea's oblate gravitational field.We find that unaccounted systematic errors are present when fitting to the combined HST and Keck datasets, even when using robust statistical techniques that can account for some types of systematics.Although the HST+Keck fit can precisely constrain Haumea's J 2 and the masses of Hi'iaka and Namaka, we reject this fit since it has unreasonably high residuals and predict physically unrealistic values for Haumea's J 2 .On the other hand, our orbit fit to only the HST data provides a better fit to the data overall.Unfortunately, this fit suffers from a degeneracy between Hi'iaka's mass and Haumea's J 2 , preventing a precise measurement of these two parameters.
For our HST only orbit fit, we detect Haumea's J 2 at ∼ 2.5σ confidence (J 2 = 0.262 +0.103  −0.112 ).Our fits are unable to discriminate between either a homogeneous or differentiated interior, but only a few additional epochs of precise astrometric observations will easily provide the precision to distinguish between these models.Our fit has also provided a measurement of Haumea's rotational pole (α p , δ p ) = (282.9+0.6 −0.7 • , −9.7 +0.6 −1.0 • ), which lies extremely close to the orbit pole of Haumea's ring (Ortiz et al. 2017).In this sense, we presume that Haumea's ring lies in Haumea's equatorial plane and is minimally perturbed by Hi'iaka and Namaka.Determining Haumea's pole allows us to place tight constraints on the inclination of the satellites w.r.t.Haumea's equator, showing that Hi'iaka and Namaka are inclined by approximately 1.01 +0.66 −0.47 • and 12.79 +1.01 −0.58 • , respectively.Both Hi'iaka and Namaka are on somewhat excited orbits, shown in both their inclination and eccentricity, hinting at past dynamical excitation ( Ćuk et al. 2013).
Using our orbit fits, we have also characterized the rotational dynamics of the Haumea system using the spin-orbit integrator SPINNY (Ragozzine et al., submitted).We find that Haumea's rotation axis precesses < 0.5 • on ∼kyr timescales, and is most strongly coupled to Hi'iaka's slow precession due to Haumea's J 2 .Hi'iaka's rotation is expected to strongly precess on decadal timescales, which should have strong effects on the evolution of its light curve.Namaka is expected to rotate extremely slowly, based on both dynamical/tidal arguments and preliminary studies of its light curve.This putative slow rotation implies that Namaka chaotically rotates due to its significantly eccentric orbit.
To enable future observations of the Haumea system we have generated a satellite ephemeris over the next decade.These observations will enable a probe of Haumea's interior, aid in understanding the spin states of Haumea's satellites, and continue to provide insights into Haumea's formation.Understanding the Haumea system as a whole is crucial for understanding large TNO formation and evolution.The production of satellites and satellite systems seems to be ubiqitous across the transneptunian region, but the processes at play are still not well-explored.Thankfully, continued observa-tions of Haumea and its satellites will enable deeper knowledge of the far reaches of our solar system.

Figure 2 .
Figure2.The normalized residuals of the best parameter set from our HST-only (non-robust) orbit fit.x and y correspond to ecliptic longitude and latitude, which is the primary coordinate system used in MultiMoon.Hi'iaka residuals are shown as circles and Namaka with triangles.The color of each point corresponds to the observation date.The circles correspond to 1, 2, and 3σ error contours.As reported in the text, this fit corresponded to a reduced χ 2 ∼ 1.10.We find that although the fit quality is worse than would be desired, the p-value of associated with the χ 2 statistic is 0.23, indicating an acceptable fit.

Figure 3 .
Figure3.The joint Haumea J2-Hi'iaka mass posterior distribution.Instead of displaying the density of sampled points as in Figure1, we show the maximum fit quality (as measured by reduced χ 2 ) in a small bin.Bins without any sampled points, indicating extremely poor quality fits, were set to the minimum bin value, although the true value is likely much worse.We find that the best fit with χ 2 ν ∼ 1.1 has a p-value of 0.23, meaning there is a 23% chance that random chance would produce a worse fit.A χ 2 ν ∼ 1.25 corresponds to a p-value of 0.05.Dashed red lines show the expected J2 values from different internal density models.The undifferentiated model assumes a homogeneous interior along with the occultation derived shape model(Ortiz et al. 2017).The differentiated model is a two layer model proposed byDunham et al. (2019).The posterior shows that both models are consistent with the data, although the differentiated model is slightly disfavored.

Figure 4 .
Figure4.Change in orbit fit quality due to the nonspherical shapes of Haumea's satellites.In this plot, ∆χ 2 = i (xi,J 2 − xi,J 2 =0)/σi,obs, where xi,J 2 and xi,J 2 =0 are synthetic astrometric measurements from models including satellite J2 and those without.σ i,obs are the true measurement uncertainties for the system as tabulated in Table1.This formulation allows us to directly determine whether the satellites' J2values produce detectable changes in the system relative astrometry, given our current observations (HST+Keck).The parameters for the models used were taken from the best fit of our orbit fits.The rotational poles of the satellites were chosen to have high obliquities (w.r.t.their Haumea-centric orbits) to enhance the effect of J2.The three separate lines show each satellite's individual contribution, as well as a comparison where both satellites have (the same) J2.For reasonable values, J2 ≲ 0.5, very little change in fit quality is found, although for extremely nonspherical shapes, Hi'iaka's J2 could begin to alter the fit.
Figure5.The spin precession of Haumea over a 1000 year integration.In the top two panels, we show the spin precession of Haumea in terms of pole ecliptic longitude and latitude.In the bottom two panels we show the precession of the orbit normal of Hi'iaka and Namaka, again in terms of ecliptic longitude and latitude.The integration parameters have been chosen to be representative of the posterior found in Table2.Similar integrations with different values for Haumea's J2 change the long-term precession period, but are qualitatively similar.Haumea's spin precession is coupled with the nodal precession of the satellites.High frequency, low amplitude variations in Haumea's pole direction are caused by Namaka's rapid precession, while low frequency, high amplitude variations are coupled with Hi'iaka's J2 precession.As can be seen, the precession of all components are strongly coupled, both through Hi'iaka-Namaka gravitational interactions and interactions with Haumea's J2.
Figure6.The precession of Hi'iaka's spin axis based on its initial obliquity.Similar to Figure5, we show the precession of Hi'iaka's spin axis in terms of ecliptic longitude and latitude.Initial integration parameters have been chosen to be representative of the posterior found in Table2.The moments of inertia of Hi'iaka were chosen to be similar to objects of similar size, although their values only change the frequency of the precession.For different initial obliquities, we find different precession periods, matching analytical theory and results in the literature (e.g.Hastings et al. 2016).Interestingly, there are small variations when Hi'iaka's spin is initially aligned with its orbit.Although initially there would be no net torque and no precession when aligned, since Hi'iaka's orbit precesses, the alignement is broken and precession begins.

Figure 8 .
Figure8.Chaotic rotation of Namaka.When initialized, Namaka rotates once every 6 days and is aligned with its orbit.After just a few decades, Namaka becomes attitude unstable and begins to tumble.Chaos is a natural consequence of Namaka's eccentric orbit and is inevitable if Namaka has been substantially despun.In this sense, Namaka is very similar to Saturn's satellite Hyperion, which chaotically rotates due to its eccentric orbit around Saturn(Wisdom et al. 1984;Klavetter 1989).