Haumea's Shape, Composition, and Internal Structure

The Kuiper Belt Object (KBO) and dwarf planet Haumea is one of the most intriguing and puzzling objects in the outer solar system. Haumea orbits beyond Pluto, with a semimajor axis of 43.2 au, and is currently near its aphelion distance of ≈51.5 au but is relatively bright at magnitude V = 17.3, due to its large size and icy surface. Haumea's mean radius is estimated to be ≈720 km (Lockwood et al. 2014) to ≈795 km (Ortiz et al. 2017), and its reflectance spectra indicate that Haumea's surface is uniformly covered by close to 100% water ice (Trujillo et al. 2007; Pinilla-Alonso et al. 2009). Haumea is the third-brightest KBO, after the dwarf planets Pluto and Makemake (Brown et al. 2006). Haumea has two small satellites, Hi'iaka and Namaka, which enable a determination of its mass, M_H = 4.006 × 10²¹ kg (Ragozzine & Brown 2009); it is the third or fourth most massive known KBO (after Pluto, Eris, and possibly Makemake). Despite its large size, it is a rapid rotator; from its light curve Haumea's rotation rate is found to be 3.91531 ± 0.00005 hr (Lellouch et al. 2010). This means that Haumea is the fastest-rotating KBO (Sheppard & Jewitt 2002) and is in fact the fastest-rotating large (>100 km) object in the solar system (Rabinowitz et al. 2006). Haumea is also associated with a collisional family (Brown et al. 2007) and is known to have a ring (Ortiz et al. 2017). Based on its rapid rotation and its collisional family, Haumea is inferred to have suffered a large collision (Brown et al. 2007), >3 Gyr ago, based on the orbital dispersion of the family members (Volk & Malhotra 2012).

Haumea is larger than other dwarf planets such as Ceres (radius = 473 km), or satellites such as Dione (radius = 561 km) or Ariel (radius = 579 km), all of which are nominally round. Despite this, Haumea exhibits a reflectance light curve with a very large peak-to-trough amplitude, Δm ≈ 0.28 in 2005 (Rabinowitz et al. 2006), Δm = 0.29 in 2007 (Lacerda et al. 2008), and Δm = 0.32 in 2009 (Lockwood et al. 2014). Since Haumea's surface is spectrally uniform, such an extreme change in brightness can only be attributed to a difference in the area presented to the observer. Haumea has been modeled as a triaxial ellipsoid with axes a > b > c, the c-axis being aligned with the rotation axis. In that case, the change in brightness from peak to trough would be given by

$$\begin{eqnarray}&&{\rm{\Delta }}m=2.5\left[{\mathrm{log}}_{10}\left(\displaystyle \frac{a}{b}\right)-{\mathrm{log}}_{10}\left(\displaystyle \frac{{r}_{1}}{{r}_{2}}\right)\right],\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&{r}_{1}={\left({a}^{2}{\cos }^{2}\phi +{c}^{2}{\sin }^{2}\phi \right)}^{1/2}\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&{r}_{2}={\left({b}^{2}{\cos }^{2}\phi +{c}^{2}{\sin }^{2}\phi \right)}^{1/2},\end{eqnarray} \tag{ 3 }$$

ϕ being the angle between the rotation axis and the line of sight (Binzel et al. 1989). If ϕ = 0°, then Δm = 0, because the same a × b ellipse would be presented to the observer. Instead, Δm is maximized when ϕ = 90°, because then the ellipse presented to the observer would vary between a × c and b × c. In that case, ${\rm{\Delta }}m=2.5\,{\mathrm{log}}_{10}(a/b)$ and the axis ratio is related directly to Δm. Assuming ϕ = 90° in 2009, when Δm = 0.32, one would derive b/a = 0.75. Taking into account the scattering properties of an icy surface, Lockwood et al. (2014) refined this to b/a = 0.80 ± 0.01. Thus, Haumea is distinctly nonspherical and is not even axisymmetric.

Haumea appears to be unique among large solar system objects in having such a distinctly nonaxisymmetric, triaxial ellipsoid shape. Based on its rapid rotation (angular velocity ω = 4.457 × 10⁻⁴ s⁻¹), Haumea is inferred to have assumed a particular shape known as a Jacobi ellipsoid. This is a class of equilibrium shapes assumed by (shearless) fluids in hydrostatic equilibrium when they rotate faster than a certain threshold (Chandrasekhar 1969). For a body with angular velocity ω and uniform density ρ, the axis ratios b/a and c/a of the ellipsoid are completely determined, and the axis ratios are single-valued functions of ω²/(πGρ). For a Jacobi ellipsoid with b/a = 0.806 and Haumea's rotation rate, the density must be ρ = 2580 kg m⁻³, and c/a = 0.520. Assuming a semiaxis of a = 960 km then yields b = 774 km and then c = 499 km; (4π/3)abc ρ exactly matches Haumea's mass. The mean radius of Haumea would be 718 km. Moreover, the cross-sectional area of Haumea would then imply a surface albedo p_V ≈ 0.71–0.84 (Rabinowitz et al. 2006; Lacerda & Jewitt 2007; Lellouch et al. 2010; Lockwood et al. 2014), consistent with the albedo of a water ice surface.

To explain Haumea's icy surface and ρ = 2580 kg m⁻³, one would have to assume that the interior of Haumea was close to 2600 kg m⁻³ in density (an interior of hydrated silicates; Desch & Neveu 2015) while its surface was a very thin ice layer (Section 2). This structure implies that Haumea suffered a giant collision in its past that may have stripped its ice mantle. For these reasons, the above axes and axis ratios were strongly favored in the literature. Other groups derived similar axes and bulk densities (Rabinowitz et al. 2006; Lacerda et al. 2008; Lellouch et al. 2010).

This model was upended by the observations of Ortiz et al. (2017) following the occultation of an 18th magnitude star by Haumea in 2017 January. The shadow of Haumea traced out an ellipse, as expected for the shadow of a triaxial ellipsoid, but the (semi)axes of the shadow ellipse were much larger than expected: b' = 569 ± 13 km by a' = 852 ± 2 km. Ortiz et al. (2017) used the shadow axes and other assumptions to derive the axes of Haumea to be a = 1161 ± 30 km, b = 852 ± 4 km, and c = 513 ± 16 km (the Appendix).

This new shape causes Haumea to look significantly different than previous models: the mean radius of Haumea is larger at 798 km, the albedo is a smaller p_V ≈ 0.51 (and would require a darkening agent in addition to water ice), and the bulk density is lower at 1885 ± 80 kg m⁻³. Moreover, the axis ratio c/a ≈ 0.44 is significantly lower than previous estimates (≈0.52) and is inconsistent with a Jacobi ellipsoid or a fluid in hydrostatic equilibrium. Ortiz et al. (2017) point out the possibility that shear stresses may be supported on Haumea by granular interparticle forces (Holsapple 2001).

