# Nested spheroidal figures of equilibrium – IV. On heterogeneous configurations

C Staelen, J-M Huré

► To cite this version:

C Staelen, J-M Huré. Nested spheroidal figures of equilibrium – IV. On heterogeneous configurations. Monthly Notices of the Royal Astronomical Society, 2024, 527 (1), pp.863-875. 10.1093/mn-ras/stad3215. hal-04294938

HAL Id: hal-04294938

https://hal.science/hal-04294938v1

Submitted on 24 Nov 2023

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Nested spheroidal figures of equilibrium – IV. On heterogeneous configurations

C. Staelen ⌜ ⚜ and J.-M. Huré ⌜ ⚜

Univ. Bordeaux, CNRS, LAB, UMR 5804, F-33600 Pessac, France

Accepted 2023 October 16. Received 2023 September 19; in original form 2023 June 7

## ABSTRACT

The theory of nested figures of equilibrium, expanded in Papers I and II, is investigated in the limit where the number of layers of the rotating body is infinite, enabling to reach full heterogeneity. In the asymptotic process, the discrete set of equations becomes a differential equation for the rotation rate. In the special case of rigid rotation (from centre to surface), we are led to an integro-differential equation (IDE) linking the ellipticity of isopycnic surfaces to the equatorial mass-density profile. In contrast with most studies, these equations are not restricted to small flattenings, but are valid for fast rotators as well. We use numerical solutions obtained from the self-consistent-field method to validate this approach. At small ellipticities (slow rotation), we fully recover Clairaut’s equation. Comparisons with Chandrasekhar’s perturbative approach and with Roberts’ work based on virial equations are successful. We derive a criterion to characterize the transition from slow to fast rotators. The treatment of heterogeneous structures containing mass-density jumps is proposed through a modified IDE.

**Key words:** gravitation – methods: analytical – planets and satellites: interiors – stars: interiors – stars: rotation.

## 1 INTRODUCTION

Unveiling the internal structure of celestial bodies is a longstanding and fundamental challenge in astrophysics. Theories have emerged three centuries ago, with a principal interest in the Earth’s interior. In the limit of slow rotation, Clairaut (1743) showed that the isopycnic surfaces are spheroids, i.e. ellipsoids of revolution. Using the theory of Maclaurin (1742) for homogeneous spheroids, he obtained a second-order, ordinary differential equation linking the flattening of isopycnics to the mass-density profile. This equation has been more recently extended by Lanzano (1962, 1974), the shape of the external surface being expanded over Legendre polynomials $P_{2n}$ up to the $n$th order. Unfortunately, Clairaut’s equation admits essentially no analytical solutions (with some exceptions, see Tisserand 1891; Marchenko 2000). Slow rotators are accessible from the ‘modified’ Lane-Emden equation in the form of series (Chandrasekhar 1933; Kovetz 1968). Besides, Clairaut’s equation is limited to small flattenings (i.e. to low rotation rates), while many systems do not belong to the category of slow rotators. This is the case of the giant planets in the Solar System. For Jupiter and Ceres, the flattening parameter $f \approx 0.07$, and this is even larger for Saturn (Tricarico 2014; Rambaux, Chambat & Castillo-Rogez 2015). Achernar represents an extreme configuration (Carciofi et al. 2008). New developments remain therefore necessary to model the structure of spinning objects, especially for moderate to fast rotation rates (e.g. Lanzano 1962; Ragazzo 2020).

The determination of the gravitational potential of rotating bodies has always demanded a high analytical effort or substantial com-

putational resources, or both (this exceeds the present context). The spheroidal shape is appealing, as its gravitational potential is known with a closed form (see e.g. Chandrasekhar 1969). Kong, Zhang & Schubert (2015) have investigated the validity of the hypothesis of spheroidal isopycnics by comparing of the ‘true’ shape obtained by numerical means to ‘perfect’ spheroids. They showed that discrepancies are small in amplitude. In fact, this remains true at moderate/large rotation, but unsurprisingly fails close to the mass-shedding limit (Hachisu 1986). Using the gravitational potential of a heterogeneous spheroid, Roberts (1963) used the tensor Virial theorem to derive equations valid for fast rotators, but no self-consistent solutions was produced.

In Huré (2022a,b, hereafter, Paper I and II, respectively), we have investigated the conditions of equilibrium of a piece-wise, heterogeneous system made of $\mathcal{L}$ homogeneous layers bounded by spheroidal surfaces. The theory of nested spheroidal figures of equilibrium (hereafter, NSFoE) assumes that these surfaces stay close to confocality (in the sense of oblate spheroidal coordinates; see Section 2) and that each layer can rotate at its own rate (Véronet 1912; Bizyaev, Borisov & Mamaev 2015). A wide range of configurations is then reachable, from quasi-spheres to very flat, disc-like objects. It must be pointed out that such solutions remain approximate, although the Virial parameter relative to the gravitational energy is very small (Staelen 2022). This is a consequence of Poincaré-Hamy theorems: a rigidly rotating body with a spheroidal stratification is not an exact figure of equilibrium.<fnref n="1" /> In this article, we investigate the solutions in the case where the number $\mathcal{L}$ of layers is infinite, which corresponds

<footnote n="1">¹Only confocal surfaces can lead to an exact equilibrium when all layers rotate in a synchronous manner. This equilibrium requires a mass-density</footnote>

<footnote>*E-mail: clement.staelen@u-bordeaux.fr</footnote>
to a fully heterogeneous body, and for a global rigid rotation. It is therefore a natural continuation of Paper II. Another motivation of the article is the case of moderate/fast rotators, characterized by a significant oblateness or flattening (larger than a per cent typically). It is therefore interesting to see to what extent classical theories, which are often limited to slow rotation, remain valid or fail. In this purpose, it is necessary to compare any analytical result with numerical solutions. In the present case, this is achieved by using the DROP code which solves the problem for a polytopic equation-of-state (EoS), various flattenings and rotation profiles (Huré & Hersant 2017; Basillais & Huré 2021).

It is obvious that the present approach is not supposed to surpass sophisticated models for stars and planets, which are dynamically and thermodynamically more complex that what the hypothesis made here allow. Stars are widely prone to mixing, transport and circulation. Planets, closer to rigid rotation, have a more simple layered structure (except at the very surface) and isopycnics surfaces are believed to be very close to spheroids, as suggested by the inversion of gravitational moments (e.g. Hubbard 2013; Nettelmann et al. 2021). After a brief summary on the theory NSFoE, we show in Section 2 how the discrete set of equations can be converted into a differential equation. For rigid, global rotation, this is equivalent to an integro-differential equation (IDE) for the ellipticity of isopycnics. In a first example, we feed this IDE with the numerical solutions obtained from the self-consistent-field (SCF) method (e.g. Hachisu 1986), and show that this approach is not only coherent but quite accurate. In Section 3, we study the behaviour of the equation in the limit of small flattenings, which happens at slow rotation. In particular, we show that the formalism is fully compatible with classical theories, namely the fundamental second-order differential equation established by Clairaut (1743), the solutions obtained by Chandrasekhar (1933) from the ‘modified’ Lane-Emden equation. We also make a comparison with the equation of Roberts (1963). The question of internal jumps is addressed in Section 4, where we derive a modified IDE and test it. In the concluding section, we propose a criterion characterizing the transition from slow to fast rotators, and give a few perspectives.

# 2 THEORY FOR HETEROGENEOUS BODIES

## 2.1 Equation set for the theory of NSFoE

We adopt the same theoretical background and same notations as in Paper I and Paper II, which can be summarized as follows. We consider $\mathcal{L}$ oblate, non-intersecting spheroidal surfaces $E_j$, $j \in [1, \mathcal{L}]$ with semimajor axis $a_j$, semiminor axis $0 < b_j \le a_j$ and eccentricity

$$ \epsilon_j = \sqrt{1 - b_j^2/a_j^2}, \quad (1) $$

as depicted in Fig. 1. These surfaces define $\mathcal{L}$ layers. We note with index $j \ge 2$ the layer bounded by $E_{j-1}$ and $E_j$ (we then have $a_{j-1} < a_j$ and $b_{j-1} < b_j$). Index 1 corresponds to the deepest layer, bounded by surface $E_1$ only. Furthermore, we assume that each layer $j$ is homogeneous, with mass density $\rho_j$, and rotates rigidly around the Z-axis at a rate $\Omega_j$. A key point in the theory of NFSoE is the possibility of asynchronous motion of layers, i.e. $\Omega_{j-1} \ne \Omega_j$. In this article, however, we will consider a subclass of configurations characterized by synchronous rotations.

inversion, which, for stability reasons, is physically not tenable (Poincaré 1888; Hamy 1890; Volterra 1903).

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Typical configuration for a heterogeneous body with finite number of homogeneous layers ($\mathcal{L}$ in total) bounded by spheroidal surfaces $E_j$.</figcaption></figure>

A fundamental parameter that controls the applicability of the theory is the ‘confocal parameter’ $c_{i,j}$, defined for each pair ($E_i, E_j$) by

$$ c_{i,j} = q_{i,j}^2 \epsilon_i^2 - \epsilon_j^2. \quad (2) $$

This parameter is positive if a surface $E_i$, interior to a surface $E_j$ is, in terms of oblate spheroidal coordinates, more oblate than layer $j$. As quoted in the introduction, only systems with confocal spheroidal surfaces (i.e. $c_{i,j} = 0$ for all pairs) correspond to an exact equilibrium (Poincaré 1888; Hamy 1890). Then, the equilibrium of any layered systems in rigid rotation with non-zero confocal parameters is necessarily approximate. As shown in Paper I and Paper II, equilibria with $|c_{i,j}| \lesssim 0.3$ typically are found to be very close to numerical simulations obtained with the DROP-code (Huré & Hersant 2017; Basillais & Huré 2021). Besides, the c-parameter is generally found to be slightly negative, meaning that isopycnic surfaces tend to be more spherical with depth in the system, or equivalently, that the ellipticity of isopycnics increases from the centre to the surface. However, note that if models of stars and planets mainly agree on such a ‘standard’ stratification, there is no argument or observational proof that definitively rules out a reversal, for some objects. This may depend on physical mechanisms at work and on the formation process. Prolate shapes can be induced by circulations or magnetic fields (e.g. Fujisawa & Eriguchi 2014, and references therein).

The starting point of the present work is (27) of Paper II, which links the properties of all layers together. This is a set of coupled, $\mathcal{L} - 1$ algebraic equations, which read,

$$ \begin{aligned} \frac{\alpha_j \tilde{\Omega}_j^2 - \tilde{\Omega}_{j+1}^2}{\alpha_j - 1} &= \sum_{i=1}^{j-1} \tilde{\rho}_{i+1}(\alpha_i - 1) \frac{\tilde{\epsilon}_i}{\epsilon_i^3} \\ &\quad \times \left[ 2 \arcsin \left( \frac{q_{i,j} \epsilon_i}{\sqrt{1+c_{i,j}}} \right) (1+c_{i,j}) \right. \\ &\quad \left. + (1-2q_{i,j}^2 \epsilon_i^2) \arcsin(q_{i,j} \epsilon_i) - q_{i,j} \epsilon_i \sqrt{1-q_{i,j}^2 \epsilon_i^2} - 2q_{i,j} \epsilon_i \tilde{\epsilon}_j \right] \\ &\quad + \sum_{i=j}^{\mathcal{L}} \tilde{\rho}_{i+1}(\alpha_i - 1) \left[ \mathcal{M}(\epsilon_i) + \frac{2}{\epsilon_i^2} (\epsilon_j^2 - \epsilon_i^2) \right] \end{aligned} \quad (3) $$

