ORBITS AND MASSES OF THE SATELLITES OF THE DWARF PLANET HAUMEA (2003 EL61)

Determining the orbits and masses of the full system requires a 15-dimensional, highly nonlinear, global χ² minimization. We found this to be impractical without a good initial guess for the orbit of Namaka to reduce the otherwise enormous parameter space, motivating the acquisition of multiple observations within a short enough timescale that Namaka's orbit is essentially Keplerian. Fitting the 2008 May HST data with a Keplerian model produced the initial guess necessary for the global minimization of the fully interacting three-point-mass model. The three-point-mass model uses 15 parameters: the masses of the Haumea, Hi'iaka, and Namaka, and, for both orbits, the osculating semimajor axis, eccentricity, inclination, longitude of the ascending node, argument of periapse, and mean anomaly at epoch HJD 2454615.0 (= May 28.5, 2008). All angles are defined in the J2000 ecliptic coordinate system. Using these orbital elements, we constructed the Cartesian locations and velocities at this epoch as an initial condition for the three-body integration. Using a sufficiently small time step (∼300 s for the final iteration), a FORTRAN 90 program integrates the system to calculate the positions relative to the primary and the positions at the exact times of observation are determined by interpolation. (Observation times were converted to Heliocentric Julian Dates, the date in the reference frame of the Sun, to account for light-travel time effects due to the motion of the Earth and Haumea, although ignoring this conversion does not have a significant effect on the solution.) Using the JPL HORIZONS ephemeris for the geocentric position of Haumea, we vectorially add the primary-centered positions of the satellites and calculate the relative astrometric on-the-sky positions of both satellites. This model orbit is then compared to the data by computing χ² in the normal fashion. We note here that this model does not include gravitational perturbations from the Sun or center-of-light/center-of-mass corrections, which are discussed below.

Like many multidimensional nonlinear minimization problems, searching for the best-fitting parameters required a global minimization algorithm to escape the ubiquitous local minima. Our algorithm for finding the global minimum starts with thousands of local minimizations, executed with the Levenberg–Marquardt algorithm mpfit³ using numerically determined derivatives. These local minimizations are given initial guesses that cover a very wide range of parameter space. Combining all the results of these local fits, the resultant parameter versus χ² plots showed the expected parabolic shape (on scales comparable to the error bars) and these were extrapolated to their minima. This process was iterated until a global minimum is found; at every step, random deviations of the parameters were added to the best-fit solutions, to ensure a full exploration of parameter space. Because many parameters are highly correlated, the ability to find the best solutions was increased significantly by adding correlated random deviations to the fit parameters as determined from the covariance matrix of the best known solutions. We also found it necessary to optimize the speed of the evaluation of χ² from the 15 system parameters; on a typical fast processor this would take a few hundredths of a second and a full local minimization would take several seconds.

To determine the error bars on the fit parameters, we use a Monte Carlo technique (B05), as suggested in Press et al. (1992). Synthetic data sets are constructed by adding independent Gaussian errors to the observed data. The synthetic data sets are then fit using our global minimization routine, resulting in 86 Monte Carlo realizations; four of the synthetic data-sets did not reach global minima and were discarded, having no significant affect on the error estimates. One-sigma one-dimensional error bars for each parameter are given by the standard deviation of global-best parameter fits from these synthetic data sets. For each parameter individually, the distributions were nearly Gaussian and were centered very nearly on the best-fit parameters determined from the actual data. The error bars were comparable to error bars estimated by calculating where χ² increased by 1 from the global minimum (see Press et al. 1992).

