© The Authors 2026

1 Introduction

Following the classical result of MacLaurin (1742), who showed that a homogeneous oblate ellipsoid of revolution is a figure of equilibrium, it was long thought that only axisymmetric configurations can describe a self-gravitating rotating fluid mass. Nearly a century later, Jacobi (1834) proved that homogeneous triaxial ellipsoids may also satisfy hydrostatic equilibrium. Jacobi’s theorem establishes that a rigidly rotating homogeneous ellipsoid with semiaxis lengths a > b > c is an equilibrium figure if and only if for given values of a and b, the rotation axis coincides with the c-axis, whose length has a unique value (see Fig. 1). Meyer (1842) further showed that the Maclaurin sequence and the Jacobi sequence share a common member at c/a ≈ 0.58272. Beyond this Meyer bifurcation point, that is, for c/a ≲ 0.58272, the Maclaurin sequence becomes unstable, and any axisymmetry-breaking perturbation distorts the Maclaurin spheroid into a Jacobi ellipsoid¹.

For heterogeneous bodies, fundamental results on the existence of barotropic triaxial figures were obtained by James (1964) and Ipser & Managan (1981). These figures have mostly been studied in the context of galactic dynamics (see, e.g., Vandervoort 1980) or tidally distorted bodies (e.g., Hachisu 1986; Even & Tohline 2009; Wahl et al. 2017).

More recently, light-curve analyses and stellar-occultation shape measurements of the trans-Neptunian objects Haumea and Quaoar have renewed interest in these heterogeneous equilibrium figures. The shape of Haumea was constrained through stellar occultation by Ortiz et al. (2017), who determined that the best-fitting ellipsoid had a = 1161 ± 30 km, b = 852 ± 4 km, and c = 513 ± 16 km. Using another equation set that took photometry into account, Kondratyev & Kornoukhov (2018) obtained the values a = 1082 ± 15 km, b = 836 ± 5 km, and c = 511 ± 13 km. For Quaoar, Kiss et al. (2024) have recently developed a thermophysical model of its surface and found that the best match to the observed thermal light curve corresponds to a triaxial figure with a = 1.29 × 10³ km, b = 1.08 × 10³ km, and c = 9.32 × 10² km. However, as illustrated in Fig. 1, the resulting axis ratios b/a and c/a for both bodies are inconsistent with Jacobi ellipsoids. To explain this discrepancy in the case of Haumea, Ortiz et al. (2017) proposed two possibilities: either Haumea is not in hydrostatic equilibrium, meaning rigid-body forces are required to properly model its shape, or it is not homogeneous, but instead composed of multiple layers. The same alternatives might be considered for Quaoar.

The likelihood that Haumea is hydrostatic was studied by Dunham et al. (2019) and Noviello et al. (2022) in the case of a two-layer planet with a silicate core embedded in an icy mantle. These authors explored how the internal differentiation enlarges the range of possible equilibrium shapes. Using the self-consistent-field (SCF) method, in particular, the threedimensional algorithm of Hachisu (1986), they obtained configurations that matched the light curve derived by Lockwood et al. (2014), but remained inconsistent with, but closer to, the shape inferred by Ortiz et al. (2017). Consequently, it is still unclear whether a differentiated hydrostatic model can fully reproduce the current observations of Haumea, or if a hydrostatic configuration is likely for Quaoar. If models like this were to be found, they would impose strong constraints on the layer compositions and relative sizes.

However, the algorithm of Hachisu (1986) uses a regular grid in spherical coordinates, which requires a high radial resolution to precisely determine the shape of the layer boundaries. This constraint leads to high computational costs and can be an obstacle to a detailed exploration of the parameter space. We aim to introduce a numerical method designed to overcome this limitation and to investigate the equilibrium figures of triaxial differentiated bodies composed of two or more homogeneous layers. In particular, we fully exploit the hypothesis of a piecewise homogeneous system to optimize the number of radial points needed (i.e., one per layer), and we thus greatly reduce the computation time compared to classical SCF methods.

The paper is organized as follows: In Sect. 2 we detail the assumptions and the equation set. In Sect. 3 we adapt the SCF method to the problem of differentiated bodies. The code BALEINES is presented in Sect. 4, and we report successful comparisons with analytical solutions from the Maclaurin and Jacobi sequences and numerical solutions from Hachisu (1986), Descamps (2015), Dunham et al. (2019), and Noviello et al. (2022). In Sect. 5 we discuss equilibrium sequences of two-layer bodies and the effect of differentiation on the axisymmetric-triaxial bifurcation, with an application to the shape of Quaoar. Our concluding remarks and some perspectives are given the last section.

Fig. 1 Possible axis ratio couples (b/a, c/a) for ellipsoidal figures of equilibrium, with the Jacobi sequence (a > b > c) as a black line and the Maclaurin sequence (a = b > c) as a gray line. Haumea and Quaoar are placed in the diagram for illustration. The data and error bars for Haumea come from the observations of Ortiz et al. (2017), and the data for Quaoar come from the thermophysical model of Kiss et al. (2024).

2 Theoretical background

2.1 Classical SCF method

The SCF method originated in quantum chemistry with the Hartree-Fock method, which provides an approximation for the wave function and the fundamental level energy of multi-atom systems. It was then adapted to the context of a self-gravitating rotating fluid by Ostriker & Mark (1968). Major improvements were brought by Hachisu (1986) to better deal with very fast rotators.

The equilibrium of a rigidly rotating and self-gravitating fluid is given by a Bernoulli-like equation, $$\int \frac{\mathrm{d}P}{\rho }-\frac{1}{2}(\mathbf{\Omega }\times \mathbf{r}{)}^2+\mathrm{\Psi }(\mathbf{r})=C,$$(1)

where P is the thermodynamic pressure, ρ is the mass density, Ω is the rotation rate vector, r is the position vector, and C is a constant to be determined, usually from a point on the free surface of the fluid. The gravitational potential, Ψ, is given by $$\mathrm{\Psi }(\mathbf{r})=-G{\iiint }_{\mathcal{V}}{\mathrm{d}}^3{\mathbf{r}}^{\prime }\frac{\rho ({\mathbf{r}}^{\prime })}{|\mathbf{r}-{\mathbf{r}}^{\prime }|},$$(2)

where V is the volume of the fluid, and G is the gravitational constant.

Thus, as remarked by Odrzywołek (2003), Eq. (1) is in fact a nonlinear Hammerstein integral equation with unknown p; this is then a fixed-point problem. When the body greatly deviates from a sphere, this problem must be solved with numerical iterative methods. After the choice of an initial guess, the SCF method can be summarized as a cycle with two main steps. First, we compute Ψ to obtain a new guess on the mass density through Eq. (1), and second, we use this new guess on ρ to compute a new gravitational potential. The procedure is repeated until the solution has converged. The main difference in the implementation of this method is the computation of the gravitational potential, which is a great problem because of the point-mass divergence in the integrand (see Eq. (2)). This problem is usually avoided either by using the multipolar expansion of the potential (Ostriker & Mark 1968; Hachisu 1986; Kiuchi et al. 2010) or by solving the Poisson equation ∇²Ψ = 4πGρ (Kadam et al. 2016; Huré & Hersant 2017).