for $j < \mathcal{L}$, where

$$ \left\{ \begin{array}{l} \tilde{\rho}_j = \rho_j / \rho_{\mathcal{L}}, \\ \alpha_j = \rho_j / \rho_{j+1}, \\ q_{i,j} = a_i / a_j, \\ \tilde{\epsilon}_j = b_j / a_j, \end{array} \right. \quad (4) $$
and

$$ \mathcal{M}(\epsilon) = (3 - 2\epsilon^2) \frac{\bar{\epsilon}}{\epsilon^3} \arcsin(\epsilon) + 3 - \frac{3}{\epsilon^2}, \quad (5) $$

is Maclaurin's function defined by (see Paper I), and

$$ \tilde{\Omega}_j = \frac{\Omega_j}{\sqrt{2\pi G \rho_{\mathcal{L}}}}. \quad (6) $$

is the dimensionless rotation rate normalized to the mass density of the uppermost layer (G is the gravitational constant). For $j = \mathcal{L}$ (the upper layer), we have

$$ \tilde{\Omega}_{\mathcal{L}}^2 = \sum_{i=1}^{\mathcal{L}} \tilde{\rho}_{i+1}(\alpha_i - 1) \frac{\bar{\epsilon}_i}{\epsilon_i^3} \left[ 2 \arcsin \left( \frac{q_{i, \mathcal{L}} \epsilon_i}{\sqrt{1 + c_{i, \mathcal{L}}}} \right) (1 + c_{i, \mathcal{L}}) \right. \\ \left. + (1 - 2q_{i, \mathcal{L}}^2 \epsilon_i^2) \arcsin(q_{i, \mathcal{L}} \epsilon_i) - q_{i, \mathcal{L}} \epsilon_i \sqrt{1 - q_{i, \mathcal{L}}^2 \epsilon_i^2} \right. \\ \left. - 2q_{i, \mathcal{L}} \epsilon_i \bar{\epsilon}_{\mathcal{L}} \right]. \quad (7) $$

## 2.2 From a discrete set of layers to a continuum

We now seek for equilibrium configurations where $\rho$ is continuous and derivable from the centre to the surface. We first consider configurations without any mass-density jumps (mass-density jumps are considered in Section 5). As each layer $j$ has its specific mass-density, and specific set of confocal parameters $c_{i,j}$, the theory of NSFoE is expected to be capable of such a prolongation, provided these confocal parameters are all 'small' enough. When $\mathcal{L}$ drastically increases, the extension of layer $j$ in the equatorial plane is $a_j - a_{j-1} \equiv \Delta a_j \to 0$. In a similar way, at the polar axis, we have $b_j - b_{j-1} \equiv \Delta b_j \to 0$. Furthermore, the difference in the mass-density between two consecutive layers is

$$ \Delta\rho_j = \rho_{j+1} - \rho_j = \rho_{\mathcal{L}}\Delta\tilde{\rho}_j \\ = \rho_{\mathcal{L}}\tilde{\rho}_{j+1}(\alpha_j - 1) \to 0. \quad (8) $$

In these conditions, (3) can be rewritten as

$$ \begin{aligned} & -\frac{\tilde{\rho}_{j+1}\tilde{\Omega}_{j+1}^2 - \tilde{\rho}_j\tilde{\Omega}_j^2}{\Delta\tilde{\rho}_j} \\ &= \sum_{i=1}^{j-1} \Delta\tilde{\rho}_i \frac{\bar{\epsilon}_i}{\epsilon_i^3} \left[ 2 \arcsin \left( \frac{q_{i,j}\epsilon_i}{\sqrt{1+c_{i,j}}} \right) (1+c_{i,j}) - 2q_{i,j}\epsilon_i\bar{\epsilon}_i \right. \\ & \qquad \left. + (1-2q_{i,j}^2\epsilon_i^2) \arcsin(q_{i,j}\epsilon_i) - q_{i,j}\epsilon_i\sqrt{1-q_{i,j}^2\epsilon_i^2} \right] \\ & \qquad + \sum_{i=j}^{\mathcal{L}} \Delta\tilde{\rho}_i \left[ \mathcal{M}(\epsilon_i) + \frac{2}{\epsilon_i^2}(\epsilon_j^2 - \epsilon_i^2) \right]. \end{aligned} \quad (9) $$

In this form, we see that (9) has the convenient form for the continuous case, as in the limit $\Delta\tilde{\rho}_i \to 0$ the sums over $i$ tend to integrals. To express these integrals, the equatorial radius of layer $j$ is rewritten in the form of the dimensionless, continuous variable

$$ \bar{\omega} \equiv \frac{a_j}{R_e}. \quad (10) $$

Equivalently, this is the semimajor axis of the isopycnic surface $E(\bar{\omega}) = E_j$, normalized to the equatorial radius of the body. In a similar manner, we associate

$$ \bar{\omega}' \equiv \frac{a_i}{R_e}, \quad (11) $$

with the equatorial radius of layer $i$. As long as $\epsilon_i \neq 1$, we have $a_1 \to 0$, otherwise a minimal radius is required. Yet, numerical solutions

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Eccentricity profile (left) prescribed for example. Functions $\kappa^{\text{in}}$ and $\kappa^{\text{out}}$ (right), defined by (A1) and (A2) respectively, versus $\bar{\omega}'$ for $\bar{\omega} \in \{0.00, 0.25, 0.50, 0.75, 1.00\}$ (labelled along the curves). Voided circles pin the value of the functions on $\bar{\omega}' = \bar{\omega}$.</figcaption></figure>

obtained from the DROP-code (see Figs 4–7) show that the deeper the layer, the smaller its flattening, namely $\nabla\epsilon > 0$, in agreement with classical theories. Thus, in this work, we will freely take $a_1 \to 0$, so that $(\bar{\omega}', \bar{\omega}) \in [0, 1]^2$.

In the perspective of a continuous mass-density profile, we must consider situations where the mass density vanishes onto the external surface. In general, the adimensioning used for the discrete theory is not appropriate and must be reconsidered. This is easily corrected. In this purpose, we choose the central mass density $\rho_c = \rho_1$ as the new reference, instead of $\rho_{\mathcal{L}}$. The main parameters of layer $j$ are then

$$ \left\{ \begin{array}{l} \hat{\rho}(\bar{\omega}) \equiv \frac{\rho_j}{\rho_c} \equiv \tilde{\rho}_j \times \frac{\rho_{\mathcal{L}}}{\rho_1}, \\ \hat{\Omega}(\bar{\omega}) \equiv \frac{\Omega_j}{\sqrt{2\pi G \rho_c}} \equiv \tilde{\Omega}_j \times \sqrt{\frac{\rho_{\mathcal{L}}}{\rho_1}}, \\ \epsilon(\bar{\omega}) \equiv \epsilon_j, \\ \bar{\epsilon}(\bar{\omega}) = \sqrt{1 - \epsilon^2(\bar{\omega})} \equiv \bar{\epsilon}_j, \end{array} \right. \quad (12) $$

where $\epsilon(\bar{\omega})$ is the eccentricity of the isopycnic surface $E(\bar{\omega})$ and $\bar{\epsilon}(\bar{\omega})$ is its axis ratio. In a similar way, the confocal parameter is now a continuous variable, namely, from (2)

$$ c(\bar{\omega}', \bar{\omega}) = \frac{\bar{\omega}'^2}{\bar{\omega}^2} \epsilon^2(\bar{\omega}') - \epsilon^2(\bar{\omega}). \quad (\equiv c_{i,j}) \quad (13) $$

With these definitions, the left-hand side (LHS) of (9) becomes

$$ \frac{\tilde{\rho}_{j+1}\tilde{\Omega}_{j+1}^2 - \tilde{\rho}_j\tilde{\Omega}_j^2}{\Delta\tilde{\rho}_j} \to \frac{d(\hat{\rho}\tilde{\Omega}^2)}{d\hat{\rho}}, \quad (14) $$

while the right-hand side (RHS) of (9), more complex, can be written in compact form as

$$ \text{RHS of (9)} \to \int_{\hat{\rho}(0)}^{\hat{\rho}(\bar{\omega})} d\hat{\rho}(\bar{\omega}') \kappa^{\text{in}}(\bar{\omega}', \bar{\omega}) \\ + \int_{\hat{\rho}(\bar{\omega})}^{\hat{\rho}(1)} d\hat{\rho}(\bar{\omega}') \kappa^{\text{out}}(\bar{\omega}', \bar{\omega}), \quad (15) $$

where the two functions $\kappa^{\text{in}}$ and $\kappa^{\text{out}}$ are explicitly given in Appendix A; see (A1) and (A2). Formally, these depend on $\bar{\omega}$, $\bar{\omega}'$, and $\epsilon$, i.e. $\kappa \equiv \kappa(\bar{\omega}', \bar{\omega}; \epsilon)$. As $\epsilon$ depends on $\bar{\omega}'$ or $\bar{\omega}$, then there are only two variables on input. Despite appearances, these functions behave very well over the integration range. Among interesting properties, we have $\kappa^{\text{in}}(\bar{\omega}, \bar{\omega}) = \kappa^{\text{out}}(\bar{\omega}, \bar{\omega})$. This is particularly important and attractive for numerical applications. We also see that $\kappa^{\text{in}}(\bar{\omega}, 0) = 0$ for $\bar{\omega} \neq 0$ ($\kappa^{\text{in}}$ is not defined for this value). The typical shape of these functions is visible in Fig. 2(b), where we have plotted $\kappa^{\text{in}}$ and $\kappa^{\text{out}}$ as functions of $\bar{\omega}'$ for fives values of $\bar{\omega}$. For this example, we have prescribed a parabolic profile for the eccentricity (see Fig. 2a), as observed in many numerical experiments (see, e.g. configuration A discussed below).
It follows from (14) and (15) that (9) reads, in the continuous limit

$$ -\frac{d(\hat{\rho}\hat{\Omega}^2)}{d\hat{\sigma}} = \frac{d\hat{\rho}}{d\hat{\sigma}} \left[ \int_{\hat{\rho}(0)}^{\hat{\rho}(\hat{\omega})} d\hat{\rho}(\hat{\omega}')\kappa^{\text{in}}(\hat{\omega}', \hat{\omega}) + \int_{\hat{\rho}(\hat{\omega})}^{\hat{\rho}(1)} d\hat{\rho}(\hat{\omega}')\kappa^{\text{out}}(\hat{\omega}', \hat{\omega}) \right]. \quad (16) $$

This equation is the main equation of the present problem. It links the eccentricity of the isopycnic surfaces and their mass density (in fact, its derivative) to the variations of the rotation rate. It enables to reach configurations where both the mass density and rotation rate vary smoothly with the equatorial radius. We can apply the same transformation to (7), and we obtain at the surface

$$ -\hat{\Omega}^2(1) = \int_{\hat{\rho}(0)}^{\hat{\rho}(1)} d\hat{\rho}(\hat{\omega}')\kappa^{\text{in}}(\hat{\omega}', 1). \quad (17) $$

We see that the RHS of (16) and (17) coincide for $\hat{\omega} = 1$, within a factor $(d\hat{\rho}/d\hat{\sigma})|_{\hat{\omega}=1}$. So, the two LHS must also coincide for $\hat{\omega} = 1$, which imposes the condition

$$ \left. \frac{d\hat{\Omega}^2}{d\hat{\sigma}} \right|_{\hat{\sigma}=1} = 0. \quad (18) $$

Thus, the squared rotation has an extremum at the surface.

## 2.3 The case of global rigid rotation: the general IDE for the eccentricity

We see that (16) is capable of modelling a wide range of situations, from rigid to differential rotation, and independently, from homogeneous to heterogeneous mass-density profiles. In this work, we focus on rigidly rotating bodies, so we have $d\hat{\Omega} = 0$. It means that (18) is naturally satisfied. Therefore, (16) becomes

$$ -\hat{\Omega}^2 = \int_{\hat{\rho}(0)}^{\hat{\rho}(\hat{\omega})} d\hat{\rho}(\hat{\omega}')\kappa^{\text{in}}(\hat{\omega}', \hat{\omega}) + \int_{\hat{\rho}(\hat{\omega})}^{\hat{\rho}(1)} d\hat{\rho}(\hat{\omega}')\kappa^{\text{out}}(\hat{\omega}', \hat{\omega}). \quad (19) $$

This equation yields the rotation rate of the body once the configuration is known through $\hat{\rho}(\hat{\omega})$ and $\epsilon(\hat{\omega})$. Alternatively, it can be used to constrain the solutions if the rotation law is prescribed in advance. We can take the derivative of (19) with respect to $\hat{\omega}$. In this purpose, we use Leibniz's integral rule, namely, for a given derivable function $g$,

$$ \frac{d}{dx} \int_{x_0}^{x} dy g(x, y) = g(x, x) + \int_{x_0}^{x} dy \frac{\partial}{\partial x} g(x, y), \quad (20) $$

where $x_0$ is a constant. In the present case, it leads to

$$ \int_{\hat{\rho}(0)}^{\hat{\rho}(\hat{\omega})} d\hat{\rho}(\hat{\omega}') \frac{\partial}{\partial \hat{\omega}} \kappa^{\text{in}}(\hat{\omega}', \hat{\omega}) + \int_{\hat{\rho}(\hat{\omega})}^{\hat{\rho}(1)} d\hat{\rho}(\hat{\omega}') \frac{\partial}{\partial \hat{\omega}} \kappa^{\text{out}}(\hat{\omega}', \hat{\omega}) = 0, \quad (21) $$

where we have used the property that $\kappa$ is continuous at $\hat{\omega}' = \hat{\omega}$ (see above). In fact, the partial derivatives can be put in the form

$$ \frac{\partial \kappa^{\text{in}}}{\partial \hat{\omega}} = 4\chi(\hat{\omega}', \hat{\omega}) - 2 \frac{d\epsilon^2}{d\hat{\sigma}} \mu(\hat{\omega}', \hat{\omega}), \quad (22) $$

and

$$ \frac{\partial \kappa^{\text{out}}}{\partial \hat{\omega}} = -2 \frac{d\epsilon^2}{d\hat{\sigma}} \nu(\hat{\omega}'), \quad (23) $$

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Function $\chi$ (left panel) defined by (A3), and functions $\mu$ and $\nu$ (right panel) defined by (A4) and (A5), respectively, versus $\hat{\omega}'$ for $\hat{\omega} \in \{0.25, 0.50, 0.75, 1.00\}$ (labelled along the curves). Voided circles pin the value of the functions on $\hat{\omega}' = \hat{\omega}$. The eccentricity profile used is the same as in Fig. 2.</figcaption></figure>

where $\chi$, $\mu$ and $\nu$ are defined in the Appendix A; see (A3), (A4) and (A5), respectively. Like $\kappa$, these functions depend on 3 quantities, $\hat{\omega}$, $\hat{\omega}'$ and $\epsilon$, but implicitly only on the 2 space variables $\hat{\omega}$ and $\hat{\omega}'$. An illustration is given in Fig. 3 for the parabolic eccentricity profile considered previously. An important property is that $\chi(\hat{\omega}, \hat{\omega}) = 0$ and $\mu(\hat{\omega}, \hat{\omega}) = \nu(\hat{\omega})$, which means that the derivatives of the $\kappa$-functions are equal at the connection, i.e.

$$ \left. \frac{\partial \kappa^{\text{in}}}{\partial \hat{\omega}} \right|_{\hat{\omega}'=\hat{\omega}} = \left. \frac{\partial \kappa^{\text{out}}}{\partial \hat{\omega}} \right|_{\hat{\omega}'=\hat{\omega}}. \quad (24) $$

This is visible in Fig. 2(b). These functions have also a relatively small amplitude, which, again, is very practical for any numerical treatment. From (22) and (23), (21) becomes

$$ 2 \int_{\hat{\rho}(0)}^{\hat{\rho}(\hat{\omega})} d\hat{\rho}(\hat{\omega}') \chi(\hat{\omega}', \hat{\omega}) = \frac{d\epsilon^2}{d\hat{\sigma}} \left[ \int_{\hat{\rho}(0)}^{\hat{\rho}(\hat{\omega})} d\hat{\rho}(\hat{\omega}') \mu(\hat{\omega}', \hat{\omega}) + \int_{\hat{\rho}(\hat{\omega})}^{\hat{\rho}(1)} d\hat{\rho}(\hat{\omega}') \nu(\hat{\omega}') \right], \quad (25) $$

which then links directly the eccentricity of the isopycnic surfaces to their mass density.

A consequence of (25) comes from the case $\hat{\omega} = 0$ (i.e. the centre of the body). Indeed, at this point, we have

$$ \left. \frac{d\epsilon^2}{d\hat{\sigma}} \right|_{\hat{\omega}=0} \times \int_{\hat{\rho}(0)}^{\hat{\rho}(1)} d\hat{\rho}(\hat{\omega}') \nu(\hat{\omega}') = 0. \quad (26) $$

As $\nu(\hat{\omega}') < 0$ $\forall \hat{\omega}'$ [this is seen from its definition; see (A5) in Appendix A] and $d\hat{\rho} < 0$ from stability consideration, (26) yields

$$ \left. \frac{d\epsilon^2}{d\hat{\sigma}} \right|_{\hat{\omega}=0} = 0, \quad (27) $$

whatever the mass density profile.

## 2.4 A note on the condition of immersion

As early quoted in this section, an important hypothesis of the theory of NSFoE is the non-intersection of the interfaces between layers. In the continuous limit, this means that the isopycnic surfaces must not cross each other. Two ellipses intersect if and only if the one with the largest major axis has also the smallest minor axis. So, if the polar radius $b$, i.e. the minor axis, is given by $b(\hat{\omega})/R_e = \hat{\omega}\bar{\epsilon}(\hat{\omega})$, we require

$$ \frac{d[\hat{\omega}\bar{\epsilon}(\hat{\omega})]}{d\hat{\omega}} > 0. \quad (28) $$
It is easily shown from the definitions of $\epsilon$ and $\bar{\epsilon}$ that

$$ \frac{d[\bar{\omega}\bar{\epsilon}(\bar{\omega})]}{d\bar{\omega}} = \bar{\epsilon} - \frac{\bar{\omega}}{2\bar{\epsilon}} \frac{d\epsilon^2}{d\bar{\omega}}. \qquad (29) $$

So, (28) can be written as

$$ \frac{d\epsilon^2}{d\bar{\omega}} < \frac{2}{\bar{\omega}} [1 - \epsilon^2(\bar{\omega})], \qquad (30) $$

which then imposes an upper limit for the eccentricity gradient. Note that requiring $\bar{\omega}(d\epsilon^2/d\bar{\omega}) = 2[1 - \epsilon^2(\bar{\omega})]\forall\bar{\omega} \in [0, 1]$ leads to $\epsilon^2(\bar{\omega}) = 1$, i.e. the body would be infinitely flat.

## 2.5 The particular case of homogeneity: Maclaurin formula recovered

A first check of the (25) is performed by considering the Maclaurin spheroid. In this case $d\hat{\rho} = 0$, and the mass density profile is

$$ \hat{\rho}(\bar{\omega}) = \mathcal{H}(1 - \bar{\omega}), \qquad (31) $$

where $\mathcal{H}$ is Heaviside's step function. In the sense of distributions, the derivative of the mass density is

$$ \frac{d\hat{\rho}}{d\bar{\omega}} = -\delta(1 - \bar{\omega}), \qquad (32) $$

where $\delta$ is Dirac distribution. We can now use (25) to obtain $\epsilon$, which is the only unknown of the problem. However, for a body where the mass density is a constant, the notion of isopycnic surfaces appears as a non-sense. Yet, the Poincaré-Wavre theorem implies that for a body where the rotation rate is constant on cylinders, isopycnic and isobaric surfaces must coincide.<fnref n="2" /> Rigid rotation has a rate which is obviously constant on cylinders, so we can consider $\epsilon$ as the eccentricity of isobaric surfaces. It can be shown that (16) and (25) become

$$ \hat{\Omega}^2 = \kappa^{\text{out}}(1, \bar{\omega}), \qquad (33) $$

and

$$ \frac{d\epsilon^2}{d\bar{\omega}} \times \nu(1) = 0, \qquad (34) $$

respectively, where we used the properties of the Dirac distribution. From (A5), we see that $\nu$ never vanishes, so the only solution to (34) is

$$ \frac{d\epsilon^2}{d\bar{\omega}} = 0 \Rightarrow \epsilon^2(\bar{\omega}) = \epsilon^2(1), \qquad (35) $$

which means that isobaric surfaces are similar spheroids. So, by explicating $\kappa^{\text{out}}$, (33) becomes

$$ \hat{\Omega}^2 = (3 - 2\epsilon_s^2) \frac{\bar{\epsilon}_s}{\epsilon_s^3} \arcsin(\epsilon_s) + 3 - \frac{3}{\epsilon_s^2} \equiv \mathcal{M}(\epsilon_s), \qquad (36) $$

where $\epsilon_s \equiv \epsilon(1)$ and $\bar{\epsilon}_s \equiv \bar{\epsilon}(1)$ are values at the surface. As expected, we fully recover the results from Maclaurin's theory.

## 2.6 Checking the IDE from a numerical reference

Unfortunately, without any prior knowledge on the mass-density profile or the eccentricity, (25) cannot solely be used to determine any internal structure. However, we can test the reliability of the

above approach. In this purpose, we find more practical to rewrite (25) in the form

$$ \frac{d\epsilon^2}{d\bar{\omega}} = \frac{2 \int_{\hat{\rho}(0)}^{\hat{\rho}(\bar{\omega})} d\hat{\rho}(\omega')\chi(\omega', \bar{\omega})}{\int_{\hat{\rho}(0)}^{\hat{\rho}(\bar{\omega})} d\hat{\rho}(\omega')\mu(\omega', \bar{\omega}) + \int_{\hat{\rho}(\bar{\omega})}^{\hat{\rho}(1)} d\hat{\rho}(\omega')\nu(\omega')} \qquad (37) $$

The test then consists in computing both sides of this expression by solutions obtained numerically from a Self-Consistent Field (SCF)-method. We use the DROP code (Huré & Hersant 2017; Basillais & Huré 2021) as the numerical reference. This code, which has been extensively used, solves the full structure of rotating, self-gravitating fluids for a wide range of flattenings, EoSs and rotation profiles.<fnref n="3" /> A fundamental ingredient is the closure relationship between pressure $p$ and mass density $\rho$. In the paper throughout, we use a polytropic EoS, namely

$$ p = K \rho^{1+1/n}, \qquad (38) $$

where $K$ and $n$ (the polytropic index) are positive constants. Once the SCF-cycle has converged, the mass-density $\rho(R, Z)$ is known. Then, we have to determine isopycnic surfaces, denoted $S_j$, from centre to surface. Clearly, isopycnics are not exact spheroids (in general, these are slightly depressed in the middle), but any equilibrium surface $S_j$ crosses the polar axis and equatorial axis respectively at points $A_j(0, b_j)$ and $B_j(a_j, 0)$ (see Fig. 1). From these two points, we can calculate a 'pseudo-eccentricity', basically from (1). This pseudo-eccentricity is then of the form $\epsilon(\bar{\omega})$. We can then use this output, together with the mass-density $\rho(\bar{\omega})$ along the equatorial plane to compute $\chi$, $\mu$, $\nu$ and the derivative of the pseudo-eccentricity $d\epsilon^2/d\bar{\omega}$, and then check (37). We will also compare our results to Clairaut's integral equation, i.e. (43), which will be discussed in Section 3, and to Roberts' result, which is given in our notations in Appendix B; see (B1).

There are four main sources of errors in this kind of numerical test. First, DROP releases numerical solutions whose accuracy depends on the resolution. Second, the determination of equilibrium surface $S_j$ (and then, points $A_j$ and $B_j$) is also not perfect. Next, the integrals in the RHS of (37) are also sensitive to the quadrature scheme, as well as the scheme for the derivative of the pseudo-eccentricities (here, we use 2nd-order schemes). Obviously, we do not expect (37) to be exactly satisfied. In turn, if both sides of this equation are very close for a broad variety of configuration, then it proves the reliability of the IDE.

## 2.7 An example

As a first illustration, we consider a rotating polytrope with surface axis-ratio $\bar{\epsilon}_s = 0.75$ and polytropic index $n = 1.5$, hereafter Configuration A. It corresponds to a fast rotator (for comparison, Achernar has an axis ratio in surface around 0.74; see Domiciano de Souza et al. 2014), which is also one of the structures given in the tables of Hachisu (1986). The mass-density, pseudo-eccentricity, and the deviation of the outermost surface to an exact spheroid are displayed in Figs 4(a)–(c) (left panels). The RHS and LHS of (37) are plotted versus $\bar{\omega}$ in Fig. 4(d). We see that the absolute deviation between these two estimates (panel e) is much less than 1 per cent for most radii, and even mess than 0.1 per cent in the outer part of the body. This agreement is already remarkable as the confocal parameters (centre and surface values) are marginally acceptable (i.e. $c(\bar{\omega}, 1)$

<footnote n="3">³While the classical version of the DROP code is typically second-order accurate in the mesh spacing, we have built an alternative version based spectral method. This enables us to reach much high precision (in a shorter computing time).</footnote>

<footnote n="2">²In fact, the theorem states equivalency between four propositions, two of which are used in this discussion; see e.g. Tassoul (1978).</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Results for configuration A ($\bar{\epsilon}_s$ = 0.75, $n$ = 1.5); see Table 1 for global quantities. Left-hand side panels: outputs from the DROP code, i.e. (a) eccentricity of the isopycnic layers (the dashed thin black line is the quadratic profile used for Figs 2 and 3); (b) radial mass density; (c) the deviation of the external surface from a spheroid. Right-hand side panels: Comparison between this work and the output from DROP, i.e. (d) $d\epsilon^2/d\sigma$ as a function of $\sigma$; (e) decimal logarithm of the gap between the analytical methods and the numerical reference for $d\epsilon^2/d\sigma$; (f) decimal logarithm of the gap between DROP and this work for the rotation rate.</figcaption></figure>

$\epsilon \in [-0.4375, 0]$). The figure also shows that the approximation is also valid for $d\epsilon^2/d\sigma$, as the discrepancy with DROP is also of order $\sim 10^{-3}$ in this case. From panel f, we see that the rotation rate $\hat{\Omega}$ deduced from (19) is not strictly a constant, as would be expected. But, we see that it varies weakly and compares greatly with the rotation rate yielded by DROP, with an error below a per cent. We see that Roberts' equation compare greatly with the numerical reference, except at short (where a divergence is seen) and large radii (with an error of $\sim 3$ per cent).

We have calculated the main global properties of the polytrope, namely the mass *M*, the volume *V* and the angular momentum *J*, the gravitational, kinetic and internal energies, *W*, *T* and *U* respectively (see Appendix C) and compared with the tables of Hachisu (1986). The results are reported in Table 1. We see that the values obtained are slightly overestimated with the present approximation. This is due to the boundary of the fluid, which is below the corresponding spheroidal surface, as seen from Fig. 4(c). Thus, the volume of the fluid, and all volume integrals following, are clearly greater than the outputs of the numerical reference. Furthermore, the value of

<table><caption>**Table 1.** Configuration A and corresponding global quantities. Results from the tables of Hachisu (1986) are reported in the first column.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Configuration A</th></tr><tr><th>Hachisu (1986)</th><th>DROP<sup>†</sup></th><th>This work</th></tr></thead><tbody><tr><td>$\bar{\epsilon}_s$</td><td>&lt;0.750</td><td>&lt;0.750</td><td>&lt;0.750</td></tr><tr><td>$n$</td><td>&lt;1.5</td><td>&lt;1.5</td><td>&lt;1.5</td></tr><tr><td>$M/[\rho_c R_e^3]$</td><td>0.430</td><td>0.43027</td><td>0.43280</td></tr><tr><td>$V/R_e^3$</td><td>3.03</td><td>3.02976</td><td>3.14159</td></tr><tr><td>$\hat{\Omega}^2 \times 2\pi$</td><td>0.227</td><td>0.22663</td><td>*0.22760</td></tr><tr><td>$J/[G\rho_c^3 R_e^{10}]^{1/2}$</td><td>0.0356</td><td>0.03556</td><td>0.03609</td></tr><tr><td>$-W/[G\rho_c^2 R_e^5]$</td><td>0.183</td><td>0.18345</td><td>0.18496</td></tr><tr><td>$T/[G\rho_c^2 R_e^5]$</td><td>0.00847</td><td>0.00846</td><td>0.00861</td></tr><tr><td>$U/[G\rho_c^2 R_e^5]$</td><td>0.167</td><td>0.16652</td><td>0.16697</td></tr><tr><td>$[VP/W]$</td><td>&lt;10<sup>-3</sup></td><td>3 &times; 10<sup>-8</sup></td><td>8 &times; 10<sup>-4</sup></td></tr></tbody></table>

Notes. $\leftarrow$ represents input data.

†SCF-method (Basillais & Huré 2021);

*Averaged, see equation (C7).

<table><caption>**Table 2.** Same legend as Table 1, but for configuration B.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Configuration B</th></tr><tr><th>Hachisu (1986)</th><th>DROP<sup>†</sup></th><th>This work</th></tr></thead><tbody><tr><td>$\bar{\epsilon}_s$</td><td>&lt;0.662</td><td>&lt;0.662</td><td>&lt;0.662</td></tr><tr><td>$n$</td><td>&lt;3.0</td><td>&lt;3.0</td><td>&lt;3.0</td></tr><tr><td>$M/[\rho_c R_e^3]$</td><td>0.0255</td><td>0.02545</td><td>0.02546</td></tr><tr><td>$V/R_e^3$</td><td>2.30</td><td>2.24640</td><td>2.77298</td></tr><tr><td>$\hat{\Omega}^2 \times 2\pi$</td><td>0.0256</td><td>0.02563</td><td>*0.02567</td></tr><tr><td>$J/[G\rho_c^3 R_e^{10}]^{1/2}$</td><td>0.00015</td><td>0.00016</td><td>0.00016</td></tr><tr><td>$-W/[G\rho_c^2 R_e^5]$</td><td>0.00140</td><td>0.00139</td><td>0.00139</td></tr><tr><td>$T/[G\rho_c^2 R_e^5]$</td><td>0.00001</td><td>0.00001</td><td>0.00001</td></tr><tr><td>$U/[G\rho_c^2 R_e^5]$</td><td>0.00137</td><td>0.00137</td><td>0.00137</td></tr><tr><td>$[VP/W]$</td><td>&lt;10<sup>-3</sup></td><td>1 &times; 10<sup>-8</sup></td><td>8 &times; 10<sup>-7</sup></td></tr></tbody></table>

the Virial parameter, i.e. $|VP/W| \approx 8 \cdot 10^{-4} \ll 1$, also validates the approximation in this case.

## 2.8 On critical rotations

We can go further in the comparison by looking at an extreme configuration, i.e. a configuration near the so-called 'critical-rotations' (Hachisu 1986), where matter at the surface is barely bounded to the system. Such objects deviate largely from spheroids and we expect the approximation to fail at this point.

We first consider configuration B, with a 'soft' EoS ($n=3$). The configuration and its global properties are reported in Table 2 and the results are plotted in Fig. 5. Surprisingly, the agreement between the spheroidal approximation and the numerical reference is very good; see Fig. 5(d) and Table 2. For $\omega > 0.2$, we see that the discrepancy is $\lesssim 10^{-3}$ in relative. For shorter radii, the gap is wider, due to the numerical precision of the derivatives, as the values themselves are 'small'; so any discrepancy is amplified. The approximation seems to stay valid at the surface, even though the deviation from a spheroid is large (see panel c). This can be explained by the mass density curve, namely panel a. Indeed, we see that, for $\omega > 0.4$, we have $\rho(\omega) \ll \rho_c$, so the contribution of this part to the gravitational potential (and thus, to the rotation rate and (25)) is negligible. So, as long as the isopycnics for $\omega < 0.4$ are close enough to spheroids, the approximation is still valid.

We also have plotted in panel (d) of Fig. 5 the upper limit of the immersion criterion, i.e. (30). Interestingly, the squared eccentricity
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Same legend as for Fig. 4, but for configuration B ($\bar{\epsilon}_s$ = 0.662, $n$ = 3). The dashed black line represents the condition for immersion, i.e. (30).</figcaption></figure>

<table><caption>**Table 3.** Same legend as Table 1, but for Configuration C.</caption><thead><tr><th rowspan="2"></th><th colspan="3">Configuration C</th></tr><tr><th>Hachisu (1986)</th><th>DROP†</th><th>this work</th></tr></thead><tbody><tr><td>ε<sub>s</sub></td><td>←0.442</td><td>←0.442</td><td>←0.442</td></tr><tr><td>n</td><td>←0.5</td><td>←0.5</td><td>←0.5</td></tr><tr><td>M/[ρ<sub>c</sub>R<sub>e</sub><sup>3</sup>]</td><td>0.767</td><td>0.77072</td><td>0.82697</td></tr><tr><td>V/R<sub>e</sub><sup>3</sup></td><td>1.59</td><td>1.56633</td><td>1.85144</td></tr><tr><td>2πΩ<sup>2</sup></td><td>0.939</td><td>0.94158</td><td>*0.93380</td></tr><tr><td>J/[Gρ<sub>c</sub><sup>3</sup>R<sub>e</sub><sup>10</sup>]<sup>1/2</sup></td><td>0.199</td><td>0.19936</td><td>0.22563</td></tr><tr><td>-W/[Gρ<sub>c</sub><sup>2</sup>R<sub>e</sub><sup>5</sup>]</td><td>0.531</td><td>0.53568</td><td>0.59836</td></tr><tr><td>T/[Gρ<sub>c</sub><sup>2</sup>R<sub>e</sub><sup>5</sup>]</td><td>0.0962</td><td>0.09672</td><td>0.10902</td></tr><tr><td>U/[Gρ<sub>c</sub><sup>2</sup>R<sub>e</sub><sup>5</sup>]</td><td>0.339</td><td>0.34230</td><td>0.35218</td></tr><tr><td>|VP/W|</td><td>&lt;10<sup>-3</sup></td><td>1 × 10<sup>-4</sup></td><td>5 × 10<sup>-2</sup></td></tr></tbody></table>

gradient seems to tend to this limit for $\bar{\omega} = 1$, i.e. at the surface. This would imply that at the critical rotation, we have $(db/da)|_{a=R_e} \to 0$, where we used the physical radii, namely the matter at the pole is crushed.

Another example of critical rotation is displayed in configuration C, where the EoS is ‘hard’ ($n = 0.5$). The configuration and its global properties are reported in Table 3 and the results are plotted in Fig. 6. Here, the agreement between the spheroidal approximation reported here and the numerical reference is not good at all, with a relative error of at least 10 per cent on $d\epsilon^2/d\bar{\omega}$ and the global properties. Only

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Same legend as for Fig. 5, but for configuration C ($\bar{\epsilon}_s$ = 0.442, $n$ = 0.5).</figcaption></figure>

the averaged rotation rate is correct, but we see from Fig. 6(f) that the rate itself is not a constant anymore (with an amplitude of, again, ~10 per cent of the mean value). This disagreement is explained by the large deviation of the external surface to a spheroid, which is not cancelled by the mass density profile, i.e. $\hat{\rho}(\bar{\omega}) \ll 1$ only very close to the surface ($\bar{\omega} = 1$). So, the deviation from a spheroid has here a real impact, as the gravitational potential arising from this mass distribution is significatively different from the one produced by a spheroidally stratified object.

However, we observe once again that the immersion criterion joins with the $d\epsilon^2/d\bar{\omega}$-curve computed from DROP at $\bar{\omega} = 1$, reinforcing our conclusion of the previous example.

# 3 THE LIMIT OF SMALL FLATTENINGS

## 3.1 The IDE at first order

The case of slowly rotating structures is of great importance in the context of planetary and stellar interiors (e.g. Chandrasekhar & Roberts 1963; Zharkov & Trubitsyn 1970). Such situations suppose that the deviation to sphericity is small, i.e. $\epsilon^2(\bar{\omega}) \ll 1$. While the Earth or the Sun can probably be considered as slow rotators, this does not seem to be the case of Jupiter and Saturn. The functions defined by (A3)–(A5) can then be expanded at first order in $\epsilon^2$. So, we obtain
$$ \chi(\omega', \omega) = \frac{\omega'^3}{2\omega^6} [\omega'^2\epsilon^2(\omega') - \omega^2\epsilon^2(\omega)] + \mathcal{O}(\epsilon^4), \quad (39) $$

$$ \mu(\omega', \omega) = -\frac{1}{3} \frac{\omega'^3}{\omega^3} \times \left[ 1 - \epsilon^2(\omega') \left( \frac{1}{2} + \frac{3}{5} \frac{\omega'^2}{\omega^2} \right) + \frac{3}{2} \epsilon^2(\omega) \right] + \mathcal{O}(\epsilon^4), \quad (40) $$

and

$$ \nu(\omega') = -\frac{1}{3} - \frac{2}{15}\epsilon^2(\omega') + \mathcal{O}(\epsilon^4), \quad (41) $$

respectively. Thus, at first order in $\epsilon^2$, (25) becomes

$$ -\frac{1}{3} \frac{d\epsilon^2}{d\omega} \left( \int_{\hat{\rho}(0)}^{\hat{\rho}(\omega)} d\hat{\rho}(\omega') \frac{\omega'^3}{\omega^3} + \int_{\hat{\rho}(\omega)}^{\hat{\rho}(1)} d\hat{\rho}(\omega') \right) \\ \approx \frac{1}{\omega^6} \int_{\hat{\rho}(0)}^{\hat{\rho}(\omega)} d\hat{\rho}(\omega') \omega'^3 [\omega'^2 \epsilon^2(\omega') - \omega^2 \epsilon^2(\omega)]. \quad (42) $$

Note that $d\epsilon^2/d\omega$ is already first order in $\epsilon^2$, so the first order terms arising from $\mu$ and $\nu$ can be neglected.

## 3.2 Clairaut's equation recovered

Except in some particular cases, the mass density vanishes continuously at the surface. By integrating (42) by parts, we obtain

$$ \frac{\omega^6}{3} \frac{d\epsilon^2}{d\omega} \hat{\rho}_m \approx \omega^5 \epsilon^2(\omega) \hat{\rho}_m \\ - \int_0^\omega d\omega' \hat{\rho}(\omega') \left[ 5\omega'^4 \epsilon^2(\omega') + \omega'^5 \frac{d\epsilon^2}{d\omega'} \right], \quad (43) $$

where

$$ \hat{\rho}_m(\omega) = \frac{3}{\omega^3} \int_0^\omega d\omega' \hat{\rho}(\omega') \omega'^2, \quad (44) $$

is classically called the mean density (e.g. Tisserand 1891; Ragazzo 2020), evaluated from the centre to the running radius. In this form, (43) is suitable to eliminate the integral by differentiation. So, we derive a second time with respect to the physical radius $a = R_e\omega$ to obtain

$$ \frac{d^2\epsilon^2}{da^2} + \frac{6}{a} \frac{\hat{\rho}}{\hat{\rho}_m} \frac{d\epsilon^2}{da} + \frac{6}{a^2} \left( \frac{\hat{\rho}}{\hat{\rho}_m} - 1 \right) \epsilon^2 \approx 0. \quad (45) $$

This result clearly recalls the fundamental equation derived by Clairaut (1743), namely

$$ \frac{d^2 f}{db^2} + \frac{6}{b} \frac{\hat{\rho}}{\langle \hat{\rho} \rangle} \frac{df}{db} + \frac{6}{b^2} \left( \frac{\hat{\rho}}{\langle \hat{\rho} \rangle} - 1 \right) f = 0. \quad (46) $$

where $f = 1 - \sqrt{1 - \epsilon^2}$ is the flattening of the isopycnic surface, $b$ is its polar radius and

$$ \langle \hat{\rho} \rangle = \frac{3}{b^3} \int_0^b db' \hat{\rho}(b') b'^2. \quad (47) $$