First, we consider a solution using only the observations from HST. Even though these are taken with different instruments (ACS, NICMOS, and mostly WFPC2), the extensive calibration of these cameras allows the direct combination of astrometry into a single data set. The best-fit parameters and errors are shown in Table 2 and the on-the-sky orientation of the orbits in 2008 March is shown in Figure 1. The reduced χ² for this model is χ²_red = 0.64 (χ² = 36.4 with 57 degrees of freedom). The data are very well fit by the three-point-mass model, as shown in Figures 2 and 3. A reduced χ² less than 1 is an indication that error bars are overestimated, assuming that they are independent; we note that using 10 separate "observations" for the 2007 February data implies that our observations are not completely independent. Even so, χ²_red values lower than 1 are typical for this kind of astrometric orbit fitting (e.g., Grundy et al. 2009). Each of the fit parameters is recovered, though the mass of Namaka is only detected with a 1.2σ significance. Namaka's mass is the hardest parameter to determine since it requires detecting minute non-Keplerian perturbations to orbit of the more massive Hi'iaka. The implications of the orbital state of the Haumea system are described in the following section. We also list in Table 3 the initial condition of the three-body integration for this solution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The HST data are sufficient to obtain a solution for Hi'iaka's orbit that is essentially the same as the orbit obtained from the initial Keck data in B05. Nevertheless, the amount and baseline of Keck NIRC2 data is useful enough to justify adding this data set to the fit. Simply combining these data sets and searching for the global minimum results in a significant degradation in the fit, going from a reduced χ² of 0.64 to a reduced χ² of ∼1.10, although we note that this is still an adequate fit. Adding the Keck data has the effect of generally lowering the error bars and subtly changing some of the retrieved parameters. Almost all of these changes are within the ∼1σ error bars of the HST only solution, except for the mass estimate of Namaka. Adding the Keck data results in a best-fit Namaka mass a factor of 10 lower than the HST data alone. The largest mass retrieved from the entire Monte Carlo suite of solutions to the HST+Keck data set is ∼8 × 10¹⁷ kg, i.e., a Namaka/Haumea mass ratio of 2 × 10⁻⁴, which is inconsistent with the brightness ratio of ∼0.02, for albedos less than 1 and densities greater than 0.3 g cm⁻³. However, this solution assumes that the Keck NIRC2 absolute astrometry (based on the solution of Ghez et al. (2008), which is not directly cross-calibrated with HST) is perfectly consistent with HST astrometry. In reality, a small difference in the relative plate scale and rotation between these two telescopes could introduce systematic errors. Adding fitted parameters that adjust the plate scale and rotation angle does not help, since this results in overfitting, as verified by trial fitting of synthetic data sets. We adopt the HST-only solution, keeping in mind that the nominal mass of Namaka may be somewhat overestimated.

Using the Monte Carlo suite of HST-only solutions, we can also calculate derived parameters and their errors. Using Kepler's Law (and ignoring the other satellite), the periods of Hi'iaka and Namaka are 49.462 ± 0.083 days and 18.2783 ± 0.0076 days, respectively, with a ratio of 2.7060 ± 0.0037, near the 8:3 resonance. The actual mean motions (and resonance occupation) will be affected by the presence of the other satellite and the nonspherical nature of the primary (discussed below).

The mass ratios of the satellite to Haumea are 0.00451 ± 0.00030 and 0.00051 ± 0.00036, respectively, and the Namaka/ Hi'iaka mass ratio is 0.116 ± 0.086. The mutual inclination of the two orbits is ϕ = 1341 ± 008, where the mutual inclination is the actual angle between the two orbits, given by cos ϕ = cos i_Hcos i_N + sin i_Hsin i_Ncos(Ω_H − Ω_N), where i and Ω are the inclination and longitude of ascending node. The origin of this significantly nonzero mutual inclination is discussed in Section 4.3.2. The mean longitude λ ≡ Ω +  ω +  M, is the angle between the reference line (J2000 ecliptic first point of Ares) and is determined well; the errors in the argument of periapse (ω) and mean anomaly (M) shown in Tabel 2 are highly anti-correlated. Our Monte Carlo results give λ_H = 153.80 ± 0.34 degrees and λ_N = 202.57 ± 0.73 degrees. Finally, under the nominal point-mass model, Namaka's argument of periapse changes by about −65 per year during the course of the observations, implying a precession period of about 55 years; the non-Keplerian nature of Namaka's orbit is detected with very high confidence.