In either case, Haumea is likely to have a rocky core surrounded by ice, but no analytical solution exists for the figure of equilibrium of a rapidly rotating, differentiated body. Therefore, it is not known whether or not Haumea is a fluid in hydrostatic equilibrium. In this paper, we attempt to reconcile the existing data from Haumea's light curve and occultation shadow, with the goal of deriving its true shape and internal structure. Besides its axes, important quantities to constrain are the ice fraction on Haumea today and the size, shape, and density of its core. A central question we can solve using these quantities is whether Haumea is a fluid in hydrostatic equilibrium or demands granular physics to support it against shear stresses. In addition, we can use the core density and ice fraction to constrain the geochemical evolution of Haumea and, by extension, other KBOs, as well as models of the origin of Haumea's collisional family.

In Section 2, we examine whether it is possible for Haumea to be a differentiated (rocky core, icy mantle) body with a Jacobi ellipsoid shape. We show that only homogeneous bodies are consistent with a Jacobi ellipsoid shape. In Section 3, we describe a code we have written to calculate the equilibrium figure of a rapidly rotating, differentiated body. In Section 4, we present the results, showing that entire families of solutions exist that allow Haumea to be a differentiated body in hydrostatic equilibrium. Some of these solutions appear consistent with observations of Haumea, particularly ≈1050 × 840 × 537 km with bulk density ≈2018 kg m⁻³. In Section 5, we discuss the implications of this solution for Haumea's structure, for the collision that created the collisional family, and for the astrobiological potential of Haumea.

Observations have suggested that Haumea has axes consistent with a Jacobi ellipsoid: for example, Lockwood et al. (2014) inferred axes of a = 960 km, b = 770 km, and c = 495 km (yielding axis ratios b/a = 0.802 and c/a = 0.516) and a uniform density of ≈2614 kg m⁻³. These axes are within 1% of the Jacobi ellipsoid solution: a Jacobi ellipsoid with Haumea's mass, rotation rate, and a = 960 km has axes b = 774.2 km and c = 498.8 km (yielding axis ratios b/a = 0.806 and c/a = 0.520) and a uniform density of 2580 kg m⁻³. Uniform density is a central assumption of the Jacobi ellipsoid solution, and yet Haumea is manifestly not uniform in density. Its reflectance spectra robustly show the existence of a uniform water ice surface (Trujillo et al. 2007; Pinilla-Alonso et al. 2009). The density of this ice, which has structure Ih at 40 K, is ≈921 kg m⁻³ (Desch et al. 2009), much lower than Haumea's mean density. Haumea, therefore, is certainly differentiated. A basic question, then, is whether a differentiated Haumea can be consistent with axes that match a Jacobi ellipsoid.

To answer this question, we have calculated the gravitational potential of a differentiated Haumea, as described in Probst (2015). We model Haumea as two nested, aligned triaxial ellipsoids. We assume that Haumea's outer surface is an ellipsoid with axes a = 960 km, b = 770 km, and c = 495 km, and we allow its core to have arbitrary density ρ_core and arbitrary axis ratios p_c = b_c/a_c and q_c = c_c/a_c. For a given density ρ_core and axis ratios p_c and q_c, the core axis a_c is chosen so that the mass of the core plus the mass of the ice mantle, with density ρ_ice = 921 kg m⁻³, equal the mass of Haumea. We then calculate whether the surface and the core–mantle boundary (CMB) are equipotential surfaces.

An equilibrium solution must have the equipotential surfaces coincident with the surface and CMB, or else vortical flows will be generated. In the absence of external and viscous forces and large internal flows, the vorticity $\overrightarrow{\omega }$ follows the equation

$$\begin{eqnarray}&&\displaystyle \frac{D\overrightarrow{\omega }}{{Dt}}=\displaystyle \frac{1}{{\rho }^{2}}\,\overrightarrow{{\rm{\nabla }}}\rho \times \overrightarrow{{\rm{\nabla }}}P,\end{eqnarray} \tag{ 4 }$$

where P is the pressure and ρ is the density. If the gradients of ρ and P are misaligned by an angle Θ, then vorticity will be generated at a rate $\sim (P/\rho )/{R}^{2}(\sin {\rm{\Theta }})$, where R is comparable to the mean radius. In a time τ, the vortical flows will circulate at rates comparable to the rotation rate, ω_H, where $\tau \,\sim {\omega }_{{\rm{H}}}\rho {R}^{2}/P{\left(\sin {\rm{\Theta }}\right)}^{-1}$. Assuming ω_H = 4.46 × 10⁻⁴ s⁻¹, ρ = 921 kg m⁻³, P = 18 MPa, and R ∼ 725 km, the timescale $\tau \sim 4/\sin {\rm{\Theta }}$ yr. Even a small mismatch, with Θ ∼ 1°, would lead to significant vortical flows within hundreds of years. Equilibrium solutions demand Θ = 0°, i.e., that $\overrightarrow{{\rm{\nabla }}}\rho$ and $\overrightarrow{{\rm{\nabla }}}P$ are parallel. In hydrostatic equilibrium, the net force is

$$\begin{eqnarray}&&\overrightarrow{F}=-\overrightarrow{{\rm{\nabla }}}P+\rho {\overrightarrow{g}}_{\mathrm{eff}},\end{eqnarray} \tag{ 5 }$$

where ${\overrightarrow{g}}_{\mathrm{eff}}$ measures the acceleration due to gravity, as well as centrifugal support due to Haumea's rotation. If the net force is zero, then it must be the case that $\overrightarrow{{\rm{\nabla }}}\rho$ is parallel to ${\overrightarrow{g}}_{\mathrm{eff}}$. In other words, the gradient in the effective gravitational potential must be parallel to the gradient in density, and surfaces of density discontinuity must be equipotential surfaces.

We solve for the effective gravitational potential by discretizing an octant of Haumea on a Cartesian grid with 60 evenly spaced zones along each axis. If a rectangular zone is entirely inside the triaxial ellipsoid defined by the core, its density is set to ρ_core; if it is entirely outside the triaxial ellipsoid defined by the surface, its density is set to zero; and if it is entirely between these two ellipsoids, its density is set to ρ_ice. For zones straddling the core and ice mantle, or straddling the ice mantle and the exterior, the density is found using a Monte Carlo method. An array of about 100 points on the surface is then generated, by creating a grid of N_θ ≈ 10 angles θ from 0 to π rad and of N_ϕ ≈ 10 angles ϕ from 0 to 2π rad. The points on the surface are defined by ${x}^{* }=a\,\sin \theta \,\cos \phi$, ${y}^{* }=b\,\sin \theta \,\sin \phi$, ${z}^{* }=c\,\cos \theta$. At each point we find the vector normal to the surface, $\overrightarrow{n}\,=(2{x}^{* }/{a}^{2}){\hat{e}}_{x}+(2{y}^{* }/{b}^{2}){\hat{e}}_{y}\,+(2{z}^{* }/{c}^{2}){\hat{e}}_{z}$, as well as the gravitational acceleration $\overrightarrow{g}$, found by summing the gravitational acceleration vectors from each zone's contribution. To this we add an additional contribution due to centrifugal support:

$$\begin{eqnarray}&&{\overrightarrow{g}}_{\mathrm{eff}}=\overrightarrow{g}+{\omega }^{2}{x}^{* }\,{\hat{e}}_{x}+{\omega }^{2}{y}^{* }\,{\hat{e}}_{y}.\end{eqnarray} \tag{ 6 }$$

Once these vectors are found at each of the surface points defined by θ and ϕ, we find the following quantity by summing over all points:

$$\begin{eqnarray}&&{ \mathcal M }=\displaystyle \frac{1}{{N}_{\theta }{N}_{\phi }}\,\displaystyle \sum _{{N}_{\theta }}\displaystyle \sum _{{N}_{\phi }}\,\displaystyle \frac{\overrightarrow{n}\cdot {\overrightarrow{g}}_{\mathrm{eff}}}{\left|\overrightarrow{n}\right|\left|{\overrightarrow{g}}_{\mathrm{eff}}\right|}.\end{eqnarray} \tag{ 7 }$$

We also define "fit angle," the mean angular deviation between the surface normal and the effective gravitational field:

$$\begin{eqnarray}&&{\rm{\Theta }}={\cos }^{-1}\,{ \mathcal M }.\end{eqnarray} \tag{ 8 }$$

If the equipotential surfaces are coincident with the surface, then ${ \mathcal M }=1$ and Θ = 0°. In a similar fashion, we define an identical metric ${ \mathcal M }$ on the CMB as well.

In Figure 1, we plot the fit angle on the surface and CMB of Haumea, as a function of the assumed core axis ratios p_c and q_c. A core density of 2700 kg m⁻³ is assumed. The core axis ratios that minimize Θ on the outer surface are p_c = b_c/a_c ≈ 0.80 and q_c = c_c/a_c ≈ 0.51. For this combination, Θ ≈ 1°. Meanwhile, the core axis ratios that minimize Θ on the CMB are p_c ≈ 0.82 and q_c ≈ 0.55. For this combination, Θ ≈ 2°. Significantly, the parameters that minimize Θ on the outer surface are not those that minimize Θ on the CMB. Given the change in Θ with changes in p_c and q_c, Θ on either the surface or the CMB must be at least a few degrees. Figure 1 shows numerous patches of higher-than-expected angles. These patches are caused by the random nature in which the grid cells straddling the CMB are populated. The coarseness of the grid leads to the code creating a bumpiness in the CMB surface, which then upwardly skews the calculation of the average fit angle in places. This is because where the surface is bumpy, the surface normal vector can be significantly different from the gravitational acceleration vector. The patches become more numerous, but smaller in magnitude, with increasing numerical resolution.

Zoom In Zoom Out Reset image size

We have repeated the analysis for other core densities of 3000 and 3300 kg m⁻³. In those cases the discrepancy between what parameters p_c and q_c minimize Θ on the surface versus what parameters minimize it on the CMB grows even larger. If Haumea has axes a = 960 km, b = 770 km, and c = 495 km and is divided into a rocky core and an icy mantle, the only way to maintain hydrostatic equilibrium on the surface and CMB is for the core density to be as close as possible to the inferred bulk density of Haumea, ≈2600 kg m⁻³, and for the core and surface axis ratios to converge. The solution is naturally driven to one of uniform density. In this case, the core composes over 96% of the mass of Haumea, and the ice thickness is <10 km on the a- and b-axes and <5 km on the c-axis. Even for this case, though, the effective equipotential surfaces fail to coincide with the surface or CMB, by several degrees. This suggests that the only way for Haumea to have axes consistent with a Jacobi ellipsoid is for it to have essentially no ice mantle, less than a few kilometers thick.

These investigations reveal two facts. First, a Haumea divided into a rocky core and icy mantle cannot have axes equal to those of a Jacobi ellipsoid. This lends some support to the finding by Ortiz et al. (2017) that Haumea's axes deviate significantly from a Jacobi ellipsoid's. Second, if Haumea cannot conform to a Jacobi ellipsoid, then it is not possible to use analytical formulae to describe its shape, and a more powerful technique must be used to derive its internal structure.

To calculate the shape of a rapidly rotating Haumea with a rocky core and icy mantle, we have written a code, named kyushu,³ that calculates the internal structure and figure of equilibrium of a differentiated body undergoing rapid, uniform rotation. Our algorithm is adapted from that of Hachisu (1986a, 1986b; hereafter H86a, H86b), for calculating the structure of stars orbiting each other in binary systems, as follows.

The Hachisu (1986a, 1986b) algorithm relies on using a governing equation derived from the Bernoulli equation, at each location on a three-dimensional grid:

$$\begin{eqnarray}&&\int {\rho }^{-1}\,{dP}+{\rm{\Phi }}-\int {{\rm{\Omega }}}^{2}{r}_{\perp }\,{{dr}}_{\perp }=C,\end{eqnarray} \tag{ 9 }$$

where the first term is the enthalpy, H, the second term is the gravitational potential energy, the third term is the rotational energy, and C is a constant. Here r_⊥ is the distance from the rotation axis.

The grid is defined in spherical polar coordinates, the variables being distance from the origin, r, the cosine of the polar angle, μ, and the azimuthal angle, ϕ. A discretized grid of r_i, with i = 1, 2, ... N_r, is defined, with r uniformly spaced between r₁ = 0 and r_Nr = R. Likewise, a discretized grid of μ_j, with j = 1, 2, ... N_μ, is defined, with μ uniformly spaced between μ₁ = 0 and μ_Nϕ = 1, and a discretized grid of ϕ_k, with k = 1, 2, ... N_ϕ, is defined, with ϕ uniformly spaced between ϕ₁ = 0 and ϕ_Nϕ = π/2. Symmetries across the equatorial plane and n = 2 symmetry about the polar axis are assumed. Quantities in the above equation, including density ρ_ijk, gravitational potential Φ_ijk, etc., are defined on the intersections of grid lines. Typical values in our calculation are R = 1300 km, N_r = 391, N_μ = 33, and N_ϕ = 33, meaning that quantities are calculated at 33 × 33 × 391 = 425,799 locations.

The enthalpy term can be calculated if the density structure and equation of state are provided. For example, if P = Kρ^γ (as for an adiabatic gas), then H = (γ)/(γ − 1)P/ρ and is immediately known as a function of the local pressure and density. For planetary materials, such as olivine, clays, or water ice, it would be more appropriate to use a Vinet (Vinet et al. 1987) or Birch-Murnaghan (Birch 1947) equation of state, with the bulk modulus and the pressure derivative of the bulk modulus specified. This equation of state can then be integrated to yield the enthalpy. The recent paper by Price & Rogers (2019) provides formulas for this. For our purposes, we neglect the self-compression of the planetary materials inside Haumea. The bulk moduli of planetary materials is typically tens of GPa, while the maximum pressure inside Haumea is < 0.4 GPa, so self-compression can be ignored. For ease of calculation, we, therefore, assume uniform densities in the ice mantle and in the rocky core and compute the enthalpy accordingly.