Let us show that (45) and (46) are fully compatible. At first order in $\epsilon^2$, we have $2f \approx \epsilon^2$ and $b \approx a(1 - \epsilon^2/2)$, so

$$ 2 \frac{df}{db} \approx \frac{d\epsilon^2}{da} \left[ 1 + \frac{1}{2} \left( \epsilon^2 + a \frac{d\epsilon^2}{da} \right) \right] \approx \frac{d\epsilon^2}{da}. \quad (48) $$

Now, as the derivatives and the function $f$ itself are already of first order in $\epsilon^2$, only the 'zeroth' order in $\langle \hat{\rho} \rangle$ is needed. At this order, we have $a \approx b$ and thus $\langle \hat{\rho} \rangle \approx \hat{\rho}_m$. Hence, we conclude that (25)

<table><caption>**Table 4.** Same legend as for Table 1 but for configuration D, which is compared with Chandrasekhar (1933) (first column).</caption><thead><tr><th rowspan="2"></th><th colspan="3">Configuration D</th></tr><tr><th>Chandra (1933)</th><th>DROP<sup>†</sup></th><th>this work</th></tr></thead><tbody><tr><td>ε<sub>s</sub></td><td>0.99008</td><td>←0.990</td><td>←0.990</td></tr><tr><td>n</td><td>←1.0</td><td>←1.0</td><td>←1.0</td></tr><tr><td>M/[ρ<sub>c</sub>R<sub>e</sub><sup>3</sup>]</td><td>1.25799</td><td>1.25807</td><td>1.25806</td></tr><tr><td>V/R<sub>e</sub><sup>3</sup></td><td>4.14641</td><td>4.14683</td><td>4.14690</td></tr><tr><td>Ω̂<sup>2</sup> × 2π</td><td>←0.01671</td><td>0.01671</td><td>*0.01671</td></tr><tr><td>J/[Gρ<sub>c</sub><sup>3</sup>R<sub>e</sub><sup>10</sup>]<sup>1/2</sup></td><td></td><td>0.04244</td><td>0.04244</td></tr><tr><td>-W/[Gρ<sub>c</sub><sup>2</sup>R<sub>e</sub><sup>5</sup>]</td><td></td><td>1.19148</td><td>1.19148</td></tr><tr><td>T/[Gρ<sub>c</sub><sup>2</sup>R<sub>e</sub><sup>5</sup>]</td><td></td><td>0.00274</td><td>0.00274</td></tr><tr><td>U/[Gρ<sub>c</sub><sup>2</sup>R<sub>e</sub><sup>5</sup>]</td><td></td><td>1.18600</td><td>1.18599</td></tr><tr><td>[VP/W]</td><td></td><td>7 × 10<sup>-12</sup></td><td>1 × 10<sup>-5</sup></td></tr></tbody></table>