The nonspherical nature of Haumea can introduce additional, potentially observable, non-Keplerian effects. The largest of these effects is due to the lowest-order gravitational moment, the quadrupole term (the dipole moment is 0 in the center of mass frame), described by J₂ (see, e.g., Murray & Dermott 2000). Haumea rotates over 100 times during a single orbit of Namaka, which orbits quite far away at ∼35 primary radii. To lowest order, therefore, it is appropriate to treat Haumea as having an "effective" time-averaged J₂. Using a code provided by E. Fahnestock, we integrated trajectories similar to Namaka's orbit around a homogeneous rotating tri-axial ellipsoid and have confirmed that the effective J₂ model deviates from the full model by less than half a milliarcsecond over three years.

The value of the effective J₂ (≡−C₂₀) for a rotating homogeneous tri-axial ellipsoid was derived by Scheeres (1994)

$$\begin{equation} J_2 R^2 = \case{1}{10}(\alpha ^2+\beta ^2-2\gamma ^2) \simeq 1.04 \times 10^{11}\; \hbox{m}^2, \end{equation} \tag{ 1 }$$

where α, β, and γ are the tri-axial radii and the numerical value corresponds to a (498 × 759 × 980) km ellipsoid as inferred from photometry (Rabinowitz et al. 2006). We note that the physical quantity actually used to determine the orbital evolution is J₂R²; in a highly tri-axial body like Haumea, it is not clear how to define R, so using J₂R² reduces confusion. If R is taken to be the volumetric effective radius, then R ≃ 652 km and the J₂ ≃ 0.244. Note that the calculation and use of J₂ implicitly requires a definition of the rotation axis, presumed to be aligned with the shortest axis of the ellipsoid.

Preliminary investigations showed that using this value of J₂R² implied a non-Keplerian effect on Namaka's orbit that was smaller, but similar to, the effect of the outer satellite. The primary observable effect of both J₂ and Hi'iaka is the precession of apses and nodes of Namaka's eccentric and inclined orbit (Murray & Dermott 2000). When adding the three relevant parameters—J₂R² and the direction of the rotational axis⁴—to our fitting procedure, we found a direct anti-correlation between J₂R² and the mass of Hi'iaka, indicating that these two parameters are degenerate in the current set of observations. The value of the reduced χ² was lowered significantly by the addition of these parameters, with the F-test returning high statistical significance. However, the best fits placed the satellites on polar orbits. We have verified that the fitted values of J₂R² and the spin pole perpendicular to the orbital poles is due to overfitting of the data. We generated simulated observations with the expected value of J₂R² and with the satellites in nearly equatorial orbits. Allowing the global fitter to vary all the parameters resulted in an overfitted solution that placed the satellites on polar orbits. Hence, allowing J₂ and the spin pole to vary in the fit is not justified; the effect of these parameters on the solution are too small and/or too degenerate to detect reliably. Furthermore, since the model without these parameters already had a reduced χ² less than 1, these additional parameters were not warranted in the first place.

It is interesting, however, to consider how including this effect would change the determination of the other parameters. We therefore ran an additional set of models with a fixed J₂R² and fixing the spin pole (more accurately, the axis by which J₂ is defined) as the mass-weighted orbital pole of Hi'iaka and Namaka; since Hi'iaka is ∼10 times more massive than Namaka, this puts Hi'iaka on a nearly equatorial orbit (i ≃ 1°), as would be expected from collisional formation. Holding J₂R² fixed at 1.04 × 10¹¹ m², we re-analyzed the HST data set using our global fitting routine. As expected, none of the parameters change by more than 1σ, except for the mass of Hi'iaka, which was reduced by almost 30% to ∼1.35 × 10¹⁹ kg. (In fits where J₂R² was allowed to vary, the tradeoff between J₂R² and Hi'iaka's mass was roughly linear, as would be expected if the sum of these effects were forced to match the observationally determined precession of Namaka.) The data were well fit by the non-point-mass model, with the forced J₂R² and spin pole solution reaching a global reduced χ² of 0.72. Since Haumea's high-amplitude light curve indicates a primary with a large quadrupole component, the nominal mass of Hi'iaka in the point-mass case is almost certainly an overestimate of its true mass. More data will be required to disentangle the degeneracy between the mass of Hi'iaka and the J₂ of Haumea. Including J₂ and/or the Keck data do not improve the estimates of Namaka's mass.