The gravitational potential term is found by numerical integration of an expansion of the gravitational potential in spherical harmonics, using Equations (2), (3), and (33)–(36) of H86b, using n = 2 symmetry and typically N_l = 16 terms in the expansion. The rotational term is defined to be Ω²Ψ, where ${\rm{\Psi }}=-{r}_{\perp }^{2}/2$, r_⊥ again being the distance from the axis. The terms H, Φ, and Ω² Ψ all spatially vary, but their sum is a constant C at all locations. The Hachisu algorithm exploits this fact by fixing two spatial points "A" and "B" to be on the boundary of the body. Point "A" lies at r = r_A, μ = 0 (in the equatorial plane), ϕ = 0, or $x=r\,{(1-{\mu }^{2})}^{1/2}\,\cos \phi ={r}_{A}$, $y=r\,{(1-{\mu }^{2})}^{1/2}\,\sin \phi =0$, z = r μ = 0. Point "B" lies at r = r_B, μ = 0 (in the equatorial plane), ϕ = π/2, or x = 0, y = r_B, z = 0. For a triaxial ellipsoid (x/a)² + (y/b)² +(z/c)² = 1, points A and B refer to the long and intermediate axes of the body in the equatorial plane, with r_A = a and r_B = b. At these two locations, H = 0, and the Hachisu algorithm then solves for the only two values of C and Ω that allow H = 0 at both these boundary points. With C and Ω defined, H is found at all locations, and the enthalpy integral $H=\int {\rho }^{-1}\,{dP}$ is inverted to find the density ρ at each location. Locations with H < 0 are assigned zero density. The farthest point with nonzero density along the z-axis can be equated with c of a triaxial ellipsoid, although of course the shape need not necessarily be a triaxial ellipsoid. After adjusting the density everywhere, the code then recalculates the gravitational potential and performs the same integrations, solving iteratively until the values of Ω and the densities ρ at all locations converge. Because the densities and the volume of the body are changed with each iteration, the mass of the object is a varying output of the model.

We apply the H86a and H86b algorithms as part of a larger iterative procedure that introduces two new variables to an equation of state: P_cmb and ρ_core. We assume that at locations within the body with pressures P < P_cmb, ρ = ρ_ice ≡ 921 kg m⁻³. At higher pressures P > P_cmb, we assume ρ = ρ_core. This divides the body into a core and an icy mantle, each with distinct densities, with the pressure equal to a uniform value P_cmb everywhere on the CMB. For P > P_cmb, H = P_cmb/ρ_ice + (P − P_cmb)/ρ_core. This form of the equation of state ignores self-compression, as is appropriate in bodies of Haumea's size made of materials like water ice, olivine, or hydrated silicates. The bulk modulus of water ice is 9.2 GPa (Shaw 1986), far higher than the likely pressures in the ice shell, <20 MPa (Section 4). Likewise, the bulk modulus of olivine is 126 GPa (Núñez-Valdez et al. 2013), and that of the hydrated silicate clay antigorite is 65 GPa (Capitani & Stixrude 2012), which is far higher than the maximum pressure inside Haumea (<300 MPa). We, therefore, expect self-compression to change the densities by <1%, and we are justified in assuming uniform densities. This equation of state allows a very simple inversion to find the density: for H < P_cmb/ρ_ice, the density is simply ρ_ice, and for H > P_cmb/ρ_ice, ρ = ρ_core.

With this definition, we iterate as follows to find P_cmb and ρ_core. In each application of the Hachisu algorithms, we initialize with a density distribution with ρ = 0 for (x/a)² +(y/b)² + (z/c)² > 1, ρ = ρ_ice for (x/a)² + (y/b)² + (z/c)² < 1, but ρ = ρ_core for (x/a)² + (y/b)² + (z/c)² < ξ². That is, we define Haumea's surface to be a triaxial ellipsoid, and we define its core to be a similar triaxial ellipsoid with aligned axes, smaller in size by a factor of ξ. The value of ξ is chosen so that the total mass of the configuration equals M_H, the mass of Haumea:

$$\begin{eqnarray}&&\xi ={\left[\displaystyle \frac{3{M}_{{\rm{H}}}/(4\pi {abc})-{\rho }_{\mathrm{ice}}}{{\rho }_{\mathrm{rock}}-{\rho }_{\mathrm{ice}}}\right]}^{1/3}.\end{eqnarray} \tag{ 10 }$$

The value of c is unknown and an output of the code, so we initialize our configuration with c = b. We apply the Hachisu algorithms several times. First, we define P_cmb ≈ 0 MPa and define ρ_core = 2700 kg m⁻³, apply the Hachisu algorithms, and calculate the mass of the body, M. Because this will not match Haumea's mass, M_H, we multiply ρ_core by a factor M_H/M and reapply the Hachisu algorithms. We do this until we have found an equilibrium configuration with Haumea's mass. Outputs of the code include ${\rho }_{\mathrm{core}}$ and Ω, and in general Ω will be smaller than Haumea's true angular frequency Ω_H = 2π/P_rot for P_cmb = 0 MPa (P_rot is the rotational period). We then repeat the procedure, finding an equilibrium configuration with Haumea's mass, having P_cmb = 40 MPa. An output of the code will be a different ρ_core and a different Ω, which in general will be >Ω_H. If the solution is bracketed, we use standard bisection techniques to find P_cmb that yields an equilibrium configuration consistent not just with Haumea's mass, but with its period as well.

Thus, inputs of the code include r_A = a and r_B = b, as well as Haumea's mass M_H and angular velocity Ω_H. The outputs of the code include the density ρ at all locations, the calculated mass M (which should comply with M = M_H), the angular frequency Ω (which should equal 2π/P_rot), and the values of c, P_cmb, and ρ_core.

To benchmark the kyushu code, we ensure that it reproduces a Jacobi ellipsoid when a homogeneous density is assumed. A Jacobi ellipsoid with Haumea's mass and rotation period and an axis a = 960 km would have axes b = 774.2 km and c = 498.8 km and a uniform density of 2579.7 kg m⁻³; this is similar to, but does not exactly equal, the solution favored by Lockwood et al. (2014), who fit Haumea's light curve assuming a = 960 km, b = 770 km, c = 495 km, and a uniform density of 2600 kg m⁻³. Running kyushu and assuming axes of a = 960 km and b = 774.2 km (b/a =0.806), the code finds an acceptable solution after about seven bracketing iterations. The mass matches Haumea's mass to within 0.03% and the rotation period to within 0.04%. The solution found is one with a very low value of P_cmb = 0.42 MPa, so that the body is essentially uniform in density, with an ice layer <5 km in thickness (the resolution of the code). The interior of the body has a uniform density of 2580.4 kg m⁻³, and the short axis has length c = 499.5 km (c/a = 0.520). The density and c-axis match a uniform-density Jacobi ellipsoid to within 0.03% and 0.14%, respectively. The code is, therefore, capable of finding the analytical solution of a uniform-density Jacobi ellipsoid, if the imposed a- and b-axes are consistent with such a solution.