is equivalent to Clairaut's differential equation in the limit of small flattenings, at first order in $\epsilon^2$. Note that some authors (e.g. Ragazzo 2020) use the mean radius $(a^2b)^{1/3}$ instead of $a$ or $b$. We can show by the same reasoning that the equations would still agree at first order.

## 3.3 An example. Comparison with Chandrasekhar's perturbative approach

To illustrate the compatibility between Clairaut's equation and (25), let us consider the numerical solution computed from DROP for a self-gravitating polytrope with $\bar{\epsilon}_s = 0.99$ and $n = 1$, hereafter configuration D; see Table 4 for the details of the configuration and the associated global quantities. We have $\epsilon_s^2 = 0.0199$, which is expected to be 'small enough' for the expansions made in the previous paragraph to be valid. We can therefore check our expansions as well as Clairaut's equation. The results are presented in Fig. 7 (same panels as for configuration A). We notice that the $\epsilon$-profile is close to a quadratic. We see, again, the excellent agreement between the present approach and Clairaut's equation. Also, we see that the global quantities obtained with the IDE are close to the one obtained with DROP, with between four to six digits shared on the values. As quoted in the introduction, this is not a surprise, as Clairaut (1743) showed that for small deviations from the sphere, i.e. small flattenings, the isopycnic surfaces are ellipses in any meridian plane.