2.2 The gravitational potential of homogeneous bodies

In homogeneous bodies, the mass density is a constant, and it might seem as if the SCF method became trivial. This is clearly not so simple because the volume of the fluid, particularly the shape of the free surface, is unknown and the classical SCF method still uses the same cycle. The main difficulty in this case is the detection of the free surface, which becomes a key point for the numerical precision. When classical grids are used, the free surface can be described precisely using a coded Freeman chain (Huré & Hersant 2017). Another approach, introduced by Ansorg et al. (2003) in the context of axisymmetric bodies, is to use the assumption of homogeneity to use the Kovalevskaya line integral for the potential (Kowalewsky 1885). As a consequence, the method does not solve for the mass density at each node of a grid, which would give the shape of the free surface, but for the shape itself. This framework then makes the problem one-dimensional and greatly reduces the computation time with respect to classical two-dimensional axisymmetric approaches. When the body has no axial symmetry, a similar transformation can be made. When we assume that ρ is constant inside the fluid, the volume integral in Eq. (2) can indeed be reduced to a surface integral (Gauss 1867; Tisserand 1891), $$\mathrm{\Psi }(\mathbf{r})=\frac{1}{2}G\rho {∯}_{\mathcal{S}}{\mathrm{d}}^2{\mathbf{r}}^{\prime }\mathbf{n}\cdot \frac{\mathbf{r}-{\mathbf{r}}^{\prime }}{|\mathbf{r}-{\mathbf{r}}^{\prime }|},$$(3)

where S is the closed surface bounding V, and n is the outward normal to S at r′. This result can be proved with the divergence theorem. A first consequence of this expression is that the kernel remains finite at coincidence, that is, for r = r′. Indeed, r′ is restricted to the surface, and as r′ → r, r - r′ thus tends to a vector tangent to the surface at r, which is clearly orthogonal to the normal n; the kernel is thus null at coincidence. As a consequence, the point-mass divergence problem is avoided, and direct integration can be made if r belongs to the surface.

2.3 Hydrostatic equilibrium of piecewise homogeneous systems

We considered a differentiated body made of L homogeneous layers. We labeled each layer with an index ℓ, where ℓ = 1 is the innermost layer and ℓ = L is the outermost one. Let S_L be the free surface of the fluid and S_ℓ be the interface between layers ℓ and ℓ + 1 (1 ≤ ℓ ≤ L - 1). As Ψ is linear with the mass density, the potential of a piecewise homogeneous mass distribution can be deduced from the superposition principle, namely $$\mathrm{\Psi }(\mathbf{r})=\frac{1}{2}G\sum _{\ell =1}^L({\rho }_{\ell }-{\rho }_{\ell +1}){∯}_{{\mathcal{S}}_{\ell }}{\mathrm{d}}^2{\mathbf{r}}^{\prime }\mathbf{n}\cdot \frac{\mathbf{r}-{\mathbf{r}}^{\prime }}{|\mathbf{r}-{\mathbf{r}}^{\prime }|},$$(4)

where ρ_ℓ is the mass density of layer ℓ. We set ρ_{L + 1} = 0 for convenience.

We assumed that the polar axis coincides with the Cartesian (Oz) axis, O being the barycenter of the fluid, while the equatorial major and minor axes coincide with (Ox) and (Oy), respectively. We worked with spherical coordinates (r, θ, ϕ), where r is the distance to the center of mass of the body, θ is the colatitude, and ϕ is the azimuthal angle. These coordinates are suitable for almost ellipsoidal figures, that is, for figures in which any ray drawn from the center of mass intersects the boundary exactly once. This is clearly not the case for toroids². With this assumption, we represented the boundary S_ℓ by a function S_ℓ(θ, ϕ), so that the equation of S_ℓ is r = S_ℓ(θ, ϕ). The potential reads $$\begin{array}{rl}& \mathrm{\Psi }(r,\theta ,\phi )\\ & =\frac{1}{2}G\sum _{\ell =1}^L({\rho }_{\ell }-{\rho }_{\ell +1}){\int }_0^{\pi }\mathrm{d}{\theta }^{\prime }{\int }_0^{2\pi }\mathrm{d}{\phi }^{\prime }\frac{\nu ({\theta }^{\prime },{\phi }^{\prime };s_{\ell };r,\theta ,\phi )}{\sigma ({\theta }^{\prime },{\phi }^{\prime };s_{\ell };r,\theta ,\phi )},\end{array}$$(5)

where $$\begin{array}{rl}\nu & ({\theta }^{\prime },{\phi }^{\prime };s_{\ell };r,\theta ,\phi )\\ & =r{s}_{\ell }^2({\theta }^{\prime },{\phi }^{\prime })\sin ({\theta }^{\prime })[\sin (\theta )\sin ({\theta }^{\prime })\cos (\phi -{\phi }^{\prime })+\cos (\theta )\cos ({\theta }^{\prime })]\\ & -r{s}_{\ell }({\theta }^{\prime },{\phi }^{\prime })\frac{\mathrm{\partial }{s}_{\ell }}{\mathrm{\partial }{\theta }^{\prime }}\sin ({\theta }^{\prime })[\sin (\theta )\cos ({\theta }^{\prime })\cos (\phi -{\phi }^{\prime })\\ & -\cos (\theta )\sin ({\theta }^{\prime })]\\ & -r{s}_{\ell }({\theta }^{\prime },{\phi }^{\prime })\frac{\mathrm{\partial }{s}_{\ell }}{\mathrm{\partial }{\phi }^{\prime }}\sin (\theta )\sin (\phi -{\phi }^{\prime })-s_{\ell }^3({\theta }^{\prime },{\phi }^{\prime })\sin ({\theta }^{\prime })\end{array}$$(6)

and $$\begin{array}{rl}{\sigma }^2& ({\theta }^{\prime },{\phi }^{\prime };s_{\ell };r,\theta ,\phi )=r^2+s_{\ell }^2({\theta }^{\prime },{\phi }^{\prime })\\ & -2r{s}_{\ell }({\theta }^{\prime },{\phi }^{\prime })[\sin (\theta )\sin ({\theta }^{\prime })\cos (\phi -{\phi }^{\prime })+\cos (\theta )\cos ({\theta }^{\prime })].\end{array}$$(7)

The hydrostatic equilibrium is governed by L Bernoulli-like equation, one per layer. Inside layer ℓ, Eq. (1) reads $$\frac{P(\mathbf{r})}{{\rho }_{\ell }}-\frac{1}{2}(\mathbf{\Omega }\times \mathbf{r}{)}^2+\mathrm{\Psi }(\mathbf{r})=C_{\ell },$$(8)

where the integral over the pressure has been readily computed, and C_ℓ is a constant.

In absence of an ambient medium, the pressure vanishes along the free surface S_L, namely $$-\frac{{\mathrm{\Omega }}^2{s}_L^2(\theta ,\phi ){\sin }^2(\theta )}{2}+\mathrm{\Psi }(s_L(\theta ,\phi ),\theta ,\phi )=C_L,$$(9)