We have run the kyushu code for 30 different combinations of a- and b-axes, with a varying from 950 to 1075 km in increments of ∼25 km and b varying from 800 to 900 km in increments of ∼25 km. In comparing a subset of these runs with runs performed with 10 km increments, we found that convergence of key outputs (axes, densities) is 0.3% or less, so we consider 25 km to be numerically converged. We find families of solutions that can conform to Haumea's mass (M = 4.006 × 10²¹ kg) and rotation period (P_rot = 3.9155 hr). Output quantities include the core density, the average (bulk) density, the outer c-axis, the shape of the CMB, and the thickness of the ice layer above the core.

In Figure 2, we plot the following quantities as functions of imposed a- and b-axes: the average (bulk) density ρ_avg, the (semi)axis c, and the thickness of the ice layer along the a-, b-, and c-axes. The density of the ice layer was imposed to be 921 kg m⁻³. In Figure 3, we plot as functions of a and b the following: the core density, ρ_core, the pressure at the CMB, P_cmb, and the semiaxes a_c, b_c, and c_c of the core. Entire families of solutions are found across the range of a and b that we explored, with the exception of simultaneous combinations of large a and large b. The input a × b combinations of 1025 × 900, 1050 × 900, 1075 × 900, and 1075 × 800 did not yield solutions because, when initializing with these parameters, kyushu was not able converge at Haumea's mass and rotation rate. Simultaneously imposing large a and b yields a low bulk density ρ_avg and large parameter ω²/(πGρ_avg) that cannot yield a c-axis consistent with Haumea's mass.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Across the explored parameter space that yielded solutions, the average density of Haumea ranges from 1905 to 2495 kg m⁻³. As might be expected, the average (bulk) density of Haumea is equally sensitive to both a and b, being inversely proportional to the volume and, therefore, to the product ab. As an example solution, we take a = 1050 km and b = 840 km, for which the average density is 2018 kg m⁻³.

For allowed solutions, the shortest (semi)axis c ranges from 504 to 546 km. The c-axis is equally sensitive to both a and b, being large when the product ab is large. For the case with a = 1050 km and b = 840 km, we find c = 537 km.

Likewise, the thickness of the ice layer increases with increasing a and b (and, therefore, c). The ice is always thickest along the a-axis, ranging from 15 to 210 km across the explored range; it is intermediate in thickness along the b-axis, ranging from 15 to 150 km; and it is thinnest along the c-axis, ranging from 5 to 80 km. For the case with a = 1050 km and b = 840 km, we find ice thicknesses of 167, 117, and 67 km along the a-, b-, and c-axes, respectively.

Across the explored parameter space that yielded solutions, the density of Haumea's core ranges from 2560 to 2740 kg m⁻³. Core density is much more sensitive to the a-axis than the b-axis, and it tends to be greater when a is greater. For the case with a = 1050 km and b = 840 km, we find a core density of 2680 kg m⁻³.

Across the explored range, the pressure at the CMB ranges from 3.2 to 36.6 MPa, with the lowest values corresponding to the thinnest ice layers and smallest values of a and b. For the case with a = 1050 km and b = 840 km, P_cmb = 30.4 MPa. Finally, the size of the core ranges considerably. In general, the core approximates a triaxial ellipsoid, with the longest axis parallel to the a-axis. For the extreme case with a = 950 km, b = 800 km, and c = 504 km (a small Haumea), we find core axes of a_c = 935 km, b_c = 785 km, and c_c = 499 km. For this configuration, b_c/a_c = 0.840 and c_c/a_c = 0.534, very close to the axis ratios for the surface, b/a = 0.842 and c/a = 0.531. The mean size of the core, relative to the mean size of the surface, is ξ = 0.991. That is, the ice layer thickness is only ∼1% of the radius of Haumea and composes 1.5% of Haumea's volume. For the opposite extreme case of a large Haumea, with a = 1050 km, b = 875 km, and c = 546 km, we find core axes of a_c = 840 km, b_c = 725 km, and c_c = 466 km. For this configuration, b_c/a_c = 0.86 and c_c/a_c = 0.55, close to the axis ratios for the surface, b/a = 0.83 and c/a = 0.52. The mean size of the core, relative to the mean size of the surface, is ξ = 0.888, meaning that the ice layer thickness is 11% of the radius of Haumea and composes 22% of its volume. For the case we consider typical, with a = 1050 km and b = 840 km, the core semiaxes are ≈883 × 723 × 470 km, which yields b_c/a_c = 0.82 and c_c/a_c = 0.53. The mean size of the core, relative to the mean size of the surface, is ξ = 0.909, and the ice layer composes 17.2% of Haumea's volume.

In general, for larger assumed sizes of Haumea, the core becomes somewhat denser as Haumea itself becomes lower in density. The core takes up a smaller fraction of the volume of Haumea. While the core remains roughly similar in shape to the ellipsoid defined by the surface, there is a tendency for the core to become slightly more spherical as Haumea's assumed size increases.

To quantitatively test whether Haumea's core and surface are both triaxial ellipsoids in our typical case, we calculated the maximum deviation from 1 of (x/a_c)² + (y/b_c)² + (z/c_c)², where x, y, and z are computed for each angle combination θ and ϕ and radius defined by $r=(r({{ir}}_{\mathrm{core}}-1)+r({{ir}}_{\mathrm{core}}))/2$. Here ir_core is the index of the first radial zone outside the core at that θ and ϕ. We find that a triaxial ellipsoid shape is consistent with the core to within 1.5% and the surface to within 0.5%. Both are within the code's resolution error. It is remarkable that the CMB solution is driven to the shape of a triaxial ellipsoid. This justifies the assumption in Section 2 that the surface and core would be triaxial ellipsoids if Haumea is differentiated.

5.1. Reconciling the Light Curve and Occultation Data Sets

5.2. Aqueous Alteration of Haumea's Core and Its Astrobiological Potential

A large range of axes a (from 950 to 1075 km) and b (from 800 to 900 km) are consistent with bodies with Haumea's mass and rotation rate and are fluid configurations in hydrostatic equilibrium. Our favored solution with a = 1050 km and b = 840 km has a mass fraction of ice of 17.2%, but this value could range from 1% to 22% across the range we explored. Across this range, however, the allowable core density varies only slightly, from ρ_core = 2560 to 2740 kg m⁻³, deviating by only a few percent from our favored density of ρ_core =2680 kg m⁻³. This is very close to the average density previously inferred for Haumea, but this appears to be coincidental. A robust result of our analysis is that Haumea's core has a density of ≈2600 kg m⁻³, overlaid by an ice mantle.