Moreover, slowly rotating polytropes have been studied by many authors, in particular by Chandrasekhar (1933). His approach is based on the Lane-Emden equation, supplemented by a small amplitude, rotational field. The equilibrium is solved in the form of series. Configurations with $n = 1$ (like configuration B) are interesting because the results arising from this theory are purely analytical and offer an interesting opportunity for comparisons. As the dimensionless rotation rate $\Omega^2$ is an input in Chandrasekhar's work (while the axis ratio $\bar{\epsilon}$ is an output), the comparison is performed by injecting the rotation rate provided by DROP into Chandrasekhar's equations. The results are reported in Table 2 (column 2). We see that the comparison is satisfactory, the agreement being much better than 1 per cent. Furthermore, the Virial quantities yielded by the spheroidal approximation are in excellent agreement with the numerical reference.

# 4 INTRODUCTION OF MASS-DENSITY JUMPS: THE MODIFIED IDE

Mass-density jumps are usually associated with a sudden change in the EoS or in the mechanism transporting matter or energy. It
is therefore interesting to render the present method as flexible as possible, and to account for such discontinuities. As often, we consider jumps as zero-thickness transitions, while, in real systems, these have always have certain spatial extension. Inspired by Section 2.5, we can easily introduce mass-density jumps in the present formalism by decomposing the mass-density profile as