According to the orbit solution, the Haumea system is currently undergoing mutual events, as reported in (Fabrycky et al. 2008). (This is also true using the other orbit solutions, e.g., HST+Keck, with or without J₂.) Using the known orbit, the angle between Namaka, Haumea, and the Earth (in the case of occultations) or the Sun (in the case of shadowing) falls well below the ∼13 mas (∼500 km) of the projected shortest axis of Haumea. Observing multiple mutual events can yield accurate and useful measurements of several system properties as shown by the results of the Pluto–Charon mutual event season (e.g., Binzel & Hubbard 1997). The depth of an event where Namaka occults Haumea leads to the ratio of albedos and, potentially, a surface albedo map of Haumea, which is known to exhibit color variations as a function of rotational phase, indicative of a variegated surface (Lacerda et al. 2008; Lacerda 2009). Over the course of a single season, Namaka will traverse several chords across Haumea allowing for a highly accurate measurement of Haumea's size, shape, and spin pole direction (e.g., Descamps et al. 2008). The precise timing of mutual events will also serve as extremely accurate astrometry, allowing for an orbital solution much more precise than reported here. We believe that incorporating these events into our astrometric model will be sufficient to independently determine the masses of all three bodies and J₂R². Our solution also predicts a satellite–satellite mutual event in 2009 July—the last such event until the next mutual event season begins around the year 2100. Our knowledge of the state of the Haumea system will improve significantly with the observation and analysis of these events. See http://web.gps.caltech.edu/~mbrown/2003EL61/mutual for up-to-date information on the Haumea mutual events. Note that both the mutual events and the three-body nature of the system are valuable for independently checking the astrometric analysis, e.g., by refining plate scales and rotations.

Using the photometry and the best-fit mass ratio found in Sections 2 and 3, we can estimate the range of plausible albedos and densities for the two satellites. The results of this calculation are shown in Figure 4. The mass and brightness ratios clearly show that the satellites must either have higher albedos or lower densities than Haumea; the difference is probably even more significant than shown in Figure 4 since the nominal masses of Hi'iaka and Namaka are probably overestimated (see Section 3). The similar spectral (Barkume et al. 2006) and photometric (Fraser & Brown 2009) properties of Haumea, Hi'iaka, and Namaka indicate that their albedos should be similar. Similar surfaces are also expected from rough calculations of ejecta exchange discussed by Stern (2009), though Benecchi et al. (2009) provide a contrary viewpoint. If the albedos are comparable, the satellite densities indicate a mostly water ice composition (ρ ≈ 1.0 g cm⁻³). This lends support to the hypothesis that the satellites are formed from a collisional debris disk composed primarily of water ice from the shattered mantle of Haumea. This can be confirmed in the future with a direct measurement of Namaka's size from mutual event photometry. Assuming a density of water ice, the estimated radii of Hi'iaka and Namaka are ∼160 km and ∼80 km, respectively.

Zoom In Zoom Out Reset image size