Comparison of the density of Haumea's core to other planetary materials provides strong clues to Haumea's history. Grain densities of ordinary and enstatite chondrites are typically >3600 kg m⁻³, and their bulk densities are typically ≈3300 kg m⁻³ because of ∼10% porosity (Wilkison et al. 2003; Consolmagno et al. 2008). Carbonaceous chondrites are marked by lower grain densities, on average 3400 kg m⁻³ (ranging from 2400 to 5700 kg m⁻³ depending on the type of chondrite), higher porosities of ≈15%–35%, and bulk densities closer to 2000 kg m⁻³ (Consolmagno et al. 2008; Macke et al. 2011). The difference is that carbonaceous chondrites are largely composed of products of aqueous alteration. In fact, the more oxidized groups of carbonaceous chondrites have higher porosity (Macke et al. 2011). Hydrated silicates typically have densities in this range. Clays such as montmorillonite, kaolinite, illite, and mica typically have densities between 2600 and 2940 kg m⁻³ (Osipov 2012). This strongly suggests that Haumea's core is composed of hydrated silicates and that Haumea's core was aqueously altered in its past.

Serpentinization is the process by which silicates typical of the dust in the solar nebula react with water on an asteroid or planet, producing new phyllosilicate minerals. An archetypal reaction would be

$$\begin{eqnarray*}&&\begin{array}{l}1.000\,{\left({\mathrm{Mg}}_{0.71}{\mathrm{Fe}}_{0.29}\right)}_{2}{\mathrm{SiO}}_{4}+1.142\,{{\rm{H}}}_{2}{\rm{O}}\,\to \\ \quad 0.474\,{\mathrm{Mg}}_{3}{\mathrm{Si}}_{2}{{\rm{O}}}_{5}{\left(\mathrm{OH}\right)}_{4}+0.193\,{\mathrm{Fe}}_{3}{{\rm{O}}}_{4}\\ \qquad +\ 0.052\,{\mathrm{SiO}}_{2}+0.194\,{{\rm{H}}}_{2}.\end{array}\end{eqnarray*}$$

In this reaction, 1 kg of olivine with fayalite content typical of carbonaceous chondrites may react with 0.129 kg of water to produce 0.826 kg of chrysotile, 0.281 kg of magnetite, 0.020 kg of silica, and 0.002 kg of hydrogen gas, which escapes the system. The total density of olivine (density = 3589 kg m⁻³) plus ice (density = 921 kg m⁻³) before the reaction is 2697 kg m⁻³. After the reaction, the mixture of chrysotile (density = 2503 kg m⁻³) plus magnetite (density = 5150 kg m⁻³) plus silica (density = 2620 kg m⁻³) has a total density of 2874 kg m⁻³ (Coleman 1971). Including a 10% porosity typical of carbonaceous chondrites, the density of the aqueously altered system would be 2612 kg m⁻³, remarkably close to the inferred density of Haumea's core.

If Haumea's core underwent significant aqueous alteration, some of this material may have dissolved in the water and ultimately found its way into the ice mantle of Haumea. In fits to Haumea's reflectance spectrum, Pinilla-Alonso et al. (2009) found that the most probable surface composition was an intimate mixture of half crystalline and half amorphous water ice, with other components composing <10% of the surface, but similar modeling by Trujillo et al. (2007) found that Haumea's surface is best fit by a mixture of roughly 81% crystalline water ice and 19% kaolinite. Kaolinite was added to the fit to provide a spectrally neutral but bluish absorber; few other planetary materials contribute to the reflectance spectrum in this way. Kaolinite is a common clay mineral [Al₂Si₂O₅(OH)₄] very similar in structure to chrysotile, produced by weathering of aluminum silicate minerals like feldspars.

A variety of phyllosilicates have been observed by the Dawn mission on the surface of Ceres (Ammannito et al. 2016), strongly suggesting aqueous alteration of silicates within a porous, permeable core or a convecting mudball (Bland et al. 2006; Travis 2017). If kaolinite can be confirmed in Haumea's mantle, this would provide strong support for the aqueous alteration of Haumea's core.

Preliminary modeling by Desch & Neveu (2015) suggests that aqueous alteration of Haumea's core is a very likely outcome. Haumea, or its pre-collision progenitor, could have been differentiated into a rocky olivine core and icy mantle. Desch & Neveu (2015) show that many factors can lead to cracking of a rocky core on small bodies. Microcracking by thermal expansion mismatch of mineral grains or by thermal expansion of pore water during heating (as the core heats by radioactive decay over the first <1 Gyr) will almost certainly introduce microfractures. These would be widened by chemical reactions and water pressurization, etc., leading to macrofractures. Cracks can heal by ductile flow of rock, but the rate is highly sensitive to temperature; below about 750 K, healing of cracks takes longer than the age of the solar system. Therefore, it is highly likely that hydrothermal circulation of water through a cracked core would ensue. Thermal modeling by Desch & Neveu (2015) suggests that Haumea's interior could be effectively fully convective, allowing water and olivine to fully react and produce phyllosilicates. Circulation of water also would help cool the core, preventing temperatures from exceeding 750 K, ensuring that fractures remain open and that the hydrated silicates would not dehydrate. Liquid water is predicted to have existed for ∼10⁸ yr, although further geochemical modeling is needed to test more proposed scenarios for Haumea's structure and evolution.

A long (∼10⁸ yr) duration of aqueous alteration suggests a period of habitability within Haumea. To develop and survive, life as we know it requires water and a long-lasting environment with little temperature variability (Davis & McKay 1996). With central temperatures approaching 750 K and surface temperatures near 40 K, a large fraction of Haumea's interior would have had intermediate temperatures consistent with liquid water (Castillo-Rogez & Lunine 2012). The origin of life is also thought to require a substrate to protect and localize biochemical reactions. Clays such as montmorillonite can act as this substrate because they can bind substantial water, are soft, and delaminate easily. Clays can also promote the assembly of RNA from nucleosides and can stimulate micelles to form vesicles (Travis 2017). The interior of Haumea may have at one point resembled regions beneath the seafloor experiencing hydrothermal circulation. These regions are conducive to life: Czaja et al. (2016) discovered archaea anaerobically metabolizing H₂S in such environments, and other studies have confirmed that microbes exist deep in fractures of hot environments (Jannasch & Mottl 1985).

5.3. Implications for the Mass of the Collisional Family

This paper presents numerical modeling designed to test three questions about the KBO Haumea: (1) Is Haumea a Jacobi ellipsoid? If it is differentiated, what is Haumea's shape? (2) Is Haumea a fluid in hydrostatic equilibrium? (3) Can Haumea's occultation and light-curve data be reconciled? We aimed to address these questions with the goal of understanding the composition and structure of Haumea to learn about its collisional history and evolution.

We have written a code, kyushu, based on the algorithms of Hachisu (H86a, H86b) to calculate the internal structure of a rapidly rotating differentiated body based on input parameters such as the semiaxes a and b. Although we did not explore all parameter space, Haumea appears to be best approximated as a differentiated triaxial ellipsoid body in hydrostatic equilibrium with axes a = 1050 km, b = 840 km, and c = 537 km. This shape fits the Ortiz et al. (2017) occultation shadow and is close to the light-curve data. As this shape, Haumea has core axes a_c = 883 km, b_c = 723 km, and c_c = 470 km, ρ_avg =2018 kg m⁻³, and ρ_core = 2680 kg m⁻³, which equates to an ice mantle composing ∼17% of Haumea's mass and ranging from 71 to 170 km in thickness. Haumea's albedo is p_v ∼ 0.66 in this case.