$$ \hat{\rho}(\omega) = \sum_{k=1}^{\mathcal{K}} [\hat{\rho}_k(\omega) - \hat{\rho}_{k+1}(\omega)] \mathcal{H}(\omega_k - \omega), \quad (49) $$

where $\mathcal{K}$ is the number of heterogeneous domains and $\hat{\rho}_k(\omega) = \rho_k(\omega)/\rho_c$ is the mass density inside domain number $k$ (we still normalize mass-densities to the central value $\rho_c$). As for the discrete case, we have set $\hat{\rho}_{\mathcal{K}+1}(\omega) = 0$ to keep a single sum, which means that the outer space is the very last domain, with index $\mathcal{K}+1$ and null mass density. There are therefore $\mathcal{K}$ jumps, located at $\omega = \omega_k$, $k \in [1, \mathcal{K}]$. Note that (49) allows for configurations with a surface discontinuity, i.e. at $\omega_{\mathcal{K}} = 1$. The derivative of this profile writes

$$ \frac{d\hat{\rho}}{d\omega} = \sum_{k=1}^{\mathcal{K}} \left[ \frac{d\hat{\rho}_k}{d\omega} - \frac{d\hat{\rho}_{k+1}}{d\omega} \right] \mathcal{H}(\omega_k - \omega) - \sum_{k=1}^{\mathcal{K}} [\hat{\rho}_k(\omega) - \hat{\rho}_{k+1}(\omega)] \delta(\omega_k - \omega). \quad (50) $$

We can thus make use of the properties of the Heaviside and Dirac distributions to generalize (25).

## 4.1 Piece-wise rotation and discontinuity in the ellipticity

Let us consider that each domain $k \in [1, \mathcal{K}]$ rotates rigidly at its own rate $\hat{\Omega}_k$. So, for a given domain $k_0$, we have $\omega \in ]\omega_{k_0-1}, \omega_{k_0}]$, and (16) becomes$^4$

$$ \begin{aligned} - \hat{\Omega}_{k_0}^2 = & \sum_{k=1}^{k_0-1} \int_{\omega_{k_0-1}}^{\omega_k} d\omega' \frac{d\hat{\rho}_k}{d\omega'} \kappa^{\text{in}}(\omega', \omega) \\ & + \int_{\omega_{k_0-1}}^{\omega} d\omega' \frac{d\hat{\rho}_k}{d\omega'} \kappa^{\text{in}}(\omega', \omega) \\ & + \int_{\omega}^{\omega_{k_0}} d\omega' \frac{d\hat{\rho}_k}{d\omega'} \kappa^{\text{out}}(\omega', \omega) \\ & + \sum_{k=k_0}^{\mathcal{K}} \int_{\omega_{k_0-1}}^{\omega_k} d\omega' \frac{d\hat{\rho}_k}{d\omega'} \kappa^{\text{out}}(\omega', \omega) \\ & - \sum_{k=1}^{k_0-1} \frac{\alpha_k - 1}{\alpha_k} \hat{\rho}_k(\omega_k) \kappa^{\text{in}}(\omega_k, \omega) \\ & - \sum_{k=k_0}^{\mathcal{K}} \frac{\alpha_k - 1}{\alpha_k} \hat{\rho}_k(\omega_k) \kappa^{\text{out}}(\omega_k, \omega), \end{aligned} \quad (51) $$

where $\alpha_k = \hat{\rho}_k(\omega_k)/\hat{\rho}_{k+1}(\omega_k)$ is the mass density jump at each interface $k$. A major question concerns the behaviour of this equation when applied to two adjacent domains. To answer this point, we write (51) at $\omega_- = \omega_{k_0} - \Delta\omega$ (inside layer $k_0$) and at $\omega_+ = \omega_{k_0} + \Delta\omega$ (inside layer $k_0 + 1$), with $\Delta\omega > 0$. In the limit where $\Delta\omega \ll 1$, the

difference in the rotation rates between $\omega_-$ and $\omega_+$ satisfies

$$ \begin{aligned} \hat{\Omega}_{k_0}^2 - \hat{\Omega}_{k_0+1}^2 = & 4 \Delta\omega \\ & \times \left\{ \sum_{k=1}^{k_0} \int_{\omega_{k_0-1}}^{\omega_k} d\omega' \frac{d\hat{\rho}_k}{d\omega'} \left[ 2\chi(\omega', \omega_{k_0}) - \frac{d\epsilon^2}{d\omega} \bigg|_{\omega_{k_0}} \mu(\omega', \omega_{k_0}) \right] \right. \\ & - \sum_{k=k_0+1}^{\mathcal{K}} \int_{\omega_{k_0-1}}^{\omega_k} d\omega' \frac{d\hat{\rho}_k}{d\omega'} \frac{d\epsilon^2}{d\omega} \bigg|_{\omega_{k_0}} \nu(\omega') \\ & - \sum_{k=1}^{k_0} \frac{\alpha_k - 1}{\alpha_k} \hat{\rho}_k(\omega_k) \left[ 2\chi(\omega_k, \omega_{k_0}) - \frac{d\epsilon^2}{d\omega} \bigg|_{\omega_{k_0}} \mu(\omega_k, \omega_{k_0}) \right] \\ & \left. + \sum_{k=k_0+1}^{\mathcal{K}} \frac{\alpha_k - 1}{\alpha_k} \hat{\rho}_k(\omega_k) \frac{d\epsilon^2}{d\omega} \bigg|_{\omega_{k_0}} \nu(\omega_k) \right\}, \end{aligned} \quad (52) $$

at first order in $\Delta\omega$. If asynchronous motion is possible, then the RHS of this expression must remain finite when $\Delta\omega \to 0$. We see from (52) that this is possible only if the eccentricity undergoes a discontinuity at $\omega_{k_0}$, namely

$$ \left. \frac{d\epsilon^2}{d\omega} \right|_{\omega_{k_0}} = \frac{\epsilon_{k_0}^2(\omega_{k_0}) - \epsilon_{k_0}^2(\omega_{k_0})}{\Delta\omega}, \quad (53) $$

where $\epsilon_{k_0}$ is the eccentricity profile in the domain $k_0$. Note that (52) cannot be used to quantify this jump, as we assumed a continuous eccentricity to arrive at this point. Indeed, if these jumps are considered from the beginning, they would cause discontinuities in the $\kappa$-functions, $\chi$, $\mu$, and $\nu$, which makes the calculations far more complex.

This 'eccentricity jump' only states that the interfaces between layers are not isopycnic surfaces. The isopycnic in the inner layer (the 'core') intersect the interface and is prolonged by another isopycnic in the outer layer (the 'envelope') whose eccentricity has no reason to be the same. This statement has two interesting consequences: i) the 'eccentricity jump' occurs not at a single value of $\omega$ but on a whole range close to any interface; ii) the potential of an incomplete Maclaurin spheroid being unknown analytically, the continuous version of the NSFoE cannot describe systems with rotational discontinuities.

## 4.2 Global, rigid rotation

By requiring $\hat{\Omega}_k = \hat{\Omega}$, $\forall k \in [1, \mathcal{K}]$, both sides in (52) are null in the limit $\Delta\omega \to 0$, meaning no eccentricity jump occurs for systems in global rotation, so that the interfaces between layers are isopycnic surfaces. As the RHS of (51) is constant, we can, as in the single-layer case, take its derivative with respect to $\omega$ inside layer $k_0$. We find

$$ \begin{aligned} & \sum_{k=1}^{k_0-1} \int_{\omega_{k_0-1}}^{\omega_k} d\omega' \frac{d\hat{\rho}_k}{d\omega'} \chi(\omega', \omega) + \int_{\omega_{k_0-1}}^{\omega} d\omega' \frac{d\hat{\rho}_k}{d\omega'} \chi(\omega', \omega) \\ & - \sum_{k=1}^{k_0-1} \frac{\alpha_k - 1}{\alpha_k} \hat{\rho}_k(\omega_k) \chi(\omega_k, \omega) = \frac{1}{2} \frac{d\epsilon^2}{d\omega} \\ & \times \left\{ \sum_{k=1}^{k_0-1} \int_{\omega_{k_0-1}}^{\omega_k} d\omega' \frac{d\hat{\rho}_k}{d\omega'} \mu(\omega', \omega) + \int_{\omega_{k_0-1}}^{\omega} d\omega' \frac{d\hat{\rho}_k}{d\omega'} \mu(\omega', \omega) \right. \\ & - \sum_{k=1}^{k_0-1} \frac{\alpha_k - 1}{\alpha_k} \hat{\rho}_k(\omega_k) \mu(\omega_k, \omega) - \sum_{k=k_0}^{\mathcal{K}} \frac{\alpha_k - 1}{\alpha_k} \hat{\rho}_k(\omega_k) \nu(\omega_k) \\ & \left. + \sum_{k=k_0+1}^{\mathcal{K}} \int_{\omega_{k_0-1}}^{\omega_k} d\omega' \frac{d\hat{\rho}_k}{d\omega'} \nu(\omega') + \int_{\omega}^{\omega_{k_0}} d\omega' \frac{d\hat{\rho}_k}{d\omega'} \nu(\omega') \right\}. \end{aligned} \quad (54) $$

<footnote>$^4$We have introduced $\omega_0 = 0$ for convenience.</footnote>
<table><caption>**Table 5.** Same legend as for Table 2, but for configuration A'.</caption><thead><tr><th></th><th colspan="2">Configuration A'</th></tr><tr><th></th><th>DROP†</th><th>This work</th></tr></thead><tbody><tr><td>ε<sub>s</sub></td><td colspan="2">← 0.750</td></tr><tr><td>n<sub>1</sub></td><td colspan="2">← 1.5</td></tr><tr><td>n<sub>2</sub></td><td colspan="2">← 3.0</td></tr><tr><td>ω<sub>1</sub>ε(ω<sub>1</sub>)</td><td colspan="2">← 0.35</td></tr><tr><td>ω<sub>1</sub></td><td colspan="2">0.35965</td></tr><tr><td>α<sub>1</sub></td><td colspan="2">← 2.0</td></tr><tr><td>M/[ρ<sub>c</sub>R<sub>e</sub><sup>3</sup>]</td><td>0.05722</td><td>0.05735</td></tr><tr><td>V/R<sub>e</sub><sup>3</sup></td><td>2.92254</td><td>3.14519</td></tr><tr><td>2πΩ<sup>2</sup></td><td>0.03765</td><td>*0.03768</td></tr><tr><td>J/[Gρ<sub>c</sub><sup>3</sup>R<sub>e</sub><sup>10</sup>]<sup>1/2</sup></td><td>0.00048</td><td>0.00048</td></tr><tr><td>-W/[Gρ<sub>c</sub><sup>5</sup>R<sub>e</sub><sup>5</sup>]</td><td>0.00631</td><td>0.00633</td></tr><tr><td>T/[Gρ<sub>c</sub><sup>2</sup>R<sub>e</sub><sup>5</sup>]</td><td>0.00004</td><td>0.00004</td></tr><tr><td>U/[Gρ<sub>c</sub><sup>2</sup>R<sub>e</sub><sup>5</sup>]</td><td>0.00622</td><td>0.00623</td></tr><tr><td>|VP/W|</td><td>8 × 10<sup>-5</sup></td><td>1 × 10<sup>-3</sup></td></tr></tbody></table>

This expression is the IDE modified by the presence of jumps. Note that it can be recast in the form of (37). As for the single-layer case, (54) cannot be solved alone as we have a single IDE for $\mathcal{K} + 1$ unknown functions, namely the mass-density profiles $\hat{\rho}_k(\omega)$ and the eccentricity $\epsilon(\omega)$. A solution requires $\mathcal{K}$ equations of state and $\mathcal{K}$ Bernoulli's equations.

## 4.3 An example

Once again, we check the self-consistency of (54) by comparison with a numerical solution from DROP; see Section 2.6. We see that (54) can be written in the form of (37), i.e. we can obtain an equation of the form $d\epsilon^2/d\sigma = g(\rho, \epsilon, \sigma)$. So, as before, we use DROP outputs to compute both sides of (54) and we then compare the results.

Configuration A' is a rotating body with surface axis ratio of 0.75, a core with polytropic index $n_1 = 1.5$ and semi-polar axis $\omega_1\epsilon(\omega_1) = 0.35$ and an envelope with polytropic index $n_2 = 3$. This system could correspond to a highly flatten object with a large convective core (whose mean radius is $\sim 40$ per cent of the star's radius) and a big radiative envelope; it may thus be considered as a very simple model for a fast-rotating high-mass star ($M \gtrsim 1.2 M_\odot$); see e.g. Maeder (2009). The global quantities are given in Table 5 and the results are plotted in Fig. 8. Again, the agreement between the spheroidal approximation reported here and the numerical reference is remarkable, within a few tenths of a per cent (except for the volume). The relative Virial parameter is also really good ($10^{-3} \ll 1$), which validates more the approach. We see that both squared eccentricity gradients compare really well to each other (see panel d), the discrepancy being around $\sim 10^{-3}$ in most of the object and around $\sim 10^{-2}$ in the neighbouring of the mass density jump, which is due to the numerical resolution in this region. Indeed, for each cylindrical radius, the interface is described by two or three points, which may not be enough to reach a good accuracy on the dynamics of the eccentricity in this region. This peak is also seen in the $\hat{\Omega}^2$ curve (panel f), where the gap to the value yielded by DROP is also about a few tenth of a per cent.