It is surprising to find the orbits in an excited state, both with nonzero eccentricities and with a rather large mutual inclination. In contrast the regular satellite systems of the gas giants, the satellites of Mars, the three satellites of Pluto (Tholen et al. 2008), and asteroid triple systems with well-known orbits (Marchis et al. 2005) are all in nearly circular and coplanar orbits. In systems of more than one satellite, perturbations between the satellites produce forced eccentricities and inclinations that will remain even with significant damping. If the excited state of the Haumea system is just a reflection of normal interactions, then there will be small free eccentricities and inclinations, which can be estimated by integrating the system for much longer than the precession timescales and computing the time average of these elements. Using this technique, and exploring the entire Monte Carlo suite of orbital solutions, we find that the free eccentricity of Hi'iaka is ∼0.07, the free eccentricity of Namaka is ∼0.21, and the time-averaged mutual inclination is ∼125. Nonzero free eccentricities and inclinations imply that the excited state of the system is not due to satellite–satellite perturbations. These integrations were calculated using the n-body code SyMBA (Levison & Duncan 1994) using the regularized mixed variable symplectic method based on the mapping by Wisdom & Holman (1991). Integration of all the Monte Carlo orbits showed for ∼2000 years showed no signs of instability, though we do note that the system chaotically enters and exits the 8:3 resonance. The orbital solutions including J₂R² ≃ 1.04 × 10¹¹ m² were generally more chaotic, but were otherwise similar to the point-mass integrations.

These integrations did not include the effect of the Sun, which adds an additional minor torque to the system that is negligible (ΔΩ ∼ 10⁻⁵ degrees) over the timescale of observations. The effects of the Sun on the satellite orbits on long timescales were not investigated. While the relative inclination between the satellite orbits and Haumea's heliocentric orbit (∼119° for Hi'iaka and ∼105° for Namaka) places this system in the regime where the Kozai effect can be important (Kozai 1962; Perets & Naoz 2008), the interactions between the satellites are strong enough to suppress weak Kozai oscillations due to the Sun, which are only active in the absence of other perturbations (Fabrycky & Tremaine 2007).

We did not include any correction to our solution for possible differences between the center of light (more precisely, the center of fitted PSF) and center of mass of Haumea. This may be important since Lacerda et al. (2008) find that Haumea's two-peaked light curve can be explained by a dark red albedo feature, which could potentially introduce an systematic astrometric error. The 2007 February Hubble data and 2008 March Keck data, both of which span a full rotation, do not require a center-of-light/center-of-mass correction for a good fit, implying that the correction should be smaller than ∼2 mas (i.e., ∼70 km). Examination of all the astrometric residuals and a low reduced-χ² confirm that center-of-light/center-of-mass corrections are not significant at our level of accuracy. For Pluto and Charon, albedo features can result in spurious orbital astrometry because Pluto and Charon are spin-locked; this is not the case for Haumea. For Keck observations where Namaka is not detected (see Table 1), it is usually because of low signal to noise and Namaka's calculated position is not near Haumea. In the cases where Namaka's light contaminates Haumea, the induced photocenter error would be less than the observed astrometric error.

All of the available evidence points to a scenario for the formation of Haumea's satellites similar to the formation of the Earth's moon: a large oblique collision created a disk of debris composed mostly of the water ice mantle of a presumably differentiated proto-Haumea. Two relatively massive moons coalesced from the predominantly water–ice disk near the Roche lobe. Interestingly, in studying the formation of Earth's Moon, about one third of the simulations of Ida et al. (1997) predict the formation of two moonlets with the outer moonlet ∼10 times more massive the inner moonlet. Although the disk accretion model used by Ida et al. (1997) made the untenable assumption that the remnant disk would immediately coagulate into solid particles, the general idea that large disks with sufficient angular momentum could result in two separate moons seems reasonable. Such collisional satellites coagulate near the Roche lobe (e.g., a distances of ∼3–5 primary radii) in nearly circular orbits and coplanar with the (new) rotational axis of the primary. For Haumea, the formation of the satellites is presumably concurrent with the formation of the family billions of years ago (Ragozzine & Brown 2007) and the satellites have undergone significant tidal evolution to reach their current orbits (B05).

4.3.1. Tidal Evolution of Semimajor Axes

The equation for the typical semimajor axis tidal expansion of a single-satellite due to primary tides is (Murray & Dermott 2000)

$$\begin{equation} \dot{a}=\frac{3k_{2p}}{Q_p} q \bigg(\frac{R_p}{a} \bigg)^5 na, \end{equation} \tag{ 2 }$$