with Ω = |Ω| being the rotation rate. Moreover, the pressure is continuous at the interface S_ℓ, which leads to $$-\frac{{\mathrm{\Omega }}^2{s}_{\ell }^2(\theta ,\phi ){\sin }^2(\theta )}{2}+\mathrm{\Psi }(s_{\ell }(\theta ,\phi ),\theta ,\phi )=\frac{{\rho }_{\ell }{C}_{\ell }-{\rho }_{\ell +1}{C}_{\ell +1}}{{\rho }_{\ell }-{\rho }_{\ell +1}}.$$(10)

Eqs. (9) and (10) merely state that S_ℓ are level surfaces, that is, surfaces at which the pressure and the sum of the centrifugal and gravitational potential are constant.

As Ψ depends only on the shape of the surfaces, the problem of the hydrostatic equilibrium of a differentiated body then consists of finding the set of functions {S_ℓ(θ,ϕ)} such that Eqs. (9) and (10) are satisfied simultaneously for all ℓ; this formulation of the problem clearly recalls the theory of equilibrium figures (e.g., Zharkov & Trubitsyn 1978; Lanzano 1982). Clearly, Eqs. (9) and (10) are nonlinear coupled integral equations for the set of functions {S_ℓ} and form a fixed-point problem. As for the classical SCF method, iterative methods are required to treat the problem in the general case.

3 The self-consistent-field method for differentiated bodies

As only the shapes of surfaces are needed at a given step, only L points are needed in the radial direction to describe the system for a given numerical resolution in (θ, ϕ). For a small number of layers, this sort of adaptive meshing reduces the computation time with respect to a classical SCF method, which needs a great number of points to precisely determine the interfaces between layers. This advantage in time is clearly reduced as the number of layers increases, so that the method probably has similar performances as the classical SCF method for a continuously heterogeneous system (i.e., L ≫ 1)³.

3.1 The dimensionless problem