# 5 DISCUSSION

## 5.1 Summary

This article investigates the condition of equilibrium of a heterogeneous system with spheroidal isopycnic surfaces (axisymmetrical

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Same legend as Fig. 4, but for configuration D ($\bar{\epsilon}_s = 0.99$, $n = 1.0$).</figcaption></figure>

case). We have derived the main IDE of the problem in the case where the rotation rate is constant onto the isopycnic surfaces, and we have deduced the corresponding IDE in the special case of rigid rotation. This IDE works for a wide range of rotation rates, not only in the slow rotating limit as often considered. Using the DROP-code as a numerical reference, we have proven the reliability of the approach for various configurations, including fast rotators; see configurations A and B. The IDE is fully compatible with Clauraut's equation in the case of slow rotation. Furthermore, we have seen a correlation between the state of critical rotation and the criterion of non-intersection of the isopycnics. As shown, mass-density jumps can be taken into account in the model as long as there are no rotational discontinuities.

## 5.2 Open questions and perspectives

(i) **Rotational discontinuities.** When rotational discontinuities are present, an eccentricity jump is mandatory, meaning the interfaces between layers were not isopycnic surfaces. The approximation of spheroidal isopycnic then fails in this case. However, if the rotational discontinuities (or equivalently, the eccentricity jumps) are small enough, it should be possible to derive an IDE for this case, as the range where the jump occurs becomes negligible. This point would merit an additional work.

(ii) **From slow to fast rotator: a criterion.** In the limit of small flattenings, our approach compares really well with the one developed by Chandrasekhar (1933) and we were able to recover Clauraut's
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Same legend as for Fig. 4, but for configuration A', which is a two-domain body (i.e. $K = 2$) with a mass-density jump at $\bar{\omega} = \bar{\omega}_1 \approx 0.36$ (marked with a vertical red dashed line; see Table 5).</figcaption></figure>