In contrast to previous studies (Rabinowitz et al. 2006; Lockwood et al. 2014), our results suggest that Haumea's ice crust amounts to a significant portion of the body. Due to the thickness of the ice, Haumea's core has a relatively high density, indicating that the composition of the core is a hydrated silicate (perhaps kaolinite). For the core to be hydrated, a long period (∼10⁸ yr) of serpentinization must have occurred during which regions of the core were potentially habitable. The thick ice crust also suggests that Haumea's collisional family (icy objects a few percent of the mass of Haumea) was produced from only a small portion of the ice Haumea started with, before Haumea suffered the collision. Insights into this type of mantle-stripping collision could be applicable to modeling metal-rich, fast-rotating triaxial ellipsoid 16 Psyche, the focus of the upcoming NASA Psyche mission (Elkins-Tanton et al. 2016).

As this study continues, we would like to expand parameter space to obtain more precise results. We can explore how Haumea would change shape or composition if we were to use different ice densities, porosities, and angles/orientations to better match the shadow, in addition to matching the light-curve amplitude/phase more precisely and using an appropriate equation of state to include the compressibility of materials. Haumea is a unique and interesting body worthy of study for its own sake, but understanding Haumea can provide insights into fundamental processes such as subsurface oceans/aqueous alteration on small bodies and dynamics of mantle-stripping collisions, acting across the solar system.

We thank Darin Ragozzine and Sarah Sonnett for helpful discussions about the collisional family and Haumea's light curve, Steve Schwartz and Viranga Perera for useful discussions about how to model Haumea using smoothed particle hydrodynamics codes, Anand Thirumualai for assisting with some calculations, Marc Neveu for conversations on Haumea's internal thermal evolution, and Leslie Rogers and Ellen Price for introducing us to the Hachisu (H86a, H86b) algorithm and for general discussions about how they implemented the Hachisu algorithm for exoplanets. Comments by an anonymous reviewer improved the clarity of the paper. We gratefully acknowledge partial support of the NASA Solar Systems Workings Program under grant 80NSSC19K0028.

Here we derive the formulae needed to calculate the axes of Haumea's shadow as it occults a star. We assume that Haumea's surface is a triaxial ellipsoid with long axis along the x-direction, with axes a > b > c, defined by those points that satisfy

$$\begin{eqnarray*}&&f(x,y,z)=\displaystyle \frac{{x}^{2}}{{a}^{2}}+\displaystyle \frac{{y}^{2}}{{b}^{2}}+\displaystyle \frac{{z}^{2}}{{c}^{2}}=1.\end{eqnarray*}$$

We assume that the star lies in a direction

$$\begin{eqnarray*}&&{\hat{e}}_{\mathrm{LOS}}=\cos \psi \,\sin \phi \,{\hat{e}}_{x}+\sin \psi \,\sin \phi \,{\hat{e}}_{y}+\cos \phi \,{\hat{e}}_{z}.\end{eqnarray*}$$

Here ϕ is the angle between the line of sight (from us through Haumea to the star) and Haumea's pole (along the z-axis). We can define two unit vectors in the plane of the sky:

$$\begin{eqnarray*}&&{\hat{e}}_{1}=-\sin \psi \,{\hat{e}}_{x}+\cos \psi \,{\hat{e}}_{y}\end{eqnarray*}$$

and

$$\begin{eqnarray*}\begin{array}{rcl}{\hat{e}}_{2} & = & {\hat{e}}_{1}\times {\hat{e}}_{\mathrm{LOS}}=+\cos \psi \cos \phi \,{\hat{e}}_{x}\\ & & +\,\sin \psi \cos \phi \,{\hat{e}}_{y}-\sin \phi \,{\hat{e}}_{z}.\end{array}\end{eqnarray*}$$

Haumea's limb is the locus of those points, defined by ${\boldsymbol{r}}$, such that the line of sight is tangent to the surface, or perpendicular to the normal, so that

$$\begin{eqnarray*}&&{\rm{\nabla }}f\,\cdot \,{\hat{e}}_{\mathrm{LOS}}=0.\end{eqnarray*}$$

All of these points satisfy

$$\begin{eqnarray*}&&z=-{c}^{2}\,\tan \phi \left(\displaystyle \frac{x\,\cos \psi }{{a}^{2}}+\displaystyle \frac{y\,\sin \psi }{{b}^{2}}\right),\end{eqnarray*}$$

which define a plane inclined to the sky. The intersection of the plane with the ellipsoid defines an ellipse, and the projection of this ellipse onto the plane of the sky—Haumea's shadow—also is an ellipse.

We project the points on Haumea's limb onto the plane of the sky by recasting ${\boldsymbol{r}}$ in the coordinate system using ${\hat{e}}_{1}$, ${\hat{e}}_{2}$, and ${\hat{e}}_{\mathrm{LOS}}$:

$$\begin{eqnarray*}\begin{array}{rcl}{\boldsymbol{r}} & = & \left({\boldsymbol{r}}\cdot {\hat{e}}_{1}\right)\,{\hat{e}}_{1}+\left({\boldsymbol{r}}\cdot {\hat{e}}_{2}\right)\,{\hat{e}}_{2}+\left({\boldsymbol{r}}\cdot {\hat{e}}_{\mathrm{LOS}}\right)\,{\hat{e}}_{\mathrm{LOS}}\\ & = & s\,{\hat{e}}_{1}+t\,{\hat{e}}_{2}+u\,{\hat{e}}_{\mathrm{LOS}},\end{array}\end{eqnarray*}$$

with

$$\begin{eqnarray*}&&s={\boldsymbol{r}}\cdot {\hat{e}}_{1}=-x\,\sin \psi +y\,\cos \psi \end{eqnarray*}$$

and

$$\begin{eqnarray*}&&t={\boldsymbol{r}}\cdot {\hat{e}}_{2}=x\,\cos \psi \,\cos \phi +y\,\sin \psi \,\cos \pi -z\,\sin \phi .\end{eqnarray*}$$

All the points on the limb have z related to x and y as above, so the boundary of the shadow, which equals the projection of the limb onto the plane of the sky, is defined by

$$\begin{eqnarray*}&&s=-x\,\sin \psi +y\,\cos \psi \end{eqnarray*}$$

and

$$\begin{eqnarray*}&&\displaystyle \frac{t}{\cos \phi }=x\,\cos \psi \left(1+\displaystyle \frac{{c}^{2}}{{a}^{2}}\,{\tan }^{2}\phi \right)+y\,\sin \psi \left(1+\displaystyle \frac{{c}^{2}}{{b}^{2}}\,{\tan }^{2}\phi \right).\end{eqnarray*}$$

Inverting, we find x, y, and z in terms of s and t for points along the limb:

$$\begin{eqnarray*}&&x=\displaystyle \frac{1}{{\rm{\Delta }}}\left[-\sin \psi \left(1+\displaystyle \frac{{c}^{2}}{{b}^{2}}\,{\tan }^{2}\phi \right)\,s+\cos \psi \,\displaystyle \frac{t}{\cos \phi }\right],\end{eqnarray*}$$

$$\begin{eqnarray*}&&y=\displaystyle \frac{1}{{\rm{\Delta }}}\left[+\cos \psi \left(1+\displaystyle \frac{{c}^{2}}{{a}^{2}}\,{\tan }^{2}\phi \right)\,s+\sin \psi \,\displaystyle \frac{t}{\cos \phi }\right],\end{eqnarray*}$$

and

$$\begin{eqnarray*}\begin{array}{rcl}z & = & -\displaystyle \frac{\tan \phi }{{\rm{\Delta }}}\left[\sin \psi \,\cos \psi \left(\displaystyle \frac{{c}^{2}}{{b}^{2}}-\displaystyle \frac{{c}^{2}}{{a}^{2}}\right)\right.\\ & & +\,\left.\left(\displaystyle \frac{{c}^{2}}{{a}^{2}}\,{\cos }^{2}\psi +\displaystyle \frac{{c}^{2}}{{b}^{2}}\,{\sin }^{2}\psi \right)\left(\displaystyle \frac{t}{\cos \phi }\right)\right],\end{array}\end{eqnarray*}$$

where

$$\begin{eqnarray*}&&{\rm{\Delta }}=1+{c}^{2}\,{\tan }^{2}\phi \left(\displaystyle \frac{{\cos }^{2}\psi }{{a}^{2}}+\displaystyle \frac{{\sin }^{2}\psi }{{b}^{2}}\right).\end{eqnarray*}$$

Substituting these expressions for x, y, and z into the equation for the ellipsoidal surface, we find that the projection of Haumea's limb onto the plane of the sky satisfies

$$\begin{eqnarray*}&&{{Ps}}^{2}+{Qst}+{{Rt}}^{2}=1,\end{eqnarray*}$$

where

$$\begin{eqnarray*}\begin{array}{rcl}P & = & \displaystyle \frac{1}{{{\rm{\Delta }}}^{2}\,{a}^{2}}\,{\sin }^{2}\psi {\left(1+\displaystyle \frac{{c}^{2}}{{b}^{2}}{\tan }^{2}\phi \right)}^{2}\\ & & +\,\displaystyle \frac{1}{{{\rm{\Delta }}}^{2}\,{b}^{2}}\,{\cos }^{2}\psi {\left(1+\displaystyle \frac{{c}^{2}}{{a}^{2}}{\tan }^{2}\phi \right)}^{2}\\ & & +\,\displaystyle \frac{{\tan }^{2}\phi }{{{\rm{\Delta }}}^{2}\,{c}^{2}}\,{\sin }^{2}\psi \,{\cos }^{2}\psi {\left(\displaystyle \frac{{c}^{2}}{{b}^{2}}-\displaystyle \frac{{c}^{2}}{{a}^{2}}\right)}^{2},\end{array}\end{eqnarray*}$$

$$\begin{eqnarray*}\begin{array}{rcl}R & = & \displaystyle \frac{1}{{{\rm{\Delta }}}^{2}\,{\cos }^{2}\phi }\left(\displaystyle \frac{{\cos }^{2}\psi }{{a}^{2}}+\displaystyle \frac{{\sin }^{2}\psi }{{b}^{2}}\right)\\ & & +\,\displaystyle \frac{{c}^{2}\,{\tan }^{2}\phi }{{{\rm{\Delta }}}^{2}\,{\cos }^{2}\phi }{\left(\displaystyle \frac{{\cos }^{2}\psi }{{a}^{2}}+\displaystyle \frac{{\sin }^{2}\psi }{{b}^{2}}\right)}^{2},\end{array}\end{eqnarray*}$$

and

$$\begin{eqnarray*}\begin{array}{rcl}Q & = & \displaystyle \frac{2\,\sin \psi \,\cos \psi }{{{\rm{\Delta }}}^{2}\,{c}^{2}\,\cos \phi }\left(\displaystyle \frac{{c}^{2}}{{b}^{2}}-\displaystyle \frac{{c}^{2}}{{a}^{2}}\right)\\ & & \times \,\left[1+{c}^{2}\,{\tan }^{2}\phi \left(\displaystyle \frac{{\cos }^{2}\psi }{{a}^{2}}+\displaystyle \frac{{\sin }^{2}\psi }{{b}^{2}}\right)\right].\end{array}\end{eqnarray*}$$

These also define an ellipse, rotated in the s–t plane.

After rotating the ellipse in the plane of the sky by an angle θ, defined by

$$\begin{eqnarray*}&&\tan 2\theta =Q/(R-P),\end{eqnarray*}$$

we find that it has axes a' and b' defined by

$$\begin{eqnarray*}&&\displaystyle \frac{1}{{\left(a^{\prime} \right)}^{2}}=\left[P\,{\cos }^{2}\theta +R\,{\sin }^{2}\theta -Q\,\sin \theta \,\cos \theta \right]\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&\displaystyle \frac{1}{{\left(b^{\prime} \right)}^{2}}=\left[P\,{\sin }^{2}\theta +R\,{\cos }^{2}\theta +Q\,\sin \theta \,\cos \theta \right].\end{eqnarray*}$$

We have written a simple code that takes a, b, c, ψ, and ϕ as inputs and solves for θ and then the semiaxes a' and b' of Haumea's shadow.

One end-member case includes ϕ = 0°, in which Haumea's pole is pointed toward the star; we derive θ = −ψ, and regardless of ψ, Haumea's shadow has axes a' = b and b' = a. Another end-member case is ϕ = 90°, in which case the line of sight to the star is parallel to Haumea's equator. The shadow will have b' = c regardless of ψ, and the other axis will be

$$\begin{eqnarray*}&&a^{\prime} ={ab}{\left[\displaystyle \frac{{\cos }^{2}\psi }{{a}^{2}}+\displaystyle \frac{{\sin }^{2}\psi }{{b}^{2}}\right]}^{1/2},\end{eqnarray*}$$

in which case a' = b if ψ = 0° (looking along the long a-axis), or a' = a if ψ = 90° (looking along the b-axis). One more end-member case is ψ = 0° but arbitrary ϕ, in which case a' = b and

$$\begin{eqnarray*}&&b^{\prime} =a\,\cos \phi {\left[1+\displaystyle \frac{{c}^{2}}{{a}^{2}}{\tan }^{2}\phi \right]}^{+1/2}.\end{eqnarray*}$$

This is the case considered by Ortiz et al. (2017). Assuming a = 1161 km, b = 852 km, c = 513 km, ψ = 0°, and ϕ = 763 (a tilt of Haumea's pole with respect to the plane of the sky of 137), we find a' = 852 km and b' = 584 km, similar to the solution found by Ortiz et al. (2017).