For the numerical treatment, it is more suitable to work with dimensionless quantities. We define $$\{\begin{array}{rl}& \hat{r}=r/a_L,\\ & {\hat{s}}_{\ell }=s_{\ell }/a_L,\\ & {\hat{\rho }}_{\ell }={\rho }_{\ell }/{\rho }_L,\\ & \hat{\mathrm{\Psi }}=\mathrm{\Psi }/(G{\rho }_L{a}_L^2),\\ & {\hat{\mathrm{\Omega }}}^2={\mathrm{\Omega }}^2/(G{\rho }_L),\\ & {\hat{C}}_{\ell }=C_{\ell }/(G{\rho }_L{a}_L^2),\end{array}$$(11)

where a_L = s_L(π/2,0) is the equatorial major semiaxis of the whole mass. In dimensionless spherical coordinates, the gravitational potential reads $$\hat{\mathrm{\Psi }}(\hat{r},\theta ,\phi )=\frac{1}{2}\sum _{\ell =1}^L({\hat{\rho }}_{\ell }-{\hat{\rho }}_{\ell +1}){\int }_0^{\pi }\mathrm{d}{\theta }^{\prime }{\int }_0^{2\pi }\mathrm{d}{\phi }^{\prime }{\kappa }_{\ell }({\theta }^{\prime },{\phi }^{\prime };\hat{r},\theta ,\phi ),$$(12)

where the kernel κ_ℓ is $${\kappa }_{\ell }({\theta }^{\prime },{\phi }^{\prime };\hat{r},\theta ,\phi )=\frac{\nu ({\theta }^{\prime },{\phi }^{\prime };s_{\ell };r,\theta ,\phi )}{\sigma ({\theta }^{\prime },{\phi }^{\prime };s_{\ell };r,\theta ,\phi )}\times \frac{1}{a_L^2}.$$(13)

For convenience, let $${\hat{\mathrm{\Psi }}}_{\ell }(\theta ,\phi )=\hat{\mathrm{\Psi }}(s_{\ell }(\theta ,\phi ),\theta ,\phi )$$(14)

be the gravitational potential along S_ℓ. Eqs. (9) and (10) can be rewritten as $${\hat{\mathrm{\Psi }}}_{\ell }(\theta ,\phi )-\frac{1}{2}{\hat{\mathrm{\Omega }}}^2{\hat{s}}_{\ell }^2(\theta ,\phi ){\sin }^2(\theta )-\frac{{\hat{\rho }}_{\ell }{\hat{C}}_{\ell }-{\hat{\rho }}_{\ell +1}{\hat{C}}_{\ell +1}}{{\hat{\rho }}_{\ell }-{\hat{\rho }}_{\ell +1}}=0,$$(15)

$\mathrm{\forall }\ell \in [[1,L]]$, where we set ${\hat{C}}_{L+1}=0$ for convenience.

3.2 Description of the cycle

For this SCF method, we started with a set of initial functions $\{{\hat{s}}_{\ell }^{(0)}(\theta ,\phi )\}$. We chose to start with a set of ellipsoidal surfaces with dimensionless semiaxis length ${\hat{a}}_{\ell }\ge {\hat{b}}_{\ell }\ge {\hat{c}}_{\ell }$, that is, $${\hat{s}}_{\ell }^{(0)}(\theta ,\phi )=\frac{1}{\sqrt{{\displaystyle \frac{{\sin }^2(\theta ){\cos }^2(\phi )}{{\hat{a}}_{\ell }^2}}+{\displaystyle \frac{{\sin }^2(\theta ){\sin }^2(\phi )}{{\hat{b}}_{\ell }^2}}+{\displaystyle \frac{{\cos }^2(\theta )}{{\hat{c}}_{\ell }^2}}}}.$$(16)

These initial ellipsoids were taken to be similar, that is, ${\hat{b}}_{\ell }/{\hat{a}}_{\ell }$ and ${\hat{c}}_{\ell }/{\hat{a}}_{\ell }$ were the same for each starting surface. The mass density ratios from one layer to the layer above, ${\hat{\rho }}_{\ell }/{\hat{\rho }}_{\ell +1},\ell \in [[1,L-1]]$, were set as inputs. The other inputs are detailed in Sect. 3.3.

Then, at iteration t > 0 of the cycle, we computed

• (i)

  the gravitational potential along each interface, ${\hat{\mathrm{\Psi }}}_{\ell }^{(t)}$, using Eq. (12);

• (ii)

  the constants ${\hat{C}}_{\ell }^{(t)}$ and the rotation rate squared, ${\hat{\mathrm{\Omega }}}^{2(t)}$, that is, the method detailed in Sect. 3.3; and

• (iii)

  a new set of surfaces $\{{\hat{s}}_{\ell }^{(t)}(\theta ,\phi )\}$ were obtained by searching for the zero of Eq. (15) for all (θ, ϕ) of the grid.

These steps were repeated until convergence, that is, when all quantities no longer varied from one step to the next within a given threshold, r∊₀. Thus, as a criterion, we used $$ϵ=max|{\hat{s}}_{\ell }^{(t+1)}(\theta ,\phi )-{\hat{s}}_{\ell }^{(t)}(\theta ,\phi )|<{ϵ}_0.$$(17)

We set ∊₀ to 10⁻¹⁴, that is, about 100 times the machine epsilon in double precision. This threshold can clearly be set higher (e.g., ∊₀ = 10⁻⁸ or ∊₀ = 10⁻⁶) to reduce the computation time; the value of 10⁻¹⁴ merely ensures that the solution is numerically stable.

3.3 Inputs and computation of the rotation rate and the constants

As mentioned in the previous paragraph, computing the constants ${\hat{C}}_{\ell }$ and the rotation rate squared Ω² is required at each iteration. We first considered a homogeneous body to facilitate the discussion. We then determined a single constant, ${\hat{C}}_1$, and the square rotation rate. Following Hachisu (1986), classical SCF methods hold two points of the free surface fixed by requiring that the Bernoulli-like equation is automatically satisfied at these points, which yields the values of ${\hat{C}}_1$ and ${\hat{\mathrm{\Omega }}}^2$. When L > 1, ${\hat{C}}_L$ and Ω² are still determined by holding two points of the free surface, S_L, fixed as previously. For a deeper layer ℓ, only ${\hat{C}}_{\ell }$ has to be determined. Thus, one point of the surface S_ℓ needs to be held fixed. As a consequence, in addition to the mass density ratios ${\hat{\rho }}_{\ell }/{\hat{\rho }}_{\ell +1}$, the position of L + 1 points, 2 on the free surface and 1 on the other boundaries, are needed to solve the problem.

Alternatively, we might impose the value of the geodetic parameter, m = Ω²R³/(GM), R being the volumetric radius and M the mass, to determine the value of Ω² at each iteration. Therefore, only one point had to be fixed on the free surface in this formulation. We might also wish to impose the volume fraction of each layer to track the evolution of a system with time. As a consequence, no point would be fixed on the inner boundaries because it would be rescaled from one iteration to the next.

The choice of inputs then depends on the physical problem that is to be solved. The implementation of all these formulations is detailed in the following paragraphs.

3.3.1 Imposing the polar axis lengths

In classical SCF methods, the usual choice is to hold the polar semiaxis length of each layer and the equatorial major semiaxis length of the free surface fixed. The latter is implicitly assumed by our choice of normalization. To do this, the rotation rate and the constants are given by $${\hat{\mathrm{\Omega }}}^2={\hat{\mathrm{\Psi }}}_L(\pi /2,0)-{\hat{\mathrm{\Psi }}}_L(0,0),$$(18)

where we used ${\hat{s}}_L^2(\pi /2,0)=1$, and $$\{\begin{array}{rl}& {\hat{C}}_L={\hat{\mathrm{\Psi }}}_L(0,0)\\ & \frac{{\hat{\rho }}_{\ell }{\hat{C}}_{\ell }-{\hat{\rho }}_{\ell +1}{\hat{C}}_{\ell +1}}{{\hat{\rho }}_{\ell }-{\hat{\rho }}_{\ell +1}}={\hat{\mathrm{\Psi }}}_{\ell }(0,0).\end{array}$$(19)

These expressions ensure that the Bernoulli-like equation is always satisfied at the points we wished to hold fixed. The inputs of this formulation are thus the polar semiaxis length of each layer, ${\hat{c}}_{\ell }$, and the mass density ratios ${\rho }_{\ell }/{\rho }_{\ell +1}$.

3.3.2 Imposing the equatorial minor axis lengths

Another possible choice is to hold the equatorial minor semiaxis of each layer and the equatorial major semiaxis length of the free surface fixed. This can be useful when the equatorial minor axis is better constrained that the polar axis, or when the bifurcation from an axisymmetric figure to a triaxial figure is to be studied, which is easier when the equatorial minor semiaxis of the surface is an input. This is described in Sect. 5.2 for two-layer systems. In this case, the expressions are very similar to the previous case. We have $${\hat{\mathrm{\Omega }}}^2=\frac{{\hat{\mathrm{\Psi }}}_L(\pi /2,0)-{\hat{\mathrm{\Psi }}}_L(\pi /2,\pi /2)}{1-{\hat{s}}_L^2(\pi /2,\pi /2)},$$(20)

where $$\{\begin{array}{rl}& {\hat{C}}_L={\hat{\mathrm{\Psi }}}_L(\pi /2,\pi /2)-\frac{1}{2}{\hat{\mathrm{\Omega }}}^2{\hat{s}}_L^2(\pi /2,\pi /2),\\ & \frac{{\hat{\rho }}_{\ell }{\hat{C}}_{\ell }-{\hat{\rho }}_{\ell +1}{\hat{C}}_{\ell +1}}{{\hat{\rho }}_{\ell }-{\hat{\rho }}_{\ell +1}}={\hat{\mathrm{\Psi }}}_{\ell }(\pi /2,\pi /2)-\frac{1}{2}{\hat{\mathrm{\Omega }}}^2{\hat{s}}_{\ell }^2(\pi /2,\pi /2).\end{array}$$(21)

The inputs of this formulation are the equatorial minor semiaxis length of each layer, ${\hat{b}}_{\ell }$ and the mass density ratios ${\rho }_{\ell }/{\rho }_{\ell +1}$.

3.3.3 Geodetic parameter

With the dimensionless quantities we used, the geodetic parameter, m, reads $$\mathtt{m}=\frac{{\hat{\mathrm{\Omega }}}^2{\hat{R}}^3}{\hat{M}},$$(22)

with $\hat{R}=R/a_L$. This value is often used in theory of figures as an input. To do this, the Bernoulli constants must be written $$\{\begin{array}{rl}& {\hat{C}}_L={\hat{\mathrm{\Psi }}}_L(\pi /2,0)-\frac{1}{2}\mathtt{m}\frac{\hat{M}}{{\hat{R}}^3}{\hat{s}}_L^2(\pi /2,0),\\ & \frac{{\hat{\rho }}_{\ell }{\hat{C}}_{\ell }-{\hat{\rho }}_{\ell +1}{\hat{C}}_{\ell +1}}{{\hat{\rho }}_{\ell }-{\hat{\rho }}_{\ell +1}}={\hat{\mathrm{\Psi }}}_{\ell }(\pi /2,0)-\frac{1}{2}\mathtt{m}\frac{\hat{M}}{{\hat{R}}^3}{\hat{s}}_{\ell }^2(\pi /2,0).\end{array}$$(23)

These expressions also hold the points at the end of the equatorial major semiaxis of each layer fixed. As the length of the polar and equatorial minor semiaxes of the free surface may vary from an iteration to another, holding the length of the equatorial major semiaxes fixed probably is the safest choice.

These equations then require the computation of the dimensionless mass, M̂, and the dimensionless volumetric radius, R̂, at each iteration of the cycle. The inputs of this formulation are the equatorial major semiaxis length of each layer, ${\hat{a}}_{\ell }$, the geodetic parameter, m, and the mass density ratios ${\rho }_{\ell }/{\rho }_{\ell +1}$.

3.3.4 Imposing the volume fraction

The volume fraction of the layers, ν_ℓ, can also be imposed to be constant from one iteration to the next. To do this, we first determined the surfaces at which pressure is continuous at each iteration, t, of the cycle. However, there is no reason that the volume of the layer ℓ determined in this way corresponds to a fraction ν_ℓ of the total volume and has to be rescaled accordingly. Let ${\hat{s}}_{\ell }^{(t+1/2)}$ be the surface equation of layer ℓ, where t + 1/2 is an intermediate step between t and t + 1, corresponding to the rescaling procedure. As the dimensionless equatorial major semiaxis length of the free surface is always set to 1, we directly have ${\hat{s}}_L^{(t+1)}={\hat{s}}_L^{(t+1/2)}$, and the total volume reads $${\hat{V}}^{(t+1)}=\frac{1}{3}{\int }_0^{\pi }\mathrm{d}\theta \sin (\theta ){\int }_0^{2\pi }\mathrm{d}\phi {[{\hat{s}}_L^{(t+1)}(\theta ,\phi )]}^3.$$(24)

Similarly, the volume bounded by the surface S_ℓ at the intermediate step t + 1/2 is given by $${\hat{V}}_{\ell }^{(t+1/2)}=\frac{1}{3}{\int }_0^{\pi }\mathrm{d}\theta \sin (\theta ){\int }_0^{2\pi }\mathrm{d}\phi {[{\hat{s}}_{\ell }^{(t+1/2)}(\theta ,\phi )]}^3.$$(25)

As a consequence, the rescaled surface equation of layer ℓ can be written $${\hat{s}}_{\ell }^{(t+1)}(\theta ,\phi )={\hat{s}}_{\ell }^{(t+1/2)}(\theta ,\phi ){\left[(1-\sum _{j=\ell +1}^L{\nu }_j)\frac{{\hat{V}}^{(t+1)}}{{\hat{V}}_j^{(t+1/2)}}\right]}^{1/3}.$$(26)

This formulation needs to be completed with one of the previous formulations to compute the rotation rate and the constants; in addition to its inputs, imposing the volume fraction requires the value of the volume fraction of each layer, ν_ℓ. This allowed us to study sequences of equilibrium in which the layers have the same volume (or equivalently, the same mass), which could model the quasi-static evolution of a hydrostatic body with a varying rotation rate. Some examples of such sequences are given in Sect. 5.1.

3.4 Numerical computation of the gravitational potential

As mentioned in the previous section, the main advantage of the surface integral in Eq. (3) is that the integrand is fully finite and allows us to compute the gravitational potential on the surface directly. However, this advantage comes with two problems.

First, the integrand is continuous and finite on the surface, but its derivatives are discontinuous at r = r′, with the consequence that the numerical precision of standard high-accuracy integration methods (Gauss-Legendre, Clenshaw-Curtis, etc.) designed for continuous and infinitely derivable functions is limited. The trapezoidal rule is not affected by this problem, but can only achieve a precision of ∼l/N², where N is the number of nodes used for the integration.

Second, it can be proved that when the potential is not computed exactly on the surface, but slightly above or below it, the kernel κ exhibits a lobe at (θ′, ϕ′) = (θ, ϕ) whose width decreases with |r - r|. Thus, for points that are very close to the surface, this lobe is poorly resolved on the numerical grid, and its contribution is slightly overestimated. This produces a discontinuity of the potential at the surface, which clearly is an artifact produced by the poorly accounted for lobe.

A possible method for solving both problems at once is to split the kernel in two, to isolate the problematic part of the integrand, and integrate it analytically (Huré 2005). When κ_ℓ,0 is the splitting kernel, the method simply reads $$\begin{array}{rl}{∯}_{{\mathcal{S}}_{\ell }}{\mathrm{d}}^2{\mathbf{r}}^{\prime }{\kappa }_{\ell }(\mathbf{r},{\mathbf{r}}^{\prime })=& {∯}_{{\mathcal{S}}_{\ell }}{\mathrm{d}}^2{\mathbf{r}}^{\prime }[{\kappa }_{\ell }(\mathbf{r},{\mathbf{r}}^{\prime })-{\kappa }_{\ell ,0}(\mathbf{r},{\mathbf{r}}^{\prime })]\\ & +{∯}_{{\mathcal{S}}_{\ell }}{\mathrm{d}}^2{\mathbf{r}}^{\prime }{\kappa }_{\ell ,0}(\mathbf{r},{\mathbf{r}}^{\prime }),\end{array}$$(27)

where the last term is expected to be analytically integrable at least once. To address the problem at coincidence at least partly, it can be proved that a the splitting kernel has to share its value and its first two derivatives at (θ, ϕ) with the original kernel. However, this makes the mathematical treatment more difficult, and most of the solutions we found for the moment were not general enough to cover a satisfying number of cases. For now, we chose to focus on the lobe, which is more problematic for the SCF-cycle. To keep a simple method that works well enough in all cases, the splitting kernel in this work was the kernel associated with a sphere, concentric and coincident with the fluid at (θ, ϕ), that is, $$\begin{array}{rl}& {\kappa }_{\ell ,0}(\hat{r},\theta ,\phi ;{\theta }^{\prime },{\phi }^{\prime })\\ & =\frac{\hat{r}{\hat{R}}_{\ell }^2\sin ({\theta }^{\prime })[\sin (\theta )\sin ({\theta }^{\prime })\cos (\phi -{\phi }^{\prime })+\cos (\theta )\cos ({\theta }^{\prime })]}{\sqrt{{\hat{r}}^2+{\hat{R}}_{\ell }^2-2\hat{r}{\hat{R}}_{\ell }[\sin (\theta )\sin ({\theta }^{\prime })\cos (\phi -{\phi }^{\prime })+\cos (\theta )\cos ({\theta }^{\prime })]}},\end{array}$$(28)

where ${\hat{R}}_{\ell }={\hat{s}}_{\ell }(\theta ,\phi )$ is the radius of the aforementioned sphere. Thus, the surface integral of κ_ℓ,0 in Eq. (27) corresponds to the gravitational potential of a sphere, namely $${∯}_{{\mathcal{S}}_{\ell }}{\mathrm{d}}^2{\mathbf{r}}^{\prime }{\kappa }_{\ell ,0}(\mathbf{r},{\mathbf{r}}^{\prime })=-\frac{2}{3}\pi \times \{\begin{array}{rlr}& (3{\hat{R}}_{\ell }^2-{\hat{r}}^2),& \hat{r}\le {\hat{R}}_{\ell },\\ & 2\frac{{\hat{R}}_{\ell }^3}{\hat{r}},& \hat{r}>{\hat{R}}_{\ell }.\end{array}$$(29)

This splitting does not enhance the accuracy of the quadratures when the field point is on the surface, but, as shown in Fig. 2, it allowed us to take the lobe better into account when the field point was slightly above or below the surface.

Fig. 2 Kernel function given by Eq. (13) for an ellipsoid with a = 1, b = 0.9, c = 0.6 at (θ, ϕ) ≈ (0.632,0.132), corresponding to a node of the Chebyshev grid, and for r = 1.001 s(θ, ϕ). The kernel, κℓ, is plotted as a function of ϕ′ for θ′ = θ and zoomed-in on the lobe. The Chebyshev grid has N + 1 = 17 nodes. The dashed line corresponds to the Chebyshev approximation of the kernel, and the dashed-and-dotted line shows the kernel of the concentric and coincident sphere, corresponding to κℓ,0 in Eq. (27).

3.5 Computation of the outputs

After convergence was reached, the outputs of the cycle itself were the dimensionless rotation rate, Ω, the shape of the surfaces, ŝ_ℓ(θ,ϕ), and the Bernoulli constants in each layer, ${\hat{C}}_{\ell }$. They can be used to compute several properties of the solution, such as the mass, M, the moment of inertia with respect to the rotation axis, I, the angular momentum, J, or the pressure at all interfaces and at the center, P_ℓ and P_c respectively. The dimensionless mass, $\hat{M}=M/({\rho }_L{a}_L^3)$, moment of inertia, $\hat{I}=I/({\rho }_L{a}_L^5)$, and angular momentum, $J/[G{\rho }_L^3{a}_L^{10}{]}^{1/2}$, are given by $$\{\begin{array}{rl}\hat{M}& =\frac{1}{3}\sum _{\ell =1}^L({\hat{\rho }}_{\ell }-{\hat{\rho }}_{\ell +1}){\int }_0^{\pi }\mathrm{d}\theta \sin (\theta ){\int }_0^{2\pi }\mathrm{d}\phi {\hat{s}}_{\ell }^3(\theta ,\phi ),\\ \hat{I}& =\frac{1}{5}\sum _{\ell =1}^L({\hat{\rho }}_{\ell }-{\hat{\rho }}_{\ell +1}){\int }_0^{\pi }\mathrm{d}\theta {\sin }^3(\theta ){\int }_0^{2\pi }\mathrm{d}\phi {\hat{s}}_{\ell }^5(\theta ,\phi ),\\ \hat{J}& =\hat{I}\hat{\mathrm{\Omega }}.\end{array}$$(30)

The pressure along surface S_ℓ is given by $$\begin{array}{rl}{\hat{P}}_{\ell }=\frac{P_{\ell }}{G{\rho }_L^2{a}_L^2}& ={\hat{\rho }}_{\ell }[{\hat{C}}_{\ell }-{\hat{\mathrm{\Psi }}}_{\ell }(0,0)]\\ & ={\hat{\rho }}_{\ell }{\hat{\rho }}_{\ell +1}\frac{{\hat{C}}_{\ell }-{\hat{C}}_{\ell +1}}{{\hat{\rho }}_{\ell }-{\hat{\rho }}_{\ell +1}},{\hat{\rho }}_{\ell }\ne {\hat{\rho }}_{\ell +1}\end{array}$$(31)

and the pressure at the center by $${\hat{P}}_{\mathrm{c}}=\frac{P_{\mathrm{c}}}{G{\rho }_L^2{a}_L^2}={\hat{\rho }}_1[{\hat{C}}_1-\hat{\mathrm{\Psi }}(0,0,0)].$$(32)

The ratio of kinetic to gravitational energies, |T/W|, often called the stability indicator in the context of SCF methods, can be computed as the ratio of the dimensionless kinetic energy, $\hat{T}=T/(G{\rho }_L^2{a}_L^5)$, $$\hat{T}=\frac{{\hat{\mathrm{\Omega }}}^2}{10}\sum _{\ell =1}^L({\rho }_{\ell }-{\rho }_{\ell +1}){\int }_0^{\pi }\mathrm{d}\theta {\sin }^3(\theta ){\int }_0^{2\pi }\mathrm{d}\phi {\hat{s}}_{\ell }^5(\theta ,\phi ),$$(33)

and the dimensionless gravitational energy, $\hat{W}=W/(G{\rho }_L^2{a}_L^5)$, $$\hat{W}=\sum _{\ell =1}^L({\hat{\rho }}_{\ell }-{\hat{\rho }}_{\ell +1}){\int }_0^{\pi }\mathrm{d}\theta \sin (\theta ){\int }_0^{2\pi }\mathrm{d}\phi \frac{{\hat{s}}_{\ell }^3(\theta ,\phi )}{5}{\hat{\mathrm{\Psi }}}_{\ell }(\theta ,\phi ),$$(34)

The proof of this last expression, which we did not find in the literature, is given in Appendix A.

4 Numerical tests

Using the equations and the cycle described in the previous sections, we built a Fortran90 code, BALEINES (for BALEINES is an algorithm looking for equilibria iteratively for nested ellipsoids using surface integrals), to solve the shape of the layers of a triaxial differentiated fluid mass at hydrostatic equilibrium. In this code, the equations of the surfaces, ${\hat{s}}_{\ell }$, are expanded on Chebyshev polynomials and the surface integrals are evaluated with Clenshaw-Curtis quadratures. The two methods share the same numerical resolution, N, corresponding to N + 1 points for θ and ϕ at the Chebyshev nodes. The planes (xOy), (yOz), and (zOx) were all assumed to be planes of symmetry⁴, so that all calculations were performed in one-eighth of the fluid mass alone, which reduced the computation time. In this section, we benchmark BALEINES with analytical solutions and other numerical codes.

4.1 Reproduction of the Maclaurin, Jacobi, and dumbbell sequences

A first test for BALEINES is the comparison of its solutions with the analytical Maclaurin and Jacobi sequences. For this test, we set L = 1, and the only input was then the meridian axis ratio c/a. The results are reported in Table 1, where the values of b/a and Ω² obtained with BALEINES are compared with the analytical solution. The deviation from the references is consistently on the same order of magnitude for c/a ≥ 0.30. Doubling the resolution decreases the error by almost an order of magnitude, from ~2 × 10⁻⁴ for N = 8 to ~3 × 10⁻⁶ for N = 32. For Clenshaw-Curtis quadratures, this is underwhelming, but is clearly due to the nonderivable character of the kernel (see the discussion in Sect. 3.4).

For c/a ≤ 0.25, the deviation of the analytical solution is greater than the aforementioned estimate. This is also visible in Fig. 3, where the solutions are reported on a angular momentum squared - rotation rate squared diagram. These equilibria, where the deviation from ellipsoids is significant, belong to the dumbbell sequence (Eriguchi et al. 1982) that branches off the Jacobi sequence at (b/a, c/a) ≈ (0.305,0.265). This is confirmed by the comparison with the solutions obtained by Hachisu (1986) and Descamps (2015) in Fig. 3.

It is interesting to note that the Jacobi sequence is unstable beyond this point, whereas the dumbbell is stable. In the same manner, the Maclaurin sequence is unstable beyond Meyer’s bifurcation point, that is, (b/a, c/a) ≈ (1,0.583), and the Jacobi sequence is stable between these two points. Any attempt to obtain an axisymmetric configuration beyond Meyer’s point first converged to a Maclaurin spheroid to a point where the convergence indicator, e, reached ∼10⁻⁷ before it bounced back up and converged to the Jacobi ellipsoid with the same c/a (see Fig. 4). This behavior is understandable by noting that the numerical precision of the quadratures acts as a perturbation of the free surface of the fluid, leading the SCF-cycle to naturally bifurcate from the unstable to the stable state. The same phenomenon occurred when we tried to obtain a Jacobi ellipsoid with c/a ≲ 0.265: BALEINES naturally converged toward the corresponding dumbbell instead. To obtain these solutions, the convergence threshold can be set higher (typically, ∊₀ = 10⁻⁶ in the example in Fig. 4) to stop the cycle before the bifurcation occurs.

Fig. 3 Comparison of numerical solutions from BALEINES (N = 16) with analytical solutions from the Maclaurin-Jacobi sequence and numerical solutions reported by Hachisu (1986) and Descamps (2015).

Fig. 4 Evolution of the convergence indicator, ∊, and equatorial axis ratio, b/a, during the cycle for an initial ellipsoid with c/a = 0.40, b/a = 1.00, and N = 16.

4.2 Comparison with the kyushu code

As a second test, we compared two-layer solutions obtained with the kyushu code (Dunham et al. 2019; Noviello et al. 2022). This code, based on the algorithm reported by Hachisu (1986), was written to study the hydrostatic structure of a planetoid with a silicate core and an icy mantle. The inputs of these models are the length of the two equatorial semiaxes, a₂ and b₂, the mantle mass density, ρ₂, the mass of the object, M and its rotation rate, Ω. The latter two quantities are not used in BALEINES for the convergence, but can be used after convergence to recover all physical quantities. We have $$\{\begin{array}{rl}{\rho }_2& =\frac{{\mathrm{\Omega }}^2}{G{\hat{\mathrm{\Omega }}}^2},\\ a_2& ={\left(\frac{GM}{{\mathrm{\Omega }}^2}\frac{{\hat{\mathrm{\Omega }}}^2}{\hat{M}}\right)}^{1/3},\end{array}$$(35)

and all other quantities are readily deduced by rescaling. We chose to use the configuration preferred by Dunham et al. (2019) for Haumea and a configuration of the tables of Noviello et al. (2022) as references. The comparisons between BALEINES and kyushu are reported in Tables 2 and 3, respectively. The resolution used in both cases was N = 16. Tests were made with N = 32 as well, but the results were the same within the number of digits reported for the reference solutions. The agreement between the two methods is satisfying, as the deviation is about a few 10⁻³ in Table 2 and below 1% in Table 3, which is the numerical precision stated by Dunham et al. (2019) and Noviello et al. (2022) for kyushu. The largest deviations are found on the mean density, with errors between 1 and 10%. This is because in kyushu, the mean density was computed as if the interfaces between layers and the free surfaces were ellipsoids. As shown in Fig. 5, the volume is overestimated because of this approximation, meaning that the mean density, ρ̄, is underestimated. In the case of the configuration reported by Dunham et al. (2019), the deviation is weak (within 3%), so the estimate of ρ is still correct; this is not the case for the configuration of Noviello et al. (2022), however.

Fig. 5 Equilibrium figures of the free surface and the core-mantle boundary obtained with BALEINES for the configurations reported in Tables 2 (top) and 3 (bottom) in the (xOz) (left) and (yOz) (right) planes. The equivalent ellipsoids are also shown as dashed lines.

5 Equilibrium sequences for differentiated bodies

5.1 Evolution of a two-layer system with increasing polar flattening

The objective of this section is to examine how a second layer alters the homogeneous sequence shown in Fig. 3. These studies were conducted in the case of polytropic fluids (e.g., Hachisu 1986), that is, fluids in which the pressure and the mass density are related by P = Kρ^1+1/n, K and n being positive constants; n is the natural parameter of the sequences for polytropes. We chose to impose the core volume fraction (CVF), that is, the ratio of the volume of the core and the total volume of the object (see Sect. 3.3). Holding this fraction fixed allowed us to keep the mass of both components constant along the sequence, which is interesting as it may be interpreted as the evolution of a two-layer system as its polar flattening increases. However, imposing the CVF has the consequence to slow down the convergence of the SCF-cycle because no point on the core-mantle boundary is fixed.

The sequences are shown in Figs. 6 and 7 for core volume fractions of 80 and 50%, respectively, with the rotation rate squared and the stability indicator as functions of the angular momentum squared. The curves obtained for differentiated bodies resemble their polytropic counterpart (see Hachisu 1986), with a shift toward higher rotations and smaller angular momentum as the body approaches homogeneity, and with open sequences, that is, there is no longer a triaxial-dumbbell bifurcation as the Jacobi-like sequence ends prematurely due to mass shedding. However, it is interesting to note that T/|W| is nearly unchanged by the presence of a deeper layer, as only the shift toward lower values of J² is clearly visible.

Fig. 6 Normalized rotation rate (top) and stability indicator (bottom) against normalized angular momentum for several equilibrium sequences of two-layer systems with a CVF of 80%. The homogeneous stable sequence of the equilibrium (Maclaurin-Jacobi-Dumbbell) is plotted in gray for reference.

5.2 Meyer bifurcation for two-layer systems. Application to Quaoar

As mentioned in the introduction and in Sect. 4, the axisymmetric-triaxial bifurcation occurs for c/a ≈ 0.58272. The bifurcation in heterogeneous bodies was studied by James (1964) and Vandervoort (1980) in the case of polytropic fluids. The former showed that no bifurcation is possible only for n ≲ 0.808 because of mass shedding, while the latter showed that Meyer’s point is shifted toward faster rotations as n increases.

The question of the bifurcation of differentiated bodies is particularly interesting for Quaoar. As mentioned in the introduction, Kiss et al. (2024) recently proposed thermophysical models of the surface of this object, and the model that matched the observed thermal light curve best was a triaxial ellipsoid with b/a ≈ 0.840 and c/a ≈ 0.724. However, the obtained axis ratios are inconsistent with a Jacobi ellipsoid, as shown in Fig. 1, and are very far from the homogeneous sequence. For Quaoar to be hydrostatic, a sequence that bifurcates from axisymmetry at c/a > 0.724 would be necessary. Thus, studying the position of the bifurcation for differentiated systems might help us to conclude on the hydrostaticity of the model of Kiss et al. (2024).

For this study, we used BALEINES in a mode in which we held the points along the equatorial minor semiaxis fixed instead of the polar semiaxis (see Sect. 3.3). We then ran BALEINES with b₂/a₂ = 0.999, b₁/b₂ ∊ [0.05,0.95] with step 0.05, and ρ₁/ρ₂ ∊ [2,9] with step 1.

Figure 8 shows that bifurcation is shifted toward slower rotations (i.e., higher c₂/a₂) for small cores (b₁/b₂ ≲ 0.55-0.65) and toward faster rotations (i.e., lower c₂/a₂) for large cores. As b₁/b₂ tends to either 0 or 1, the fluid is closer to homogeneity and the bifurcation then tends to Meyer’s point as expected. However, the value of c₂/a₂ at the bifurcation deviates by less than 10% from Meyer’s point, with a maximum value of 0.604. This is still clearly not close enough to the required (minimum) value of 0.724. We note that ρ₁/ρ₂ = 9 is already unrealistic for a trans-Neptunian object: this value is higher than the density ratio of iron and water. Even considering an absurdly high value of ρ₁/ρ₂ = 35, the maximum value of c₂/a₂ only reaches 0.626. We can then conclude that the triaxial shape derived by Kiss et al. (2024) is inconsistent with an object at hydrostatic equilibrium. The study of the bifurcation point is clearly not the only method to obtain this conclusion (see, e.g., Kondratyev 2023).

Fig. 7 Same as Fig. 6, but for a CVF of 50%.

Fig. 8 Polar to major equatorial semiaxes of the free surface, c2/a2, corresponding to the axisymmetric-triaxial bifurcation as a function of b1/b2 for several values of ρ1/ρ2. The Jacobi ellipsoid with b/a = 0.999, which has c/a ≈ 0.5824, is plotted as well for reference. Cubic splines were used to draw the plain lines.

6 Discussion

We reported a new numerical code, BALEINES, for studying the interior structure of a rotating triaxial mass composed of homogeneous layers in hydrostatic equilibrium. The efficiency of the method was ensured by an adaptive meshing that optimizes the number of radial points needed, that is, only one per layer. The gravitational potential was obtained via a sum of surface integrals. These surface integrals do not present the point-mass divergence of the classical triple integral, which facilitates the numerical integration. The performances of the code agree very well with analytical solutions and the kyushu code (Dunham et al. 2019; Noviello et al. 2022), which was used to study the interior of the dwarf planet Haumea. We then used BALEINES to study the axisymmetric-triaxial bifurcation of two-layer systems and showed that the thermophysical model of the surface preferred by Kiss et al. (2024) is inconsistent with a mass at hydrostatic equilibrium.

A short-term perspective clearly is to use BALEINES to study whether Haumea is a differentiated body at hydrostatic equilibrium. The adaptive meshing would allow us to scan the parameter space in the two-layer case in a more refined way to search for hydrostatic solutions close to observations or previous models. If the two-layer models were encouraging, it would be a thrilling possibility to explore the three-layer parameter space to constrain the internal structure of Haumea based on its triaxial shape alone.

An external potential in the form of a multipolar expansion, for example, could be introduced in Ψ to study the effect of tides on the figure. This might be of interest for modeling the hydrostatic shape of a differentiated satellite in the context of the JUICE mission (Grasset et al. 2013) toward the moons of Jupiter. In any case, the system has to be tidally locked for the problem to be still independent of time in the corotating frame. This condition is mandatory for using the Bernoulli equation. We might also consider a whole fluid instead of a prescribed potential and study the interaction between the two masses.

Acknowledgements

We are grateful to N. Rambaux and F. Chambat for stimulating discussions. We thank the anonymous referee for bibliographic inputs and comments to improve the paper.

References

Appendix A Gravitational energy of a homogeneous body

The gravitational energy of a self-gravitating body reads $$W=\frac{1}{2}{\iiint }_{\mathcal{V}}{\mathrm{d}}^3\mathbf{r}\rho (\mathbf{r})\mathrm{\Psi }(\mathbf{r}),$$(A.1)

which can be also written as $$W=-{\iiint }_{\mathcal{V}}{\mathrm{d}}^3\mathbf{r}\rho (\mathbf{r})\mathbf{r}\cdot \mathrm{\nabla }\mathrm{\Psi }.$$(A.2)

First, we assume that the body is fully homogeneous and we integrate by parts: $$W=-\rho {\iiint }_{\mathcal{V}}{\mathrm{d}}^3\mathbf{r}\mathrm{\nabla }\cdot (\mathbf{r}\mathrm{\Psi })+3\rho {\iiint }_{\mathcal{V}}{\mathrm{d}}^3\mathbf{r}\mathrm{\Psi }(\mathbf{r}).$$(A.3)

The second term of the right-hand side is clearly 6W. So, using the divergence theorem, we finally obtain $$W=\frac{1}{5}\rho {∯}_{\mathcal{S}}{\mathrm{d}}^2\mathbf{r}(\mathbf{n}\cdot \mathbf{r})\mathrm{\Psi }(\mathbf{r}).$$(A.4)

This expression is attractive in our framework, as it only requires the values of the gravitational potential on the surface of the mass.

All Tables

All Figures

Fig. 1 Possible axis ratio couples (b/a, c/a) for ellipsoidal figures of equilibrium, with the Jacobi sequence (a > b > c) as a black line and the Maclaurin sequence (a = b > c) as a gray line. Haumea and Quaoar are placed in the diagram for illustration. The data and error bars for Haumea come from the observations of Ortiz et al. (2017), and the data for Quaoar come from the thermophysical model of Kiss et al. (2024).

In the text

Fig. 2 Kernel function given by Eq. (13) for an ellipsoid with a = 1, b = 0.9, c = 0.6 at (θ, ϕ) ≈ (0.632,0.132), corresponding to a node of the Chebyshev grid, and for r = 1.001 s(θ, ϕ). The kernel, κℓ, is plotted as a function of ϕ′ for θ′ = θ and zoomed-in on the lobe. The Chebyshev grid has N + 1 = 17 nodes. The dashed line corresponds to the Chebyshev approximation of the kernel, and the dashed-and-dotted line shows the kernel of the concentric and coincident sphere, corresponding to κℓ,0 in Eq. (27).

In the text

	Fig. 3 Comparison of numerical solutions from BALEINES (N = 16) with analytical solutions from the Maclaurin-Jacobi sequence and numerical solutions reported by Hachisu (1986) and Descamps (2015).
In the text

	Fig. 4 Evolution of the convergence indicator, ∊, and equatorial axis ratio, b/a, during the cycle for an initial ellipsoid with c/a = 0.40, b/a = 1.00, and N = 16.
In the text

	Fig. 5 Equilibrium figures of the free surface and the core-mantle boundary obtained with BALEINES for the configurations reported in Tables 2 (top) and 3 (bottom) in the (xOz) (left) and (yOz) (right) planes. The equivalent ellipsoids are also shown as dashed lines.
In the text

	Fig. 6 Normalized rotation rate (top) and stability indicator (bottom) against normalized angular momentum for several equilibrium sequences of two-layer systems with a CVF of 80%. The homogeneous stable sequence of the equilibrium (Maclaurin-Jacobi-Dumbbell) is plotted in gray for reference.
In the text

	Fig. 7 Same as Fig. 6, but for a CVF of 50%.
In the text

	Fig. 8 Polar to major equatorial semiaxes of the free surface, c2/a2, corresponding to the axisymmetric-triaxial bifurcation as a function of b1/b2 for several values of ρ1/ρ2. The Jacobi ellipsoid with b/a = 0.999, which has c/a ≈ 0.5824, is plotted as well for reference. Cubic splines were used to draw the plain lines.
In the text