equation at first order in $\epsilon^2$. This addresses the question of the limit between slow rotators (well described by Clairaut's theory) and fast rotators, which can be roughly answered as follows. Let us develop $\nu$ at second order in $\epsilon^2$ (for convenience, we use $\nu$ instead of $\chi$ or $\mu$ as it is a function of a single variable). From (A5), we directly obtain

$$ \nu(\bar{\omega}) = -\frac{1}{3} - \frac{2}{15}\epsilon^2(\bar{\omega}) - \frac{8}{105}\epsilon^4(\bar{\omega}) + \mathcal{O}(\epsilon^6). \quad (55) $$

Now, let $\eta$ be the ratio of the fourth-order term to the second-order term. We have

$$ \eta = \frac{4}{7}\epsilon^2(\bar{\omega}). \quad (56) $$

So, roughly, the error in the quantity $d\epsilon^2/d\bar{\omega}$ made by using Clairaut's equation, i.e. (45), is of the order of $\eta$. The corresponding axis ratio at the surface is then

$$ \bar{\epsilon}_s \ge \sqrt{1 - \frac{7}{4}\eta}. \quad (57) $$

We see that for configuration D (see Fig. 7d), which has an axis ratio of 0.99 at the surface, the maximum error is of order $10^{-2}$ (we do not take into account the part $\bar{\omega} < 0.2$, which is dominated by the errors of the finite-difference scheme). This then corresponds to the criterion (57).

(iii) Can we expand the IDE at higher orders? As shown, expanding the IDE at first order in $\epsilon^2$ leads to Clairaut's equation. It would then be interesting to derive a second-order Clairaut equation, basically by expanding the IDE to a second-order expansion in $\epsilon^2$. This would be another approach to the expansion of Clairaut's equation than Lanzano (1962, 1974) who has performed a multipolar expansion of the shape of the object. However, preliminary calculations indicate that the problem might not be any easier than the equation set reported here. This point is still under investigation.

(iv) Do exact solutions to the IDE exist? As it is well known, analytical solutions are always powerful tools for making models and diagnosis tools, regarding observations. The existence of analytical solutions to the IDE in the form $\rho(\epsilon)$ would be very interesting, and it already represents an exciting perspective. Clairaut's equation is known to have a few analytical solution (e.g. Tisserand 1891; Marchenko 2000). Given the complexity of the IDE, we expect any analytical solution to be only approximate. Solutions via a series expansion or linearization for example would be interesting to seek for.

(v) Towards 2D structures? As quoted, (25) is not sufficient in itself to derive models for interiors of rotating bodies; it is the case of Clairaut's equation as well. The IDE has to be combined with an EoS and to Bernoulli's equation. However, the IDE enables to reduce the number of dimensions of the problem, from two to one, through the relationship $\rho(\epsilon)$. The computation of the gravitational potential is skipped in this process (in fact, it is already incorporated in the IDE). This is very attractive, in particular in terms of computing time if a large number of structures have to be computed (see below). We are currently preparing an article dealing with the structure of spheroidal stars and planets from an SCF-method (Hachisu 1986) through this dimension reduction.

(vi) Inverse problems. Planets like Jupiter and Saturn do probably not belong to the category of slow rotators. The IDE could therefore be of great help in generating fast internal 2D structures (with appropriate EoS), under the conditions of the hypothesis of the NSFoE. Next, it would be easy to compute the gravitational moments and to isolate solutions that match the values 'measured' by space probes. Yet, as pictured by for example Miguel & Vazan (2023), high-order gravitational moments mostly describe the outer layers of the object, which is the most poorly described zone by the theory reported here; see also Basillais & Huré (2023) (Paper III) and references therein. As such, we expect only the first two moments to be accurate enough. Furthermore, as quoted by Nettelmann et al. (2021), the concentric Maclaurin spheroid (CMS) method by (Hubbard 2013) has high computational needs, meaning that a scan of a given parameter space is tedious. With a very fast algorithm, it could be possible to identify places in the parameter space compatible with the measured $J_{2n}$, which could be further studied with more sophisticated algorithms (e.g. the CMS method). Obviously, in the case of gaseous planets, the presence of complex winds at the very surface is not strictly compatible with the NSFoE (the 3D structure of a gaseous planet with zonal winds has been studied by Kong, Zhang & Schubert 2016). This is worst in stars where meridional circulations are present (see e.g. Zahn 1992).

ACKNOWLEDGEMENTS

We are grateful to A. Albouy, G. Boué, and M. Serrero for stimulating discussions during our visit at IMCCE in May 2023. We thank the anonymous referee for fruitful comments relative to the inverse problem met in planetary science.
DATA AVAILABILITY

All data are incorporated into the article.

REFERENCES

Basillais B., Huré J.-M., 2021, MNRAS, 506, 3773

Basillais B., Huré J.-M., 2023, MNRAS, 520, 1504

Bizyaev I. A., Borisov A. V., Mamaev I. S., 2015, Celest. Mech. Dyn. Astron., 122, 1

Carciofi A. C., Domiciano de Souza A., Magalhães A. M., Bjorkman J. E., Vakili F., 2008, ApJ, 676, L41

Chandrasekhar S., 1933, MNRAS, 93, 390

Chandrasekhar S., 1969, Ellipsoidal Figures of Equilibrium. Yale Univ. Press, New Haven

Chandrasekhar S., Roberts P. H., 1963, ApJ, 138, 801

Clairaut A. C., 1743, Théorie de la figure de la Terre tirée des principes de l'hydrostatique. David Fils, Paris

Domiciano de Souza A. et al., 2014, A&A, 569, A10

Fujisawa K., Eriguchi Y., 2014, MNRAS, 438, L61

Hachisu I., 1986, ApJS, 61, 479

Hamy M., 1890, J. math. pures appl. 4e sér., 6, 69

Hubbard W. B., 2013, ApJ, 768, 43

Huré J.-M., 2022a, MNRAS, 512, 4031 (Paper I)

Huré J.-M., 2022b, MNRAS, 512, 4047 (Paper II)

Huré J.-M., Hersant F., 2017, MNRAS, 464, 4761

Kong D., Zhang K., Schubert G., 2015, Phys. Earth Planet. Inter., 249, 43

Kong D., Zhang K., Schubert G., 2016, ApJ, 826, 127

Kovetz A., 1968, ApJ, 154, 999

Lanzano P., 1962, Icarus, 1, 121

Lanzano P., 1974, Ap&SS, 29, 161

Maclaurin C., 1742, A Treatise of Fluxions. In Two Books. 1, T.W and T. Ruddimans, Edinburgh

Maeder A., 2009, Physics, Formation and Evolution of Rotating Stars. Springer Berlin, Heidelberg

Marchenko A. N., 2000, Astron. School's Rep., 1, 34

Miguel Y., Vazan A., 2023, Remote Sens., 15, 681

Nettelmann N. et al., 2021, Planet. Sci. J., 2, 241

Poincaré H., 1888, C. R. Seances Acad. Sci., 106, 1571

Ragazzo C., 2020, São Paulo J. Math. Sci., 14, 1

Rambaux N., Chambat F., Castillo-Rogez J. C., 2015, A&A, 584, A127

Roberts P. H., 1963, ApJ, 138, 809

Staelen C., 2022, MSc dissertation, Univ. Bordeaux

Tassoul J.-L., 1978, Theory of rotating stars. Princeton Univ. Press, Princeton

Tisserand F., 1891, Traité de mécanique céleste - II. Théorie de la figure des corps célestes et de leur mouvement de rotation. Gauthier-Villars et fils, Paris

Tricarico P., 2014, ApJ, 782, 99

Véronet A., 1912, J. Math. pures appl. 6e sér., 8, 331

Volterra V., 1903, Acta Math., 27, 105

Zahn J. P., 1992, A&A, 265, 115

Zharkov V. N., Trubitsyn V. P., 1970, SvA, 13, 981

APPENDIX A: KERNEL FUNCTIONS

Let us write explicitly the kernel functions of the integrals of the main equations of this work.

$$ \forall \sigma \in [0, 1], \forall \sigma' \in [0, \sigma[, \text{ we have} $$

$$ \kappa^{\text{in}}(\sigma', \sigma) = \frac{\bar{\epsilon}(\sigma')}{\epsilon^3(\sigma')} \left\{ \begin{aligned} & \left[ 1 - 2 \frac{\sigma'^2 \epsilon^2(\sigma')}{\sigma^2} \right] \arcsin \left( \frac{\sigma' \epsilon(\sigma')}{\sigma} \right) \\ & - \frac{\sigma' \epsilon(\sigma')}{\sigma} \left( 2 \bar{\epsilon}(\sigma) + \sqrt{1 - \frac{\sigma'^2 \epsilon^2(\sigma')}{\sigma^2}} \right) \\ & + 2 \arcsin \left( \frac{\sigma' \epsilon(\sigma')}{\sigma \sqrt{1 + c(\sigma', \sigma)}} \right) \end{aligned} \right\}. \quad (A1) $$

$$ \forall \sigma \in [0, 1], \forall \sigma' \in [\sigma, 1], \text{ we have} $$

$$ \kappa^{\text{out}}(\sigma', \sigma) = [3 - 2\epsilon^2(\sigma)] \times \left[ \frac{\bar{\epsilon}(\sigma')}{\epsilon^3(\sigma')} \arcsin(\epsilon(\sigma')) - \frac{1}{\epsilon^2(\sigma')} \right] + 1. \quad (A2) $$

One could be worried by the multiple divergences in $\kappa^{\text{in}}$ at $\sigma = 0$. Yet, in this case, we see that the first integral in (16) vanishes and the divergences are then never taken into account in the calculations.

The property $\kappa^{\text{in}}(\sigma, \sigma) = \kappa^{\text{out}}(\sigma, \sigma)$ is easily proven by remembering $c(\sigma, \sigma) = 0$.

The derivative functions of $\kappa^{\text{in}}$ and $\kappa^{\text{out}}$ written in equations (22) and (23) are given by

$$ \chi(\sigma', \sigma) = \frac{\sigma'^2 \bar{\epsilon}(\sigma')}{\sigma^3 \epsilon(\sigma')} \left[ \arcsin \left( \frac{\sigma' \epsilon(\sigma')}{\sigma} \right) - \arcsin \left( \frac{\sigma' \epsilon(\sigma')}{\sigma \sqrt{1 + c(\sigma', \sigma)}} \right) \right], \quad (A3) $$

$$ \mu(\sigma', \sigma) = \frac{\bar{\epsilon}(\sigma')}{\epsilon^3(\sigma')} \left[ \arcsin \left( \frac{\sigma' \epsilon(\sigma')}{\sigma \sqrt{1 + c(\sigma', \sigma)}} \right) - \frac{\sigma' \epsilon(\sigma')}{\sigma \bar{\epsilon}(\sigma')} \right], \quad (A4) $$

$$ \nu(\sigma') = \frac{\bar{\epsilon}(\sigma')}{\epsilon^3(\sigma')} \arcsin(\epsilon(\sigma')) - \frac{1}{\epsilon^2(\sigma')}, \quad (A5) $$

where $\chi$ and $\mu$ are defined for $\sigma \in [0, 1]$, $\sigma' \in [0, \sigma[$ and $\nu$ is defined for $\sigma' \in [0, 1]$.

Once again, we can easily prove that $\chi(\sigma, \sigma) = 0$ and $\mu(\sigma, \sigma) = \nu(\sigma, \sigma)$, leading to the continuity in $\sigma' = \sigma$ of the derivative of the $\kappa$-functions.

APPENDIX B: ROBERTS’ EQUATION

Equation (3.23) of Roberts (1963) reads

$$ \begin{aligned} \frac{\mathrm{d}e^2}{\mathrm{d}\sigma} & \left\{ \frac{16\pi}{45} \hat{\Omega}^2 \sigma'^3 \bar{\epsilon}(\sigma) - \frac{4}{15} \hat{M}(\sigma) - \frac{2}{5} \hat{D}(\sigma) \right. \\ & \left. + \frac{3 - 2\epsilon^2(\sigma)}{2\epsilon^5(\sigma)} \bar{\epsilon}(\sigma) \left\{ [3 - 2\epsilon^2(\sigma)] \arcsin(\epsilon(\sigma)) - 3\epsilon(\sigma)\bar{\epsilon}(\sigma) \hat{D}(\sigma) \right\} \right\} \\ & = \frac{3 - 2\epsilon^2(\sigma)}{\sigma \epsilon^3(\sigma)} \bar{\epsilon}(\sigma) \left\{ [3 - 2\epsilon^2(\sigma)] \arcsin(\epsilon(\sigma)) - 3\epsilon(\sigma)\bar{\epsilon}(\sigma) \hat{D}(\sigma) \right\}, \end{aligned} \quad (B1) $$

where

$$ \hat{M}(\sigma) = 4\pi \int_{0}^{\sigma} d\sigma' \frac{\sigma' \hat{\rho}(\sigma')}{\bar{\epsilon}(\sigma')} \left[ 1 - \epsilon^2(\sigma) - \frac{\sigma'}{6} \frac{\mathrm{d}\epsilon^2}{\mathrm{d}\sigma'} \right] \quad (B2) $$

and

$$ \hat{D}(\sigma) = 4\pi \int_{0}^{\sigma} d\sigma' \frac{\sigma' \hat{\rho}(\sigma')}{\bar{\epsilon}(\sigma')} \left\{ \begin{aligned} & [1 - \epsilon^2(\sigma')] \left[ 1 - \frac{5}{3} \frac{\sigma'^2 \epsilon^2(\sigma')}{\sigma^2 \epsilon^2(\sigma)} \right] \\ & - \frac{\sigma'}{6} \frac{\mathrm{d}\epsilon^2}{\mathrm{d}\sigma'} \left[ 1 + \frac{\sigma'^2[2 - 3\epsilon^2(\sigma')]}{\sigma^2 \epsilon^2(\sigma)} \right] \end{aligned} \right\}. \quad (B3) $$
APPENDIX C: ON THE VOLUME INTEGRALS

The quantities calculated in Tables 1–4 write

$$ egin{cases} M = \int dV \rho(\boldsymbol{\omega}), \\ V = \int dV \left( = \frac{4}{3}\pi R_e^3 \bar{\epsilon}(1) \right), \\ J = \int dV \rho(\boldsymbol{\omega}) \Omega(\boldsymbol{\omega}) R^2, \\ W = \frac{1}{2} \int dV \rho(\boldsymbol{\omega}) \Psi(\boldsymbol{\omega}, \theta), \\ U = \frac{3}{2} \int dV p(\boldsymbol{\omega}), \\ T = \frac{1}{2} \int dV \rho(\boldsymbol{\omega}) \Omega^2(\boldsymbol{\omega}) R^2. \end{cases} \tag{C1} $$

So we need to express the volume element $dV$ and the cylindrical radius $R$ as functions of the spherical polar angle $\theta$ and $\varpi$. Along an isopycnic surface, the spherical radius $r$ reads

$$ \frac{r}{R_e} = \frac{\varpi \sqrt{1 - \epsilon^2(\varpi)}}{\sqrt{1 - \epsilon^2(\varpi) \sin^2(\theta)}}. \tag{C2} $$

The Jacobian matrix $\mathfrak{J}$ of the transformation from the Cartesian coordinates to an isopycnic coordinate system $(\varpi, \theta, \varphi)$, where $\varphi$ is the spherical azimuthal angle, reads

$$ \mathfrak{J} = \begin{pmatrix} \frac{\partial r}{\partial \varpi} & \sin(\theta) \cos(\varpi) & \frac{r \cos(\theta) \cos(\varpi)}{1 - \epsilon^2(\varpi) \sin^2(\theta)} & -r \sin(\theta) \sin(\varpi) \\ \frac{\partial r}{\partial \varpi} & \sin(\theta) \sin(\varpi) & \frac{r \cos(\theta) \sin(\varpi)}{1 - \epsilon^2(\varpi) \sin^2(\theta)} & r \sin(\theta) \cos(\varpi) \\ \frac{\partial r}{\partial \varpi} & \cos(\theta) & \frac{r \sin(\theta)}{1 - \epsilon^2(\varpi) \sin^2(\theta)} & 0 \end{pmatrix}, \tag{C3} $$

where $(\partial r / \partial \varpi)$ reads

$$ \frac{1}{R_e} \frac{\partial r}{\partial \varpi} = \frac{\sqrt{1 - \epsilon^2(\varpi)}}{\sqrt{1 - \epsilon^2(\varpi) \sin^2(\theta)}} - \frac{\varpi}{2} \frac{d\epsilon^2}{d\varpi} \frac{\cos^2(\theta)}{\sqrt{1 - \epsilon^2(\varpi)[1 - \epsilon^2(\varpi) \sin^2(\theta)]^{3/2}}}. \tag{C4} $$

So, the volume element is then given by

$$ \begin{aligned} \frac{dV}{R_e^3} &= \frac{\det(\mathfrak{J})}{R_e^3} d\varpi d\theta d\varphi \\ &= \frac{[1 - \epsilon^2(\varpi)][1 - \epsilon^2(\varpi) \sin^2(\theta)] - \frac{\varpi}{2} \frac{d\epsilon^2}{d\varpi} \cos^2(\theta)}{[1 - \epsilon^2(\varpi) \sin^2(\theta)]^{5/2}} \\ & \quad \times \varpi^2 \sqrt{1 - \epsilon^2(\varpi)} \sin(\theta) d\varpi d\theta d\varphi. \end{aligned} \tag{C5} $$

The cylindrical radius is given by $R = r \sin(\theta)$, namely

$$ \frac{R}{R_e} = \frac{\varpi \sqrt{1 - \epsilon^2(\varpi) \sin(\theta)}}{\sqrt{1 - \epsilon^2(\varpi) \sin^2(\theta)}}. \tag{C6} $$

All the integrals of (C1) can now be computed numerically (via a trapezoidal rule for instance).

For the rotation rate, as it is not exactly a constant due to the spheroidal approximation, we can obtain an mean value by integrating over the moment of inertia, namely

$$ \langle \Omega \rangle = \frac{\int dV \rho(\varpi) \Omega(\varpi) R^2}{\int dV \rho(\varpi) R^2} \cdot \left( = \frac{J}{I} \right). \tag{C7} $$

This paper has been typeset from a TeX/Latex file prepared by the author.