where k_2p is the second-degree Love number of the primary, Q_p is the primary tidal dissipation parameter (see, e.g., Goldreich & Soter 1966), q is the mass ratio, R_p is the primary radius, a is the satellite semimajor axis, and $n\equiv \sqrt{\frac{GM_{\rm tot}}{a^3}}$ is the satellite mean motion. As pointed out by B05, applying this equation to Hi'iaka's orbit (using the new-found q = 0.0045) indicates that Haumea must be extremely dissipative: Q_p ≃ 17, averaged over the age of the solar system, more dissipative than any known object except ocean tides on the present-day Earth. This high dissipation assumes an unrealistically high k₂ ≃ 1.5, which would be achieved only if Haumea were perfectly fluid. Using the strength of an rocky body and the Yoder (1995) method of estimating k₂, Haumea's estimated k_2p is ∼0.003. Such a value of k_2p would imply a absurdly low and physically implausible Q_p ≪ 1. Starting Hi'iaka on more distant orbits, e.g., the current orbit of Namaka, does not help much in this regard since the tidal expansion at large semimajor axes is the slowest part of the tidal evolution. There are three considerations, however, that may mitigate the apparent requirement of an astonishingly dissipative Haumea. First, if tidal forcing creates a tidal bulge that lags by a constant time (as in Mignard 1980), then Haumea's rapid rotation (which hardly changes throughout tidal evolution) would naturally lead to a significant increase in tidal evolution. In other words, if Q_p is frequency dependent (as it seems to be for solid bodies, see Efroimsky & Lainey 2007), then an effective Q of ∼16 may be equivalent to an object with a one-day rotation period maintaining an effective Q of 100, the typically assumed value for icy solid bodies. That is, Haumea's higher-than-expected dissipation may be related to its fast rotation. Second, the above calculation used the volumetric radius R_p ≃ 650 km, in calculating the magnitude of the tidal bulge torque, where we note that Equation (2) assumes a spherical primary. A complete calculation of the actual torque caused by tidal bulges on a highly nonspherical body is beyond the scope of this paper, but it seems reasonable that since tidal bulges are highly distance-dependent, the volumetric radius may lead to an underestimate in the tidal torque and resulting orbital expansion. Using R_p ≃ 1000 km, likely an overestimate, allows k_2p/Q_p to go down to 7.5 × 10⁻⁶, consistent both with k_2p ≃ 0.003 and Q_p ≃ 400. Clearly, a re-evaluation of tidal torques and the resulting orbital change for satellites around nonspherical primaries is warranted before making assumptions about the tidal properties of Haumea.

The third issue that affects Hi'iaka's tidal evolution is Namaka. While generally the tidal evolution of the two satellites is independent, if the satellites form a resonance it might be possible to boost the orbital expansion of Hi'iaka via angular momentum transfer with the more tidally affected Namaka. (Note that applying the semimajor axis evolution questions to Namaka requires a somewhat less dissipative primary, i.e., Q_p values ∼8 times larger than discussed above.) Even outside of resonance, forced eccentricities can lead to higher dissipation in both satellites, somewhat increasing the orbital expansion rate. Ignoring satellite interactions and applying Equation (2) to each satellite results in the expected relationship between the mass ratio and the semimajor axis ratio of (Canup et al. 1999; Murray & Dermott 2000)

$$\begin{equation} \frac{m_1}{m_2} \simeq \bigg(\frac{a_1}{a_2} \bigg)^{13/2}, \end{equation} \tag{ 3 }$$

where evaluating the right-hand side using the determined orbits implies that $\frac{m_1}{m_2} \simeq 75.4 \pm 0.4$. This mass ratio is highly inconsistent with a brightness ratio of ∼3.7 for the satellites, implying that the satellites have not reached the asymptotic tidal end-state. Equation (3) is also diagnostic of whether the tidally evolving satellites are on converging or diverging orbits: that the left-hand side of Equation (3) is greater than the right-hand side implies that the satellites are on convergent orbits, i.e., the ratio of semimajor axes is increasing and the ratio of the orbital periods is decreasing (when not in resonance).

4.3.2. Tidal Evolution of Eccentricities and Inclinations