# Approaching the structure of rotating bodies from dimension reduction

C. Staelen, J.-M. Huré

► To cite this version:

C. Staelen, J.-M. Huré. Approaching the structure of rotating bodies from dimension reduction. Astronomy and Astrophysics - A&A, 2024, 684, pp.A59. 10.1051/0004-6361/202348590. hal-04558440

HAL Id: hal-04558440

https://hal.science/hal-04558440

Submitted on 24 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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.
# Approaching the structure of rotating bodies from dimension reduction

Clément Staelen, Jean-Marc Huré

► To cite this version:

Clément Staelen, Jean-Marc Huré. Approaching the structure of rotating bodies from dimension reduction. Astronomy & Astrophysics - A&A, 2024, 684, pp.A59. 10.1051/0004-6361/202348590. hal-04558440

HAL Id: hal-04558440

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

Submitted on 24 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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 Attribution 4.0 International License
# Approaching the structure of rotating bodies from dimension reduction ⋆

C. Staelen and J.-M. Huré

Univ. Bordeaux, CNRS, LAB, UMR 5804, 33600 Pessac, France
e-mail: clement.staelen@u-bordeaux.fr

Received 13 November 2023 / Accepted 25 January 2024

## ABSTRACT

We show that the two-dimensional structure of a rigidly rotating self-gravitating body is accessible with relatively good precision by assuming a purely spheroidal stratification. With this hypothesis, the two-dimensional problem becomes one-dimensional, and consists in solving two coupled fixed-point equations in terms of equatorial mass density and eccentricity of isopycnics. We propose a simple algorithm of resolution based on the self-consistent field method. Compared to the full unconstrained-surface two-dimensional problem, the precision in the normalized enthalpy field is better than $10^{-3}$ in absolute, and the computing time is drastically reduced. In addition, this one-dimensional approach is fully appropriate to fast rotators, works for any density profile (including any barotropic equation of state), and can account for mass density jumps in the system, including the existence of an ambient pressure. Several tests are given.

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

## 1. Introduction

The structure of static polytropic stars, in the classical sense, is traditionally described by the Lane-Emden equation, which admits a wide variety of solutions (e.g., Srivastava 1962; Sharma 1977; Seidov & Kuzakhmedov 1978; Liu 1996; Horedt 2004; Mach 2012; Tohline 2021). With rotation, the shape and the internal isobars deviate from sphericity, and it becomes difficult to anticipate precisely the topology of the gravitational field lines and to make direct use of the Gauss theorem. The context of slow rotation is attractive as it enables us to perform various kinds of analytical expansions, for instance in the form of series (Chandrasekhar 1933; Chandrasekhar & Roberts 1963; Roxburgh & Strittmatter 1966; Kovetz 1968). Unfortunately, the treatment required for moderate and fast rotations is much harder; there is almost no analytical way (Roberts 1963). The Poisson equation must be solved numerically, while the fluid boundary is not known in advance. Obviously, numerical methods offer a better range of options; for instance, they can model almost any kind of rotation profile, flattening, and equation of state (EOS) (e.g., James 1964; Ostriker & Mark 1968; Hachisu 1986a); they include magnetic fields (e.g., Tomimura & Eriguchi 2005; Lander & Jones 2009), ambient pressure (e.g., Huré et al. 2018), mass density jumps (Kiuchi et al. 2010; Kadam et al. 2016; Basillais & Huré 2021); and they can reach three-dimensional configurations and multiplicity (Hachisu 1986b; Even & Tohline 2009). In this last context, the self-consistent field (SCF) method of resolution, which consists in finding a fixed point for the mass density $ρ(r)$ from the pertinent equation set, is very efficient when appropriately initialized and scaled, and it has largely been used to model stars, binaries, and rings (Hachisu 1986a; Odrzywołek 2003). However, it is known from classical theories that, at slow rotations,

the equilibrium configurations remain very close to ellipsoids of revolution (e.g., Véronet 1912; Chandrasekhar 1933), and take a slightly sub-elliptical shape in a meridian plane between the pole and the equator (see also, e.g., Cisneros Parra et al. 2015). This assessment also applies for fast rotations, except close to the mass-shedding limit (Staelen & Huré 2024).

Any assumption made upon the mass density structure or on the symmetry is expected to reduce the mathematical complexity of the problem, but it diverts from the exact problem. This is the case when isobaric or isopycnic surfaces are locked to spheroids (i.e., ellipses in the meridional plane). Actually, under axial symmetry it is possible to benefit from Newton's and Maclaurin's theories and to use the closed-form for the potential of the homogeneous spheroid (Chandrasekhar 1969; Binney & Tremaine 1987), and subsequently construct heterogeneous bodies by piling-up coaxial homogeneous spheroids (see, e.g., Abramyan & Kaplan 1974; Montalvo et al. 1983). In a series of papers (Huré 2022a,b; Staelen & Huré 2024, hereafter Papers I, II, and IV, respectively), we developed a theory that solves the equilibrium of a heterogeneous body made of $\mathcal{L}$ homogeneous layers bounded by spheroids with different eccentricities, and in asynchronous rotational motion. As shown, it is possible to determine a relationship between the rotation rate $\Omega_i$ of each layer $i \in [1, \mathcal{L}]$, the parameter set of the spheroids $E_i$, and the mass-densities $\rho_i$ involved. This theory is approximate, but works very well provided the interfaces between layers are close enough to be confocal with each other. In the limit $\mathcal{L} \to \infty$, we showed that, for global rigid rotation, the eccentricity $\epsilon = [1 - b^2/a^2]^{1/2}$ of isopycnic surfaces (with semimajor axis $a$ and semiminor axis $b$) and the mass density profile $\rho$ obey a general integro-differential equation (IDE) (see Eq. (19) in Paper IV and Eq. (3) below). This IDE works very well for a wide range of rotation rates (or flattenings). It encompasses Clairaut's fundamental equation in the limit of small flattenings, and behaves correctly even close to the mass-shedding

<footnote>* A minimum driver program is available at https://github.com/clstaelen/ssba</footnote>
# Approaching the structure of rotating bodies from dimension reduction ⋆

C. Staelen and J.-M. Huré

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

e-mail: clement.staelen@u-bordeaux.fr

Received 13 November 2023 / Accepted 25 January 2024

## ABSTRACT

We show that the two-dimensional structure of a rigidly rotating self-gravitating body is accessible with relatively good precision by assuming a purely spheroidal stratification. With this hypothesis, the two-dimensional problem becomes one-dimensional, and consists in solving two coupled fixed-point equations in terms of equatorial mass density and eccentricity of isopycnics. We propose a simple algorithm of resolution based on the self-consistent field method. Compared to the full unconstrained-surface two-dimensional problem, the precision in the normalized enthalpy field is better than $10^{-3}$ in absolute, and the computing time is drastically reduced. In addition, this one-dimensional approach is fully appropriate to fast rotators, works for any density profile (including any barotropic equation of state), and can account for mass density jumps in the system, including the existence of an ambient pressure. Several tests are given.

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

## 1. Introduction

The structure of static polytropic stars, in the classical sense, is traditionally described by the Lane-Emden equation, which admits a wide variety of solutions (e.g., Srivastava 1962; Sharma 1977; Seidov & Kuzakhmedov 1978; Liu 1996; Horedt 2004; Mach 2012; Tohline 2021). With rotation, the shape and the internal isobars deviate from sphericity, and it becomes difficult to anticipate precisely the topology of the gravitational field lines and to make direct use of the Gauss theorem. The context of slow rotation is attractive as it enables us to perform various kinds of analytical expansions, for instance in the form of series (Chandrasekhar 1933; Chandrasekhar & Roberts 1963; Roxburgh & Strittmatter 1966; Kovetz 1968). Unfortunately, the treatment required for moderate and fast rotations is much harder; there is almost no analytical way (Roberts 1963). The Poisson equation must be solved numerically, while the fluid boundary is not known in advance. Obviously, numerical methods offer a better range of options; for instance, they can model almost any kind of rotation profile, flattening, and equation of state (EOS) (e.g., James 1964; Ostriker & Mark 1968; Hachisu 1986a); they include magnetic fields (e.g., Tomimura & Eriguchi 2005; Lander & Jones 2009), ambient pressure (e.g., Huré et al. 2018), mass density jumps (Kiuchi et al. 2010; Kadam et al. 2016; Basillais & Huré 2021); and they can reach three-dimensional configurations and multiplicity (Hachisu 1986b; Even & Tohline 2009). In this last context, the self-consistent field (SCF) method of resolution, which consists in finding a fixed point for the mass density $ρ(r)$ from the pertinent equation set, is very efficient when appropriately initialized and scaled, and it has largely been used to model stars, binaries, and rings (Hachisu 1986a; Odrzywołek 2003). However, it is known from classical theories that, at slow rotations,

the equilibrium configurations remain very close to ellipsoids of revolution (e.g., Véronet 1912; Chandrasekhar 1933), and take a slightly sub-elliptical shape in a meridian plane between the pole and the equator (see also, e.g., Cisneros Parra et al. 2015). This assessment also applies for fast rotations, except close to the mass-shedding limit (Staelen & Huré 2024).

Any assumption made upon the mass density structure or on the symmetry is expected to reduce the mathematical complexity of the problem, but it diverts from the exact problem. This is the case when isobaric or isopycnic surfaces are locked to spheroids (i.e., ellipses in the meridional plane). Actually, under axial symmetry it is possible to benefit from Newton's and Maclaurin's theories and to use the closed-form for the potential of the homogeneous spheroid (Chandrasekhar 1969; Binney & Tremaine 1987), and subsequently construct heterogeneous bodies by piling-up coaxial homogeneous spheroids (see, e.g., Abramyan & Kaplan 1974; Montalvo et al. 1983). In a series of papers (Huré 2022a,b; Staelen & Huré 2024, hereafter Papers I, II, and IV, respectively), we developed a theory that solves the equilibrium of a heterogeneous body made of $L$ homogeneous layers bounded by spheroids with different eccentricities, and in asynchronous rotational motion. As shown, it is possible to determine a relationship between the rotation rate $\Omega_i$ of each layer $i \in [1, L]$, the parameter set of the spheroids $E_i$, and the mass-densities $\rho_i$ involved. This theory is approximate, but works very well provided the interfaces between layers are close enough to be confocal with each other. In the limit $L \to \infty$, we showed that, for global rigid rotation, the eccentricity $\epsilon = [1 - b^2/a^2]^{1/2}$ of isopycnic surfaces (with semimajor axis $a$ and semiminor axis $b$) and the mass density profile $\rho$ obey a general integro-differential equation (IDE) (see Eq. (19) in Paper IV and Eq. (3) below). This IDE works very well for a wide range of rotation rates (or flattenings). It encompasses Clairaut's fundamental equation in the limit of small flattenings, and behaves correctly even close to the mass-shedding

<footnote>* A minimum driver program is available at https://github.com/clstaelen/ssba</footnote>
limit. This is important in order to provide tools appropriate to fast rotators, which may not be rare entities in the Universe (Rampalli et al. 2023).

The assumption of spheroidal isopycnics is undoubtedly strong, but it is not only motivated by mathematical simplifications coming out. It is widely supported by the numerical experiments. Actually, the true surface of a rigidly rotating fluid is generally very close to spheroids for a wide range of flattenings, even in the presence of mass density jumps. For instance, for a polytropic gas with index $n = 1$ and a polar-to-equatorial axis ratio $R_p/R_e = 0.95$, the shapes of isopycnics deviate in altitude from ellipses by less than $10^{-3}$ in relative from the center to the surface; the volumes differ by a similar value. For some problems, the deviations are too large and the hypothesis of spheroidal stratification must be abandoned (e.g., Zharkov & Trubitsyn 1970; Hubbard 2013; Cisneros Parra et al. 2015; Nettelmann 2017; Debras & Chabrier 2018). However, there are many situations where, on the contrary, such a precision is sufficient, for instance for the construction of mass-radius relationships or for dynamical studies (e.g., Hadjifotinou 2000; Kluźniak & Rosińska 2013; Mishra & Vaidya 2015; Venditti et al. 2020).

In this article we propose a direct exploitation of the IDE derived in Paper IV, looking for self-consistent solutions. As the underlying equation set is reduced to the Bernoulli equation combined with the Poisson equation, the solution necessarily has a limited range of applications, as quoted above. Our approach cannot compete with sophisticated models for the stellar structures for instance (see, e.g., Rieutord et al. 2006, 2016; Espinosa Lara & Rieutord 2013; Houdayer & Reese 2023, for the state of art). Our main aim is to analyze the performance of the spheroidally stratified barotrope (SSB) approximation, which to our knowledge has never been reported. By construction, the IDE concerns only equatorial values. When it is combined with the centrifugal equilibrium along the polar axis, we obtain a type of equatorial projection of the original bidimensional equation set. This opens the possibility to determine the two-dimensional structure of a fully heterogeneous rigidly rotating body from a one-dimensional approach. This is explained in Sect. 2. An attractive point is that the IDE is valuable regardless of the EOS, leaving a certain flexibility in terms of applications. We discuss in Sect. 3 a simple iterative algorithm to solve this set of projected equations. The cycle is based on the classical self-consistent field (SCF) method. In the present case, we have to solve two coupled fixed-point problems, one for the eccentricity profile $\epsilon(a)$ and one for the mass density profile $\rho(a)$, in the range $a \in [0, R_e]$. The full 2D structure is reconstructed by unfolding the equatorial solution $\{\rho(a), \epsilon(a)\}$. Several tests are proposed in Sect. 4, including static, slowly rotating, and highly rotating configurations, as well as Hachisu's equilibrium sequences (Hachisu 1986a). For this purpose, a numerical reference is needed to compare the approximate solution to the real one. We use the DROP code, which computes the equilibrium configurations in full 2D (Huré & Hersant 2017; Basillais & Huré 2021), the Poisson being solved by the multigrid-method from finite-difference equation. Using the spectral version of the code would not help much as we are dealing with deviations that are not tiny, typically on the order of $10^{-4}$. In Sect. 5 we show how to account for mass density jumps, which enables us to model multi-domain bodies. Two examples are discussed, namely the pressurized $n = 5$ polytrope and the Earth's interior, by using the nonrotating Preliminary Reference Earth Model (Dziewonski & Anderson 1981). A summary is found in the last section.

## 2. Theoretical background

Basically, a star is an equilibrium state between gas pressure, gravity forces, and centrifugation. In the barotropic approximation, the relevant equation set is

$$ H + \Phi + \Psi = C^{\text{te}}, \qquad (1a) $$

$$ \nabla^2 \Psi = 4\pi G\rho, \qquad (1b) $$

$$ f(H, \rho) = 0, \qquad (1c) $$

where $\Psi$ is the gravitational potential, $\Phi = -\int dR \Omega^2 R$ is the centrifugal potential, $\Omega$ is the local rotation rate, $P$ is the pressure, $\rho$ is mass density of the fluid, and $H = \int dP/\rho$ is the enthalpy (the self-gravitating flow is isentropic). The link between $H$ and $\rho$ in Eq. (1c) is usually done via the EOS, namely $P(\rho)$. This set represents a conservative form for the flow, is independent of time, ignores viscosity, and supposes that rotation is constant along cylinders, coaxial with the rotation axis (e.g., Tassoul 1978; Amendt et al. 1989), which means that $\Phi$ depends on the cylindrical radius $R$ only. In the present article we consider axially symmetric configurations in the framework of the nested spheroidal figure of equilibrium (NSFoE) reported in Papers I and II, which assumes that, in a layered system, all the surfaces bounding the homogeneous layers ($\mathcal{L}$ in total, $\rho_i$ is the mass density of layer $i$) are perfect spheroids $E_i(a_i, b_i)$. In the meridional plane, these surfaces are ellipses with equation

$$ \frac{R^2}{a_i^2} + \frac{Z^2}{b_i^2} = 1, \qquad (2) $$

where $a_i \in [0, R_e]$, $b_i \in [0, R_p]$ and $R_e$ and $R_p$ are the equatorial and polar radii of the body, respectively. In these conditions, the equilibrium (if it exists) is perfectly determined as a function of the set $\{\rho_i, E_i\}_{i \in [1, \mathcal{L}]}$.

### 2.1. Dimension reduction (equatorial projection)

In Paper IV we showed that, in the limit where the number $\mathcal{L}$ of layers is infinite, Eqs. (1a) and (1b) can be formally solved regardless of any EOS. In the special case of rigid rotation considered here, we obtain an IDE linking the eccentricity profile $\epsilon(a)$ and the mass density profile $\rho(a)$. In order to render the problem scale-free, the semimajor axis $a$ of a given isopycnic surface is expressed in units of the equatorial radius $\bar{\omega} = a/R_e \in [0, 1]$, and $\hat{\rho} = \rho/\rho_c$ is the dimensionless mass density ($\rho_c$ is the central mass density). In these conditions, the IDE reads

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

where $\chi$ and $\mu$ are continuous smooth functions in the whole domain $\bar{\omega} \in [0, 1]$. These are defined in Paper IV (see also Appendix A, respectively Eqs. (A.3) and (A.4)). As explicitly stated, $\chi$ and $\mu$ depend on the radius $\bar{\omega}$ and on the eccentricity profile $\epsilon(\bar{\omega})$. We can safely replace Eqs. (1a) and (1b) by Eq. (3), and express Eq. (1a) along the polar axis. This is possible as the gravitational potential is known. As the constant in the RHS, we can take the value of the LHS at the pole, which is the most straightforward. After some algebra (see Appendix C for a proof), we find for the dimensionless enthalpy

$$ \hat{H}(\bar{\omega}) = \hat{H}(1) + 2\pi \int_{\hat{\rho}(0)}^{\hat{\rho}(1)} d\hat{\rho}(\bar{\omega}') [\eta(\bar{\omega}', 1; \epsilon) - \eta(\bar{\omega}', \bar{\omega}; \epsilon)], $$
limit. This is important in order to provide tools appropriate to fast rotators, which may not be rare entities in the Universe (Rampalli et al. 2023).

The assumption of spheroidal isopycnics is undoubtedly strong, but it is not only motivated by mathematical simplifications coming out. It is widely supported by the numerical experiments. Actually, the true surface of a rigidly rotating fluid is generally very close to spheroids for a wide range of flattenings, even in the presence of mass density jumps. For instance, for a polytropic gas with index $n = 1$ and a polar-to-equatorial axis ratio $R_p/R_e = 0.95$, the shapes of isopycnics deviate in altitude from ellipses by less than $10^{-3}$ in relative from the center to the surface; the volumes differ by a similar value. For some problems, the deviations are too large and the hypothesis of spheroidal stratification must be abandoned (e.g., Zharkov & Trubitsyn 1970; Hubbard 2013; Cisneros Parra et al. 2015; Nettelmann 2017; Debras & Chabrier 2018). However, there are many situations where, on the contrary, such a precision is sufficient, for instance for the construction of mass-radius relationships or for dynamical studies (e.g., Hadjifotinou 2000; Kluźniak & Rosińska 2013; Mishra & Vaidya 2015; Venditti et al. 2020).

In this article we propose a direct exploitation of the IDE derived in Paper IV, looking for self-consistent solutions. As the underlying equation set is reduced to the Bernoulli equation combined with the Poisson equation, the solution necessarily has a limited range of applications, as quoted above. Our approach cannot compete with sophisticated models for the stellar structures for instance (see, e.g., Rieutord et al. 2006, 2016; Espinosa Lara & Rieutord 2013; Houdayer & Reese 2023, for the state of art). Our main aim is to analyze the performance of the spheroidally stratified barotrope (SSB) approximation, which to our knowledge has never been reported. By construction, the IDE concerns only equatorial values. When it is combined with the centrifugal equilibrium along the polar axis, we obtain a type of equatorial projection of the original bidimensional equation set. This opens the possibility to determine the two-dimensional structure of a fully heterogeneous rigidly rotating body from a one-dimensional approach. This is explained in Sect. 2. An attractive point is that the IDE is valuable regardless of the EOS, leaving a certain flexibility in terms of applications. We discuss in Sect. 3 a simple iterative algorithm to solve this set of projected equations. The cycle is based on the classical self-consistent field (SCF) method. In the present case, we have to solve two coupled fixed-point problems, one for the eccentricity profile $\varepsilon(a)$ and one for the mass density profile $\rho(a)$, in the range $a \in [0, R_e]$. The full 2D structure is reconstructed by unfolding the equatorial solution $\{\rho(a), \varepsilon(a)\}$. Several tests are proposed in Sect. 4, including static, slowly rotating, and highly rotating configurations, as well as Hachisu's equilibrium sequences (Hachisu 1986a). For this purpose, a numerical reference is needed to compare the approximate solution to the real one. We use the DROP code, which computes the equilibrium configurations in full 2D (Huré & Hersant 2017; Basillais & Huré 2021), the Poisson being solved by the multigrid-method from finite-difference equation. Using the spectral version of the code would not help much as we are dealing with deviations that are not tiny, typically on the order of $10^{-4}$. In Sect. 5 we show how to account for mass density jumps, which enables us to model multi-domain bodies. Two examples are discussed, namely the pressurized $n = 5$ polytrope and the Earth's interior, by using the nonrotating Preliminary Reference Earth Model (Dziewonski & Anderson 1981). A summary is found in the last section.

## 2. Theoretical background

Basically, a star is an equilibrium state between gas pressure, gravity forces, and centrifugation. In the barotropic approximation, the relevant equation set is

$$ \left\{ \begin{array}{l} H + \Phi + \Psi = C^{\text{te}}, \\ \nabla^2 \Psi = 4\pi G \rho, \\ f(H, \rho) = 0, \end{array} \right. \qquad (1a) $$

$$ \nabla^2 \Psi = 4\pi G \rho, \qquad (1b) $$

$$ f(H, \rho) = 0, \qquad (1c) $$

where $\Psi$ is the gravitational potential, $\Phi = -\int dR \Omega^2 R$ is the centrifugal potential, $\Omega$ is the local rotation rate, $P$ is the pressure, $\rho$ is mass density of the fluid, and $H = \int dP/\rho$ is the enthalpy (the self-gravitating flow is isentropic). The link between $H$ and $\rho$ in Eq. (1c) is usually done via the EOS, namely $P(\rho)$. This set represents a conservative form for the flow, is independent of time, ignores viscosity, and supposes that rotation is constant along cylinders, coaxial with the rotation axis (e.g., Tassoul 1978; Amendt et al. 1989), which means that $\Phi$ depends on the cylindrical radius $R$ only. In the present article we consider axially symmetric configurations in the framework of the nested spheroidal figure of equilibrium (NSFoE) reported in Papers I and II, which assumes that, in a layered system, all the surfaces bounding the homogeneous layers ($\mathcal{L}$ in total, $\rho_i$ is the mass density of layer $i$) are perfect spheroids $E_i(a_i, b_i)$. In the meridional plane, these surfaces are ellipses with equation

$$ \frac{R^2}{a_i^2} + \frac{Z^2}{b_i^2} = 1, \qquad (2) $$

where $a_i \in [0, R_e]$, $b_i \in [0, R_p]$ and $R_e$ and $R_p$ are the equatorial and polar radii of the body, respectively. In these conditions, the equilibrium (if it exists) is perfectly determined as a function of the set $\{\rho_i, E_i\}_{i \in [1, \mathcal{L}]}$.

### 2.1. Dimension reduction (equatorial projection)

In Paper IV we showed that, in the limit where the number $\mathcal{L}$ of layers is infinite, Eqs. (1a) and (1b) can be formally solved regardless of any EOS. In the special case of rigid rotation considered here, we obtain an IDE linking the eccentricity profile $\varepsilon(a)$ and the mass density profile $\rho(a)$. In order to render the problem scale-free, the semimajor axis $a$ of a given isopycnic surface is expressed in units of the equatorial radius $\bar{\omega} = a/R_e \in [0, 1]$, and $\hat{\rho} = \rho/\rho_c$ is the dimensionless mass density ($\rho_c$ is the central mass density). In these conditions, the IDE reads

$$ \frac{d\varepsilon^2}{d\bar{\omega}} = \frac{2 \int_{\hat{\rho}(0)}^{\hat{\rho}(\bar{\omega})} d\hat{\rho}(\bar{\omega}') \chi(\bar{\omega}', \bar{\omega}; \varepsilon)}{\int_{\hat{\rho}(0)}^{\hat{\rho}(1)} d\hat{\rho}(\bar{\omega}') \mu(\bar{\omega}', \bar{\omega}; \varepsilon)}, \qquad (3) $$

where $\chi$ and $\mu$ are continuous smooth functions in the whole domain $\bar{\omega} \in [0, 1]$. These are defined in Paper IV (see also Appendix A, respectively Eqs. (A.3) and (A.4)). As explicitly stated, $\chi$ and $\mu$ depend on the radius $\bar{\omega}$ and on the eccentricity profile $\varepsilon(\bar{\omega})$. We can safely replace Eqs. (1a) and (1b) by Eq. (3), and express Eq. (1a) along the polar axis. This is possible as the gravitational potential is known. As the constant in the RHS, we can take the value of the LHS at the pole, which is the most straightforward. After some algebra (see Appendix C for a proof), we find for the dimensionless enthalpy

$$ \hat{H}(\bar{\omega}) = \hat{H}(1) + 2\pi \int_{\hat{\rho}(0)}^{\hat{\rho}(1)} d\hat{\rho}(\bar{\omega}') [\eta(\bar{\omega}', 1; \varepsilon) - \eta(\bar{\omega}', \bar{\omega}; \varepsilon)], $$
<figure><img src="image_1.png" /><figcaption>**Fig. 1.** Schematic of the unfolding process. The one-dimensional equation set is first solved for the mass density $\rho$ and eccentricity $\varepsilon$ along the equatorial axis in the full interval $[0, R_e]$. Then, the two-dimensional structure is reconstructed by propagating, for each radius $(a, 0)$, the solution (enthalpy, mass density, or pressure) along the ellipse $E(a, b)$ up to the polar axis at $(0, b)$ (see Sect. 2.2).</figcaption></figure>

where $\hat{H} = H/G\rho_c R_e^2$, and $\eta$ is reported in Appendix A (see Eq. (A.5)). It follows that Eq. (1), in the framework of the theory of NSFoE, becomes

$$ \left\{ \begin{array}{l} \hat{H}(\omega) = \hat{H}(1) + \\ \quad 2\pi \int_{\hat{\rho}(0)}^{\hat{\rho}(1)} d\hat{\rho}(\omega') [\eta(\omega', 1; \varepsilon) - \eta(\omega', \omega; \varepsilon)] \end{array} \right. \qquad (4a) $$

$$ \frac{d\varepsilon^2}{d\omega} = \frac{2 \int_{\hat{\rho}(0)}^{\hat{\rho}(\omega)} d\hat{\rho}(\omega') \chi(\omega', \omega; \varepsilon)}{\int_{\hat{\rho}(0)}^{\hat{\rho}(1)} d\hat{\rho}(\omega') \mu(\omega', \omega; \varepsilon)}, \qquad (4b) $$

$$ f(\hat{H}(\omega), \hat{\rho}(\omega)) = 0, \qquad (4c) $$

where the last equation is the EOS. We note that, if the mass density profile happens to be known in advance (e.g., deduced from observational data) or prescribed, Eq. (4a) is not needed anymore (see, e.g., Sect. 5.3). We immediately see that the problem is now fully one-dimensional. It means that we can compute the 2D structure of a rotating gaseous body from a 1D approach. This dimension reduction thus relies on a SSB approximation.

## 2.2. Unfolding and global quantities

Once the 1D solution {$\hat{\rho}(\omega), \varepsilon(\omega)$} is calculated, we can reconstruct any bidimensional map by unfolding the mass density, at each radius $a$, along the corresponding ellipse $E(a, b)$ up to the polar axis at $(0, b)$, as depicted in Fig. 1. In spherical polar coordinates $(r, \theta, \varphi)$, we actually have

$$ \hat{\rho}(r, \theta) = \hat{\rho}(\omega), \qquad (5) $$

where the two coordinates are easily determined by using the parametric formula

$$ r^2 = R_e^2 \omega^2 [1 - \varepsilon^2(\omega) \sin^2(\beta)], \qquad (6a) $$

$$ \tan(\theta) = \frac{1}{\bar{\varepsilon}(\omega)} \cot(\beta), \qquad (6b) $$

with $\theta \in [0, \pi/2]$ and $\beta \in [0, \pi/2]$, and where

$$ \bar{\varepsilon}(\omega) = [1 - \varepsilon^2(\omega)]^{1/2} \qquad (7) $$

is the axis ratio of the isopycnic $E(a, b)$ considered. The conversion to cylindrical coordinates is straightforward:

$R = R_e \omega \cos(\beta)$ and $Z = R_e \omega \bar{\varepsilon}(\omega) \sin(\beta)$. We can subsequently deduce all global quantities (e.g., mass, volume).

In $(\omega, \beta, \varphi)$ coordinates, the volume element writes

$$ \frac{dV}{R_e^3} = 2\pi \frac{\omega^2 \cos(\beta)}{\bar{\varepsilon}(\omega)} \left[ \bar{\varepsilon}^2(\omega) - \frac{1}{2} \sin^2(\beta) \frac{d\varepsilon^2}{d\omega} \right] d\omega d\beta d\varphi. \qquad (8) $$

For most global properties of the system, the integrations over $\beta$ are trivial. For instance, the mass $M$ and the moment of inertia $I$ are reduced to a single integral over $\omega$ and we have

$$ \left\{ \begin{array}{l} \frac{M}{\rho_c R_e^3} = 4\pi \int_0^1 d\omega \frac{\omega^2 \hat{\rho}(\omega)}{\bar{\varepsilon}(\omega)} \left[ \bar{\varepsilon}^2(\omega) - \frac{1}{6} \frac{d\varepsilon^2}{d\omega} \right], \\ \frac{I}{\rho_c R_e^5} = \frac{8}{3\pi} \int_0^1 d\omega \frac{\omega^4 \hat{\rho}(\omega)}{\bar{\varepsilon}(\omega)} \left[ \bar{\varepsilon}^2(\omega) - \frac{1}{10} \frac{d\varepsilon^2}{d\omega} \right]. \end{array} \right. \qquad (9) $$

According to Paper IV, the rotation rate $\Omega$ is given by

$$ \hat{\Omega}^2 = -2\pi \int_{\hat{\rho}(0)}^{\hat{\rho}(1)} d\hat{\rho}(\omega') \kappa(\omega', \omega; \varepsilon), \qquad (10) $$

where $\hat{\Omega} = \Omega/\sqrt{G\rho_c}$ and $\kappa$ is defined by Eq. (A.6) (see again Paper IV). It is worth recalling that, as the theory of NSFoE is an approximate theory, Eq. (10) is expected to exhibit slight variations with the radius (see below), and should rigorously be denoted $\hat{\Omega}^2(\omega)$.

## 2.3. Comments

As can be seen in Paper IV, the three functions $\chi$, $\mu$, and $\eta$ have small amplitude, but all take real, negative values in the range of $\omega$ of interest. It is therefore clear that, if the integral in Eq. (4a) happens to be essentially negative, then the mass density can become negative, in which case the computed solution cannot be compatible with a physical solution. Density inversions are eventually acceptable, but $\rho \ge 0$ is a firm condition. A similar remark holds for Eq. (10). For Eq. (4b), things are less restrictive. Actually, negative values of $\varepsilon^2$ correspond to prolate spheroids. This is not a problem because there is a mathematical continuity in the gravitational potential when the eccentricity becomes a purely imaginary number (i.e., when $\varepsilon = 0^+ \rightarrow i0^+$; Binney & Tremaine 1987). From a numerical point of view, however, the requires a dedicated treatment, and it seems preferable to consider the axis ratio, namely $\bar{\varepsilon}(\omega) \le 1$ instead of $\varepsilon$ (see below, and Appendix B).

# 3. Solution with a SCF-algorithm. Implementation and example

## 3.1. Cycle and convergence criterion

It is well known that Eqs. (1a) and (1b) correspond to a fixed-point problem: $\rho = f(\rho)$ in terms of the mass density profile. Clearly, the equation set Eq. (4) has a similar structure, but we have two coupled fixed-point problems of the form

$$ \hat{\rho}^{(t)} \leftarrow f_1(\hat{\rho}^{(t-1)}, \bar{\varepsilon}^{(t-1)}), \qquad (11a) $$

$$ \bar{\varepsilon}^{(t)} \leftarrow f_2(\hat{\rho}^{(t-1)}, \bar{\varepsilon}^{(t-1)}), \qquad (11b) $$

It is therefore natural to proceed in the same manner as for the standard SCF method: we guess the profiles for the mass density and the axis ratio, namely $\hat{\rho}^{(0)}(\omega)$ and $\bar{\varepsilon}^{(0)}(\omega)$, and let the profiles
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Schematic of the unfolding process. The one-dimensional equation set is first solved for the mass density $\rho$ and eccentricity $\varepsilon$ along the equatorial axis in the full interval $[0, R_e]$. Then, the two-dimensional structure is reconstructed by propagating, for each radius $(a, 0)$, the solution (enthalpy, mass density, or pressure) along the ellipse $E(a, b)$ up to the polar axis at $(0, b)$ (see Sect. 2.2).</figcaption></figure>

where $\hat{H} = H/G\rho_c R_e^2$, and $\eta$ is reported in Appendix A (see Eq. (A.5)). It follows that Eq. (1), in the framework of the theory of NSFoE, becomes

$$ \left\{ \begin{aligned} \hat{H}(\omega) = \hat{H}(1) + & \quad (4a) \\ 2\pi \int_{\hat{\rho}(0)}^{\hat{\rho}(1)} d\hat{\rho}(\omega') [\eta(\omega', 1; \varepsilon) - \eta(\omega', \omega; \varepsilon)] \end{aligned} \right. $$

$$ \frac{d\varepsilon^2}{d\omega} = \frac{2 \int_{\hat{\rho}(0)}^{\hat{\rho}(\omega)} d\hat{\rho}(\omega') \chi(\omega', \omega; \varepsilon)}{\int_{\hat{\rho}(0)}^{\hat{\rho}(1)} d\hat{\rho}(\omega') \mu(\omega', \omega; \varepsilon)}, \quad (4b) $$

$$ f(\hat{H}(\omega), \hat{\rho}(\omega)) = 0, \quad (4c) $$

where the last equation is the EOS. We note that, if the mass density profile happens to be known in advance (e.g., deduced from observational data) or prescribed, Eq. (4a) is not needed anymore (see, e.g., Sect. 5.3). We immediately see that the problem is now fully one-dimensional. It means that we can compute the 2D structure of a rotating gaseous body from a 1D approach. This dimension reduction thus relies on a SSB approximation.

## 2.2. Unfolding and global quantities

Once the 1D solution {$\hat{\rho}(\omega), \varepsilon(\omega)$} is calculated, we can reconstruct any bidimensional map by unfolding the mass density, at each radius $a$, along the corresponding ellipse $E(a, b)$ up to the polar axis at $(0, b)$, as depicted in Fig. 1. In spherical polar coordinates $(r, \theta, \varphi)$, we actually have

$$ \hat{\rho}(r, \theta) = \hat{\rho}(\omega), \quad (5) $$

where the two coordinates are easily determined by using the parametric formula

$$ r^2 = R_e^2 \omega^2 [1 - \varepsilon^2(\omega) \sin^2(\beta)], \quad (6a) $$

$$ \tan(\theta) = \frac{1}{\bar{\varepsilon}(\omega)} \cot(\beta), \quad (6b) $$

with $\theta \in [0, \pi/2]$ and $\beta \in [0, \pi/2]$, and where

$$ \bar{\varepsilon}(\omega) = [1 - \varepsilon^2(\omega)]^{1/2} \quad (7) $$

is the axis ratio of the isopycnic $E(a, b)$ considered. The conversion to cylindrical coordinates is straightforward:

$R = R_e \omega \cos(\beta)$ and $Z = R_e \omega \bar{\varepsilon}(\omega) \sin(\beta)$. We can subsequently deduce all global quantities (e.g., mass, volume).

In $(\omega, \beta, \varphi)$ coordinates, the volume element writes

$$ \frac{dV}{R_e^3} = 2\pi \frac{\omega^2 \cos(\beta)}{\bar{\varepsilon}(\omega)} \left[ \bar{\varepsilon}^2(\omega) - \frac{1}{2} \sin^2(\beta) \frac{d\varepsilon^2}{d\omega} \right] d\omega d\beta d\varphi. \quad (8) $$

For most global properties of the system, the integrations over $\beta$ are trivial. For instance, the mass $M$ and the moment of inertia $I$ are reduced to a single integral over $\omega$ and we have

$$ \left\{ \begin{aligned} \frac{M}{\rho_c R_e^3} &= 4\pi \int_0^1 d\omega \frac{\omega^2 \hat{\rho}(\omega)}{\bar{\varepsilon}(\omega)} \left[ \bar{\varepsilon}^2(\omega) - \frac{1}{6} \frac{d\varepsilon^2}{d\omega} \right], \\ \frac{I}{\rho_c R_e^5} &= \frac{8}{3\pi} \int_0^1 d\omega \frac{\omega^4 \hat{\rho}(\omega)}{\bar{\varepsilon}(\omega)} \left[ \bar{\varepsilon}^2(\omega) - \frac{1}{10} \frac{d\varepsilon^2}{d\omega} \right]. \end{aligned} \right. \quad (9) $$

According to Paper IV, the rotation rate $\Omega$ is given by

$$ \hat{\Omega}^2 = -2\pi \int_{\hat{\rho}(0)}^{\hat{\rho}(1)} d\hat{\rho}(\omega') \kappa(\omega', \omega; \varepsilon), \quad (10) $$

where $\hat{\Omega} = \Omega / \sqrt{G\rho_c}$ and $\kappa$ is defined by Eq. (A.6) (see again Paper IV). It is worth recalling that, as the theory of NSFoE is an approximate theory, Eq. (10) is expected to exhibit slight variations with the radius (see below), and should rigorously be denoted $\hat{\Omega}^2(\omega)$.

## 2.3. Comments

As can be seen in Paper IV, the three functions $\chi$, $\mu$, and $\eta$ have small amplitude, but all take real, negative values in the range of $\omega$ of interest. It is therefore clear that, if the integral in Eq. (4a) happens to be essentially negative, then the mass density can become negative, in which case the computed solution cannot be compatible with a physical solution. Density inversions are eventually acceptable, but $\rho \ge 0$ is a firm condition. A similar remark holds for Eq. (10). For Eq. (4b), things are less restrictive. Actually, negative values of $\varepsilon^2$ correspond to prolate spheroids. This is not a problem because there is a mathematical continuity in the gravitational potential when the eccentricity becomes a purely imaginary number (i.e., when $\varepsilon = 0^+ \rightarrow i0^+$; Binney & Tremaine 1987). From a numerical point of view, however, the requires a dedicated treatment, and it seems preferable to consider the axis ratio, namely $\bar{\varepsilon}(\omega) \le 1$ instead of $\varepsilon$ (see below, and Appendix B).

# 3. Solution with a SCF-algorithm. Implementation and example

## 3.1. Cycle and convergence criterion

It is well known that Eqs. (1a) and (1b) correspond to a fixed-point problem: $\rho = f(\rho)$ in terms of the mass density profile. Clearly, the equation set Eq. (4) has a similar structure, but we have two coupled fixed-point problems of the form

$$ \hat{\rho}^{(t)} \leftarrow f_1(\hat{\rho}^{(t-1)}, \bar{\varepsilon}^{(t-1)}), \quad (11a) $$

$$ \bar{\varepsilon}^{(t)} \leftarrow f_2(\hat{\rho}^{(t-1)}, \bar{\varepsilon}^{(t-1)}), \quad (11b) $$

It is therefore natural to proceed in the same manner as for the standard SCF method: we guess the profiles for the mass density and the axis ratio, namely $\hat{\rho}^{(0)}(\omega)$ and $\bar{\varepsilon}^{(0)}(\omega)$, and let the profiles
evolve until stabilization. This can be accomplished following the iterative scheme, for $t \ge 1$:

1. $\bar{\epsilon}^{(t)}$ is obtained from Eq. (7), after integrating Eq. (4b);

2. $\hat{H}^{(t)}$ is computed according to Eq. (4a);

3. $\hat{\rho}^{(t)}$ is obtained from Eq. (4c).

We see that there are two other options (with no significant impact on the performance of the cycle), depending on the order of assignments; in other words, $\bar{\epsilon}$ in Eq. (11b) can be computed using the $\rho$-profile updated from Eq. (11a), or $\rho$ in Eq. (11a) can be computed from the $\bar{\epsilon}$-profile updated from Eq. (11b). At convergence, two successive profiles must be numerically unchanged, in which case the cycle ends. As we traditionally use double precision computers, we take$^{1,2}$ $\delta^{(t)} \le 10^{-14}$, where

$$ \delta^{(t)} = \max(\Delta^{(t)}\hat{\rho}, \Delta^{(t)}\bar{\epsilon}), \quad (12) $$

where we define still for $t \ge 1$ and $\bar{\omega} \in [0, 1]$

$$ \Delta^{(t)}\hat{\rho} = \sup |\hat{\rho}^{(t)}(\bar{\omega}) - \hat{\rho}^{(t-1)}(\bar{\omega})|, \quad (13) $$

and use a similar definition for $\Delta^{(t)}\bar{\epsilon}$.

## 3.2. Boundary conditions

In the standard (single-domain) case, and in the absence of external pressure (see Sect. 5.2), we can take Dirichlet boundary conditions (BCs) at the outer boundary $\bar{\omega} = 1$

$$ \begin{cases} \hat{\rho}(1) = 0 \\ \bar{\epsilon}(1) = \bar{\epsilon}_s \end{cases}, \quad (14) $$

where $\bar{\epsilon}_s = R_p/R_c$ is the axis ratio of the outermost surface. At the center of the body, we have

$$ \left\{ \begin{array}{l} \hat{\rho}(0) = 1 \\ \left. \frac{d\bar{\epsilon}^2}{d\bar{\omega}} \right|_{\bar{\omega}=0} = 0. \end{array} \right. \quad (15) $$

Unfortunately, $\bar{\epsilon}(0)$ cannot be easily deduced.

## 3.3. Method and implementation (quadrature schemes, $\bar{\omega}$-grid, and seeds)

We see that Eq. (4) involves derivatives and quadratures, and there are different ways to estimate $\rho(\bar{\omega})$ and $\bar{\epsilon}(\bar{\omega})$ for given BCs. Here we decided to recast Eq. (4b) in integral form, i.e.,

$$ \bar{\epsilon}^2(\bar{\omega}) = \bar{\epsilon}_s^2 - \int_{\bar{\omega}}^{1} \left( \frac{d\bar{\epsilon}^2}{d\bar{\omega}} \right) d\bar{\omega}, \quad (16) $$

where the term in parentheses is simply the right-hand side of Eq. (4b) within a minus sign. So, we have only to deal with quadratures. The eccentricity being unknown at the center, we chose to integrate downward, from the surface to the center, with Eq. (14) as Dirichlet BCs. The computational grid is made of $N+1$ equally spaced nodes: $\bar{\omega}_i \in \{0, \frac{1}{N}, \dots, \bar{\omega}_i = \frac{i}{N}, \dots, 1\}$. We used the trapezoidal rule as the quadrature rule, which is second-order in the step size (more efficient schemes can be used

at this level). According to Eq. (14), still assuming null ambient pressure, we take

$$ \left\{ \begin{array}{l} \hat{H}^{(0)}(\bar{\omega}) = \hat{H}_c(1 - \bar{\omega}^2), \\ \bar{\epsilon}^{(0)}(\bar{\omega}) = 1 - (1 - \bar{\epsilon}_s)\bar{\omega}^2 \end{array} \right. \quad (17a) $$

$$ \quad (17b) $$

as the starting guess. Initially, we thus have $\hat{H}^{(0)}(0) = \hat{H}_c$ at the center, and $\bar{\epsilon}^{(0)}(0) = 1$. These quadratic profiles seem appropriate for a wide range of configurations, although the solutions generally have a nonzero central eccentricity. Obviously, the BCs must be applied at each step in the cycle.

## 3.4. A note about the equation of state

The EOS is the fundamental ingredient. Without loss of generality, we consider a polytropic gas where the pressure is a power law of the mass density, which leads to the expression $H = K(n+1)\rho^{1/n}$, where K is a positive constant and the polytropic index $n > 0$ is finite. In this case the relationship between $\hat{H}$ and $\hat{\rho}$ is

$$ \frac{\hat{H}}{\hat{H}_c} = \hat{\rho}^{1/n}, \quad (18) $$

where $\hat{H}_c = K(n+1)/G\rho_c^{1-1/n}R_c^2$ (this relation assumes $\hat{\rho} \ge 0$). The mass density $\hat{\rho}^{(t)}$ inside the cycle, including the seed, is deduced from the EOS through Eq. (4c).

## 3.5. Example of cycle convergence

The first example is for $\bar{\epsilon}_s = 0.9$ and $n = 3$, and the grid has $N+1 = 257$ equally spaced nodes. This configuration corresponds, for instance, to a radiation-pressure-dominated ideal gas or to a white dwarf with fully degenerate extremely relativistic electrons. With this parameter, we are already beyond slow rotators. We ran the code based on the SSB approximation. Figure 2 shows the evolution of $\hat{\rho}(\bar{\omega})$ and $\bar{\epsilon}^2(\bar{\omega})$ during the cycle, as well as the deviations $\Delta\hat{\rho}$ and $\Delta\bar{\epsilon}^2$ from one step to the other. Figure 2f gives the convergence parameter $\delta^{(t)}$ defined by Eq. (12) from the beginning to convergence. We see that $\delta^{(t)}$ decreases exponentially with the step $t$ and the algorithm converges quickly on a solution. Convergence is reached after 70 cycles with the current criterion. After step 20, we already have $\delta^{(t)} \sim 10^{-5}$, which is on the order of the accuracy of the numerical scheme (i.e., $1/N^2 \sim 10^{-5}$ with 257 nodes). The next steps are necessary to reach the threshold of $10^{-14}$. We show in Fig. 2c the rotation rate computed from Eq. (10). Unsurprisingly, there is a slight variation with the radius, on the order of $10^{-2}$ in relative. This is due to the approximate nature of the theory of NSFoE. The function $\kappa$ involved when computing $\Omega^2$ changes sign for most pairs $(\bar{\omega}', \bar{\omega})$, which in some cases results in subtracting two quantities close to each other, thereby amplifying errors. This effect can be reduced by increasing the numerical resolution. Actually, a test with $N = 1024$ shows a mean value $\langle \hat{\Omega}^2 \rangle \approx 1.3014 \times 10^{-2}$ and an amplitude around $10^{-3}$ in relative, which is more reasonable.

We list in Table 1 the main global quantities at equilibrium, namely the mass $M$, the moment of inertia $I$, the volume $V$, the angular momentum $J$, the gravitational energy $W$, the kinetic energy $T$, the internal energy $U$, and the virial parameter $VP = W + 2T + U$. We see that $VP/|W|$ is on the order of $10^{-4}$ in absolute, which is very good. We also list the values output by the DROP code, which solves Eq. (1) in full 2D

<footnote n="1">¹ In practice, this level of stability is easily reached. However, it needs to be raised if the numerical resolution is very high.</footnote>

<footnote n="2">² This may appear useless or excessive. In some cases, however, the cycle first converges toward one state and then “bounces” to reach a completely different one. Reaching machine precision offers better confidence in the solution.</footnote>
evolve until stabilization. This can be accomplished following the iterative scheme, for $t \ge 1$:

1. $\bar{\epsilon}^{(t)}$ is obtained from Eq. (7), after integrating Eq. (4b);

2. $\hat{H}^{(t)}$ is computed according to Eq. (4a);

3. $\hat{\rho}^{(t)}$ is obtained from Eq. (4c).

We see that there are two other options (with no significant impact on the performance of the cycle), depending on the order of assignments; in other words, $\bar{\epsilon}$ in Eq. (11b) can be computed using the $\rho$-profile updated from Eq. (11a), or $\rho$ in Eq. (11a) can be computed from the $\bar{\epsilon}$-profile updated from Eq. (11b). At convergence, two successive profiles must be numerically unchanged, in which case the cycle ends. As we traditionally use double precision computers, we take$^{1,2}$ $\delta^{(t)} \le 10^{-14}$, where

$$ \delta^{(t)} = \max(\Delta^{(t)}\hat{\rho}, \Delta^{(t)}\bar{\epsilon}), \quad (12) $$

where we define still for $t \ge 1$ and $\bar{\omega} \in [0, 1]$

$$ \Delta^{(t)}\hat{\rho} = \sup |\hat{\rho}^{(t)}(\bar{\omega}) - \hat{\rho}^{(t-1)}(\bar{\omega})|, \quad (13) $$

and use a similar definition for $\Delta^{(t)}\bar{\epsilon}$.

## 3.2. Boundary conditions

In the standard (single-domain) case, and in the absence of external pressure (see Sect. 5.2), we can take Dirichlet boundary conditions (BCs) at the outer boundary $\bar{\omega} = 1$

$$ \begin{cases} \hat{\rho}(1) = 0 \\ \bar{\epsilon}(1) = \bar{\epsilon}_s \end{cases}, \quad (14) $$

where $\bar{\epsilon}_s = R_p/R_c$ is the axis ratio of the outermost surface. At the center of the body, we have

$$ \begin{cases} \hat{\rho}(0) = 1 \\ \left. \frac{d\bar{\epsilon}^2}{d\bar{\omega}} \right|_{\bar{\omega}=0} = 0 \end{cases}, \quad (15) $$

Unfortunately, $\bar{\epsilon}(0)$ cannot be easily deduced.

## 3.3. Method and implementation (quadrature schemes, $\bar{\omega}$-grid, and seeds)

We see that Eq. (4) involves derivatives and quadratures, and there are different ways to estimate $\rho(\bar{\omega})$ and $\bar{\epsilon}(\bar{\omega})$ for given BCs. Here we decided to recast Eq. (4b) in integral form, i.e.,

$$ \bar{\epsilon}^2(\bar{\omega}) = \bar{\epsilon}_s^2 - \int_{\bar{\omega}}^{1} \left( \frac{d\bar{\epsilon}^2}{d\bar{\omega}} \right) d\bar{\omega}, \quad (16) $$

where the term in parentheses is simply the right-hand side of Eq. (4b) within a minus sign. So, we have only to deal with quadratures. The eccentricity being unknown at the center, we chose to integrate downward, from the surface to the center, with Eq. (14) as Dirichlet BCs. The computational grid is made of $N + 1$ equally spaced nodes: $\bar{\omega}_i \in \{0, \frac{1}{N}, \dots, \bar{\omega}_i = \frac{i}{N}, \dots, 1\}$. We used the trapezoidal rule as the quadrature rule, which is second-order in the step size (more efficient schemes can be used

at this level). According to Eq. (14), still assuming null ambient pressure, we take

$$ \begin{cases} \hat{H}^{(0)}(\bar{\omega}) = \hat{H}_c(1 - \bar{\omega}^2), \\ \bar{\epsilon}^{(0)}(\bar{\omega}) = 1 - (1 - \bar{\epsilon}_s)\bar{\omega}^2 \end{cases}, \quad (17a) $$

$$ \begin{cases} \hat{H}^{(0)}(\bar{\omega}) = \hat{H}_c(1 - \bar{\omega}^2), \\ \bar{\epsilon}^{(0)}(\bar{\omega}) = 1 - (1 - \bar{\epsilon}_s)\bar{\omega}^2 \end{cases}, \quad (17b) $$

as the starting guess. Initially, we thus have $\hat{H}^{(0)}(0) = \hat{H}_c$ at the center, and $\bar{\epsilon}^{(0)}(0) = 1$. These quadratic profiles seem appropriate for a wide range of configurations, although the solutions generally have a nonzero central eccentricity. Obviously, the BCs must be applied at each step in the cycle.

## 3.4. A note about the equation of state

The EOS is the fundamental ingredient. Without loss of generality, we consider a polytropic gas where the pressure is a power law of the mass density, which leads to the expression $H = K(n+1)\rho^{1/n}$, where K is a positive constant and the polytropic index $n > 0$ is finite. In this case the relationship between $\hat{H}$ and $\hat{\rho}$ is

$$ \frac{\hat{H}}{\hat{H}_c} = \hat{\rho}^{1/n}, \quad (18) $$

where $\hat{H}_c = K(n+1)/G\rho_c^{1-1/n}R_c^2$ (this relation assumes $\hat{\rho} \ge 0$). The mass density $\hat{\rho}^{(t)}$ inside the cycle, including the seed, is deduced from the EOS through Eq. (4c).

## 3.5. Example of cycle convergence

The first example is for $\bar{\epsilon}_s = 0.9$ and $n = 3$, and the grid has $N + 1 = 257$ equally spaced nodes. This configuration corresponds, for instance, to a radiation-pressure-dominated ideal gas or to a white dwarf with fully degenerate extremely relativistic electrons. With this parameter, we are already beyond slow rotators. We ran the code based on the SSB approximation. Figure 2 shows the evolution of $\hat{\rho}(\bar{\omega})$ and $\bar{\epsilon}^2(\bar{\omega})$ during the cycle, as well as the deviations $\Delta\hat{\rho}$ and $\Delta\bar{\epsilon}^2$ from one step to the other. Figure 2f gives the convergence parameter $\delta^{(t)}$ defined by Eq. (12) from the beginning to convergence. We see that $\delta^{(t)}$ decreases exponentially with the step $t$ and the algorithm converges quickly on a solution. Convergence is reached after 70 cycles with the current criterion. After step 20, we already have $\delta^{(t)} \sim 10^{-5}$, which is on the order of the accuracy of the numerical scheme (i.e., $1/N^2 \sim 10^{-5}$ with 257 nodes). The next steps are necessary to reach the threshold of $10^{-14}$. We show in Fig. 2c the rotation rate computed from Eq. (10). Unsurprisingly, there is a slight variation with the radius, on the order of $10^{-2}$ in relative. This is due to the approximate nature of the theory of NSFoE. The function $\kappa$ involved when computing $\Omega^2$ changes sign for most pairs $(\bar{\omega}', \bar{\omega})$, which in some cases results in subtracting two quantities close to each other, thereby amplifying errors. This effect can be reduced by increasing the numerical resolution. Actually, a test with $N = 1024$ shows a mean value $\langle \hat{\Omega}^2 \rangle \approx 1.3014 \times 10^{-2}$ and an amplitude around $10^{-3}$ in relative, which is more reasonable.

We list in Table 1 the main global quantities at equilibrium, namely the mass $M$, the moment of inertia $I$, the volume $V$, the angular momentum $J$, the gravitational energy $W$, the kinetic energy $T$, the internal energy $U$, and the virial parameter $VP = W + 2T + U$. We see that $VP/|W|$ is on the order of $10^{-4}$ in absolute, which is very good. We also list the values output by the DROP code, which solves Eq. (1) in full 2D

<footnote n="1">¹ In practice, this level of stability is easily reached. However, it needs to be raised if the numerical resolution is very high.</footnote>

<footnote n="2">² This may appear useless or excessive. In some cases, however, the cycle first converges toward one state and then “bounces” to reach a completely different one. Reaching machine precision offers better confidence in the solution.</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 2.** Radial profiles for $\hat{\rho}$ (panel a) and $\bar{\epsilon}$ (panel b) during the SCF cycle for a few selected steps $t$ (labeled on the curves) of the SCF cycle, and absolute deviations (panels d and e) between steps $t-1$ and $t$. Also shown are the square of the rotation rate (panel c) computed at convergence from Eq. (10) and its mean value (orange), and the convergence parameter $\delta$ computed from Eq. (12) as a function of step $t$ (panel f).</figcaption></figure>

(Huré & Hersant 2017; Basillais & Huré 2021). We see that the agreement between the two methods is quite good, with deviations on the order of a few $10^{-3}$ in relative, while the resolution is moderate.

### 3.6. Example of two-dimensional reconstruction (deprojection)

Once a solution in the form of a pair of profiles $\hat{\rho}(\omega)$ and $\epsilon(\omega)$ is known, we can reconstruct the 2D structure by using Eq. (2). For this second example, we ran the SCF code under the same conditions as above, but for $\bar{\epsilon}_s = 0.75$ and $n = 1.5$, which could correspond to a fast-rotating fully convective star. We note that $\bar{\epsilon}_s \approx 0.74$ for Achernar (Domiciano de Souza et al. 2014). We can compare the $\rho$-map obtained in this way to the field output by the DROP code. This is shown in Fig. 3. We see that the largest differences are mainly located close to the surface, as the true solution is slightly sub-elliptical. The deviations do not exceed about $3 \times 10^{-3}$ in absolute (this value is on the order of the virial

<table><caption>**Table 1.** Global quantities computed for configuration A discussed in Sect. 3.3 with the SSB approximation for $N = 256$ (last column).</caption><thead><tr><th>Configuration A</th><th>DROP<sup>(†)</sup> code</th><th>This work</th></tr></thead><tbody><tr><td>$\bar{\epsilon}_s$</td><td>← 0.900</td><td>← 0.900</td></tr><tr><td>$n$</td><td>← 3.0</td><td>← 3.0</td></tr><tr><td>$M/(\rho_c R_e^3)$</td><td>5.9999 × 10<sup>-2</sup></td><td>6.0062 × 10<sup>-2</sup></td></tr><tr><td>$I/(\rho_c R_e^5)$</td><td>3.9410 × 10<sup>-3</sup></td><td>3.9475 × 10<sup>-3</sup></td></tr><tr><td>$V/R_e^3$</td><td>3.7436</td><td>3.7699</td></tr><tr><td>$\hat{\Omega}^2$</td><td>1.3006 × 10<sup>-2</sup></td><td>1.3024 × 10<sup>-2</sup><br>(*)1.3014 × 10<sup>-2</sup></td></tr><tr><td>$J/(G\rho_c^3 R_e^{10})^{1/2}$</td><td>4.4945 × 10<sup>-4</sup></td><td>4.5051 × 10<sup>-4</sup></td></tr><tr><td>$-W/(G\rho_c^2 R_e^5)$</td><td>5.8487 × 10<sup>-3</sup></td><td>5.8599 × 10<sup>-3</sup></td></tr><tr><td>$T/(G\rho_c^2 R_e^5)$</td><td>2.5629 × 10<sup>-5</sup></td><td>2.5707 × 10<sup>-5</sup></td></tr><tr><td>$U/(G\rho_c^2 R_e^5)$</td><td>5.7976 × 10<sup>-3</sup></td><td>5.8077 × 10<sup>-3</sup></td></tr><tr><td>$[VP/W]$</td><td>4 × 10<sup>-5</sup></td><td>1 × 10<sup>-4</sup></td></tr></tbody></table>

**Notes.** The rotation rate is a mean value (see text). Also listed (second column) are values output from the DROP code at a resolution 257 × 257. $\leftarrow$ input data; $^{(†)}$ see Huré & Hersant (2017) and Basillais & Huré (2021); $^{(*)}$ value obtained for $N = 1024$.

<table><caption>**Table 2.** Same caption and same conditions as for Table 1, but for a fast rotator.</caption><thead><tr><th>Configuration B</th><th>DROP<sup>(†)</sup> code</th><th>This work</th></tr></thead><tbody><tr><td>$\bar{\epsilon}_s$</td><td>← 0.750</td><td>← 0.750</td></tr><tr><td>$n$</td><td>← 1.5</td><td>← 1.5</td></tr><tr><td>$M/(\rho_c R_e^3)$</td><td>4.3026 × 10<sup>-1</sup></td><td>4.3397 × 10<sup>-1</sup></td></tr><tr><td>$I/(\rho_c R_e^5)$</td><td>7.4701 × 10<sup>-2</sup></td><td>7.5961 × 10<sup>-2</sup></td></tr><tr><td>$V/R_e^3$</td><td>3.0297</td><td>3.1415</td></tr><tr><td>$\hat{\Omega}^2$</td><td>2.2662 × 10<sup>-1</sup></td><td>2.2808 × 10<sup>-1</sup><br>(*)2.2807 × 10<sup>-1</sup></td></tr><tr><td>$J/(G\rho_c^3 R_e^{10})^{1/2}$</td><td>3.5561 × 10<sup>-2</sup></td><td>3.6278 × 10<sup>-2</sup></td></tr><tr><td>$-W/(G\rho_c^2 R_e^5)$</td><td>1.8344 × 10<sup>-1</sup></td><td>1.8584 × 10<sup>-1</sup></td></tr><tr><td>$T/(G\rho_c^2 R_e^5)$</td><td>8.4644 × 10<sup>-3</sup></td><td>8.6630 × 10<sup>-3</sup></td></tr><tr><td>$U/(G\rho_c^2 R_e^5)$</td><td>1.6651 × 10<sup>-1</sup></td><td>1.6790 × 10<sup>-1</sup></td></tr><tr><td>$[VP/W]$</td><td>1 × 10<sup>-5</sup></td><td>3 × 10<sup>-3</sup></td></tr></tbody></table>

parameter, see below). This result is remarkable, in particular because the surface is the place where the mass density is small and vanishes. The best agreement is observed at the center, at the pole, and at the equator, with absolute deviations of about $10^{-5}$. The main quantities are reported in Table 2. We see that most global quantities are correctly reproduced within a percent, which is satisfying for a rotator this fast.

### 3.7. Varying the surface axis ratio and the polytropic index: Computing vs. precision

We performed similar comparisons by varying the surface axis ratio $\bar{\epsilon}_s$ and the polytropic index $n$, again with a moderate numerical resolution corresponding to $N = 256$. The results are gathered in Table F.1, where we list the mass, the rotation rate, the relative virial parameter and the root mean square (RMS) difference between the two structures (mass, density). There are three series. For the first the index and the resolution are held fixed ($n = 1.5, N = 256$) and the surface eccentricity increases to the critical rotation, at about $\bar{\epsilon}_s \approx 0.62$. In the second series the surface eccentricity and the resolution are fixed ($\bar{\epsilon}_s = 0.95$,
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Radial profiles for $\hat{\rho}$ (panel a) and $\bar{\varepsilon}$ (panel b) during the SCF cycle for a few selected steps $t$ (labeled on the curves) of the SCF cycle, and absolute deviations (panels d and e) between steps $t-1$ and $t$. Also shown are the square of the rotation rate (panel c) computed at convergence from Eq. (10) and its mean value (orange), and the convergence parameter $\delta$ computed from Eq. (12) as a function of step $t$ (panel f).</figcaption></figure>

(Huré & Hersant 2017; Basillais & Huré 2021). We see that the agreement between the two methods is quite good, with deviations on the order of a few $10^{-3}$ in relative, while the resolution is moderate.

### 3.6. Example of two-dimensional reconstruction (deprojection)

Once a solution in the form of a pair of profiles $\hat{\rho}(\omega)$ and $\varepsilon(\omega)$ is known, we can reconstruct the 2D structure by using Eq. (2). For this second example, we ran the SCF code under the same conditions as above, but for $\bar{\varepsilon}_s = 0.75$ and $n = 1.5$, which could correspond to a fast-rotating fully convective star. We note that $\bar{\varepsilon}_s \approx 0.74$ for Achernar (Domiciano de Souza et al. 2014). We can compare the $\rho$-map obtained in this way to the field output by the DROP code. This is shown in Fig. 3. We see that the largest differences are mainly located close to the surface, as the true solution is slightly sub-elliptical. The deviations do not exceed about $3 \times 10^{-3}$ in absolute (this value is on the order of the virial

<table><caption>**Table 1.** Global quantities computed for configuration A discussed in Sect. 3.3 with the SSB approximation for $N = 256$ (last column).</caption><thead><tr><th>Configuration A</th><th>DROP<sup>(†)</sup> code</th><th>This work</th></tr></thead><tbody><tr><td>$\bar{\varepsilon}_s$</td><td>$\leftarrow 0.900$</td><td>$\leftarrow 0.900$</td></tr><tr><td>$n$</td><td>$\leftarrow 3.0$</td><td>$\leftarrow 3.0$</td></tr><tr><td>$M/(\rho_c R_e^3)$</td><td>$5.9999 \times 10^{-2}$</td><td>$6.0062 \times 10^{-2}$</td></tr><tr><td>$I/(\rho_c R_e^5)$</td><td>$3.9410 \times 10^{-3}$</td><td>$3.9475 \times 10^{-3}$</td></tr><tr><td>$V/R_e^3$</td><td>$3.7436$</td><td>$3.7699$</td></tr><tr><td>$\hat{\Omega}^2$</td><td>$1.3006 \times 10^{-2}$</td><td>$1.3024 \times 10^{-2}$<br>($\ast$)$1.3014 \times 10^{-2}$</td></tr><tr><td>$J/(G\rho_c^3 R_e^{10})^{1/2}$</td><td>$4.4945 \times 10^{-4}$</td><td>$4.5051 \times 10^{-4}$</td></tr><tr><td>$-W/(G\rho_c^2 R_e^5)$</td><td>$5.8487 \times 10^{-3}$</td><td>$5.8599 \times 10^{-3}$</td></tr><tr><td>$T/(G\rho_c^2 R_e^5)$</td><td>$2.5629 \times 10^{-5}$</td><td>$2.5707 \times 10^{-5}$</td></tr><tr><td>$U/(G\rho_c^2 R_e^5)$</td><td>$5.7976 \times 10^{-3}$</td><td>$5.8077 \times 10^{-3}$</td></tr><tr><td>$[VP/W]$</td><td>$4 \times 10^{-5}$</td><td>$1 \times 10^{-4}$</td></tr></tbody></table>

**Notes.** The rotation rate is a mean value (see text). Also listed (second column) are values output from the DROP code at a resolution $257 \times 257$. $\leftarrow$ input data; $^{(†)}$ see Huré & Hersant (2017) and Basillais & Huré (2021); $^{(\ast)}$ value obtained for $N = 1024$.

<table><caption>**Table 2.** Same caption and same conditions as for Table 1, but for a fast rotator.</caption><thead><tr><th>Configuration B</th><th>DROP<sup>(†)</sup> code</th><th>This work</th></tr></thead><tbody><tr><td>$\bar{\varepsilon}_s$</td><td>$\leftarrow 0.750$</td><td>$\leftarrow 0.750$</td></tr><tr><td>$n$</td><td>$\leftarrow 1.5$</td><td>$\leftarrow 1.5$</td></tr><tr><td>$M/(\rho_c R_e^3)$</td><td>$4.3026 \times 10^{-1}$</td><td>$4.3397 \times 10^{-1}$</td></tr><tr><td>$I/(\rho_c R_e^5)$</td><td>$7.4701 \times 10^{-2}$</td><td>$7.5961 \times 10^{-2}$</td></tr><tr><td>$V/R_e^3$</td><td>$3.0297$</td><td>$3.1415$</td></tr><tr><td>$\hat{\Omega}^2$</td><td>$2.2662 \times 10^{-1}$</td><td>$2.2808 \times 10^{-1}$<br>($\ast$)$2.2807 \times 10^{-1}$</td></tr><tr><td>$J/(G\rho_c^3 R_e^{10})^{1/2}$</td><td>$3.5561 \times 10^{-2}$</td><td>$3.6278 \times 10^{-2}$</td></tr><tr><td>$-W/(G\rho_c^2 R_e^5)$</td><td>$1.8344 \times 10^{-1}$</td><td>$1.8584 \times 10^{-1}$</td></tr><tr><td>$T/(G\rho_c^2 R_e^5)$</td><td>$8.4644 \times 10^{-3}$</td><td>$8.6630 \times 10^{-3}$</td></tr><tr><td>$U/(G\rho_c^2 R_e^5)$</td><td>$1.6651 \times 10^{-1}$</td><td>$1.6790 \times 10^{-1}$</td></tr><tr><td>$[VP/W]$</td><td>$1 \times 10^{-5}$</td><td>$3 \times 10^{-3}$</td></tr></tbody></table>

parameter, see below). This result is remarkable, in particular because the surface is the place where the mass density is small and vanishes. The best agreement is observed at the center, at the pole, and at the equator, with absolute deviations of about $10^{-5}$. The main quantities are reported in Table 2. We see that most global quantities are correctly reproduced within a percent, which is satisfying for a rotator this fast.

### 3.7. Varying the surface axis ratio and the polytropic index: Computing vs. precision

We performed similar comparisons by varying the surface axis ratio $\bar{\varepsilon}_s$ and the polytropic index $n$, again with a moderate numerical resolution corresponding to $N = 256$. The results are gathered in Table F.1, where we list the mass, the rotation rate, the relative virial parameter and the root mean square (RMS) difference between the two structures (mass, density). There are three series. For the first the index and the resolution are held fixed ($n = 1.5, N = 256$) and the surface eccentricity increases to the critical rotation, at about $\bar{\varepsilon}_s \approx 0.62$. In the second series the surface eccentricity and the resolution are fixed ($\bar{\varepsilon}_s = 0.95$,
<figure><img src="image_3.png" /><figcaption>**Fig. 3.** Color-coded mass density map (log. scale) computed from Eq. (4) that takes advantage of dimension reduction (top) and absolute difference (log. scale) with the reference DROP code (bottom), for configuration B. The parameters are $\bar{\epsilon}_s = 0.75$ and $n = 1.5$, with 257 nodes per direction (see also Table 2).</figcaption></figure>

$N = 256$) and $n$ varies from 0.5 to 4. In the last series, the configuration is fixed ($\bar{\epsilon}_s = 0.95$, $n = 1.5$) and we increase the resolution from 32 to 2048. We see that the agreement is good: the maximum RMS value is $4 \times 10^{-3}$ and the mass and rotation rate are well reproduced within a percent. From the virial parameters and the RMS, two trends are seen. First, the method is less and less accurate as the axis ratio decreases, especially for a “hard” EOS. This is expected because the true surface deviates more from a perfect spheroid as the rotation becomes faster. Second, the method is increasingly accurate as the polytropic index increases (“soft” EOS), even close to critical rotations. As $n$ increases, the density becomes peaked at the center, and the contribution to the mass (and potential) of the “wings” becomes small to negligible. This is visible in Fig. 2a, for instance. For $n = 3$ the mass density vanishes quickly toward the surface ($\hat{\rho} \lesssim 10^{-3}$ at $\bar{\omega} \gtrsim 0.7$).

The computing time reported in the table is obtained on a standard laptop, without any specification optimization. For $n = 1.5$ and $\bar{\epsilon}_s = 0.95$, the number of iterations is about 30 and is not sensitive to $N$. We find $\sim (N/581.295)^{1.930}$ s for convergence with the SSB approximation, to be compared to the full 2D problem $\sim (N/77.824)^{3.089}$ s.

## 4. More tests

In order to better see the power of the method and its flexibility, we present in this section several tests, including static and rotating configurations (see next section for systems with mass

<figure><img src="image_4.png" /><figcaption>**Fig. 4.** Mass density profile obtained by solving Eq. (4) with the SCF algorithm ($N = 256$) in the static case ($\bar{\epsilon}_s = 1$) and for $n = 1$, compared to the corresponding exact solution of the Lane-Emden equation (left panel). Also shown is the absolute deviation between these two solutions, for several values of $N$ (right panel); see also Table 3.</figcaption></figure>

<table><caption>**Table 3.** Mass and moment of inertia for the static polytropic system with $n = 1$ (second column), and values computed from the SCF method, using the 2D code (third column) and after dimension reduction (fourth column).</caption><thead><tr><th></th><th>Exact</th><th>DROP code</th><th>This work</th></tr></thead><tbody><tr><td>$M/(\rho_c R_c^3)$</td><td>1.27324</td><td>1.27322</td><td>1.27323</td></tr><tr><td>$I/(\rho_c R_c^5)$</td><td>0.33280</td><td>0.33279</td><td>0.33279</td></tr><tr><td>Steps for convergence</td><td></td><td>29</td><td>20</td></tr><tr><td>Virial parameter</td><td>0</td><td>$4 \times 10^{-6}$</td><td>$2 \times 10^{-5}$</td></tr></tbody></table>

density jumps). The convergence criterion and the numerical resolution are (unless stated otherwise) the same as in Sect. 3.5.

### 4.1. Static case with index unity

We consider here a static gas with polytropic index $n = 1$ (see Sect. 1). All isobars are spherical and $\epsilon = 0$ for any $\bar{\omega}$. In this case the solution of the Lane-Emden equation is analytical$^3$, namely

$$ \hat{\rho}(\bar{\omega}) = \frac{\sin(\pi\bar{\omega})}{\pi\bar{\omega}}, \qquad (19) $$

which enables a direct comparison. The theory of NSFoE is exact for spherical configuration because the confocal parameters are null, which is therefore the case of Eq. (4). If we inject $\epsilon(\bar{\omega}) = 0$ in the four functions $\kappa$, $\chi$, $\mu$, and $\eta$, an expansion is required in the limit $\epsilon \to 0$. We find that $\epsilon = 0$ is in fact a regular singular point, and it follows that $\kappa = \chi = 0$ (see Appendix D). Thus, we recover that the body is nonrotating and that all isopycnics are spherical. We have applied the SSB approximation. The results are shown in Fig. 4. The main output data are listed in Table 3. We see that the deviation from the analytical solution is of order $10^{-5}$, and it decreases as $N$ increases. This occurs because deviations are directly linked to the order of the quadra-ture scheme in this case.

<footnote>$^3$ This is also the case for $n = 0$ and $n = 5$, but $n = 0$ represents a homogeneous object and this test was already considered in Paper IV. The case with $n = 5$ has an infinite radial extent and mass; see Sect. 5.2.</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Color-coded mass density map (log. scale) computed from Eq. (4) that takes advantage of dimension reduction (top) and absolute difference (log. scale) with the reference DROP code (bottom), for configuration B. The parameters are $\bar{\epsilon}_s = 0.75$ and $n = 1.5$, with 257 nodes per direction (see also Table 2).</figcaption></figure>

$N = 256$) and $n$ varies from 0.5 to 4. In the last series, the configuration is fixed ($\bar{\epsilon}_s = 0.95$, $n = 1.5$) and we increase the resolution from 32 to 2048. We see that the agreement is good: the maximum RMS value is $4 \times 10^{-3}$ and the mass and rotation rate are well reproduced within a percent. From the virial parameters and the RMS, two trends are seen. First, the method is less and less accurate as the axis ratio decreases, especially for a “hard” EOS. This is expected because the true surface deviates more from a perfect spheroid as the rotation becomes faster. Second, the method is increasingly accurate as the polytropic index increases (“soft” EOS), even close to critical rotations. As $n$ increases, the density becomes peaked at the center, and the contribution to the mass (and potential) of the “wings” becomes small to negligible. This is visible in Fig. 2a, for instance. For $n = 3$ the mass density vanishes quickly toward the surface ($\hat{\rho} \lesssim 10^{-3}$ at $\bar{\omega} \gtrsim 0.7$).

The computing time reported in the table is obtained on a standard laptop, without any specification optimization. For $n = 1.5$ and $\bar{\epsilon}_s = 0.95$, the number of iterations is about 30 and is not sensitive to $N$. We find $\sim (N/581.295)^{1.930}$ s for convergence with the SSB approximation, to be compared to the full 2D problem $\sim (N/77.824)^{3.089}$ s.

## 4. More tests

In order to better see the power of the method and its flexibility, we present in this section several tests, including static and rotating configurations (see next section for systems with mass

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Mass density profile obtained by solving Eq. (4) with the SCF algorithm ($N = 256$) in the static case ($\bar{\epsilon}_s = 1$) and for $n = 1$, compared to the corresponding exact solution of the Lane-Emden equation (left panel). Also shown is the absolute deviation between these two solutions, for several values of $N$ (right panel); see also Table 3.</figcaption></figure>

<table><caption>**Table 3.** Mass and moment of inertia for the static polytropic system with $n = 1$ (second column), and values computed from the SCF method, using the 2D code (third column) and after dimension reduction (fourth column).</caption><thead><tr><th></th><th>Exact</th><th>DROP code</th><th>This work</th></tr></thead><tbody><tr><td>$M/(\rho_c R_c^3)$</td><td>1.27324</td><td>1.27322</td><td>1.27323</td></tr><tr><td>$I/(\rho_c R_c^5)$</td><td>0.33280</td><td>0.33279</td><td>0.33279</td></tr><tr><td>Steps for convergence</td><td></td><td>29</td><td>20</td></tr><tr><td>Virial parameter</td><td>0</td><td>$4 \times 10^{-6}$</td><td>$2 \times 10^{-5}$</td></tr></tbody></table>

density jumps). The convergence criterion and the numerical resolution are (unless stated otherwise) the same as in Sect. 3.5.

### 4.1. Static case with index unity

We consider here a static gas with polytropic index $n = 1$ (see Sect. 1). All isobars are spherical and $\epsilon = 0$ for any $\bar{\omega}$. In this case the solution of the Lane-Emden equation is analytical$^3$, namely

$$ \hat{\rho}(\bar{\omega}) = \frac{\sin(\pi\bar{\omega})}{\pi\bar{\omega}}, \qquad (19) $$

which enables a direct comparison. The theory of NSFoE is exact for spherical configuration because the confocal parameters are null, which is therefore the case of Eq. (4). If we inject $\epsilon(\bar{\omega}) = 0$ in the four functions $\kappa, \chi, \mu$, and $\eta$, an expansion is required in the limit $\epsilon \to 0$. We find that $\epsilon = 0$ is in fact a regular singular point, and it follows that $\kappa = \chi = 0$ (see Appendix D). Thus, we recover that the body is nonrotating and that all isopycnics are spherical. We have applied the SSB approximation. The results are shown in Fig. 4. The main output data are listed in Table 3. We see that the deviation from the analytical solution is of order $10^{-5}$, and it decreases as $N$ increases. This occurs because deviations are directly linked to the order of the quadra-ture scheme in this case.

<footnote>$^3$ This is also the case for $n = 0$ and $n = 5$, but $n = 0$ represents a homogeneous object and this test was already considered in Paper IV. The case with $n = 5$ has an infinite radial extent and mass; see Sect. 5.2.</footnote>
<figure><img src="image_5.png" /><figcaption>Fig. 5. Mass density profile and squared eccentricity profile (left panels) obtained by solving Eq. (4) with the SCF algorithm, and comparison with the exact nonrotating first-order solution from Clairaut’s fundamental equation (right panels).</figcaption></figure>

## 4.2. The case of slow rotation

In the limit of slow rotation, various approximations can be found (see, e.g., Chandrasekhar 1933). Of particular interest is Clairaut’s theory, which is first-order accurate in the square of the eccentricity, $\epsilon^2 \ll 1$. It happens that there are a few closed-form solutions to Clairaut’s second-order differential equation compatible with physically realistic BCs. Among them, Tisserand (1891) and Marchenko (2000) have a Legendre-Laplace solution for $n = 1$: the solution is the same as for the nonrotating case (i.e., a sine cardinal for the mass density)$^4$. For Eq. (19), the eccentricity profile reads

$$ \bar{\epsilon}^2(\omega) = 1 - (1 - \bar{\epsilon}_s^2) \frac{(\pi^2 \omega^2 - 3) \sin(\pi\omega) + 3\pi\omega \cos(\pi\omega)}{3\omega^2 [\pi\omega \cos(\pi\omega) - \sin(\pi\omega)]} \quad (20) $$

We ran the code for $n = 1$ and $\bar{\epsilon}_s = 0.99$. The results are displayed in Fig. 5. We see that the mass density profile of the rotating case departs only slightly from a sine cardinal (the deviation is of order $10^{-3}$). More importantly, the actual method yields an eccentricity profile that is very close to Eq. (20), within a few $10^{-5}$ in absolute, which is highly satisfying. This confirms the efficiency of the SSB approximation at slow rotation. Figure 6 displays the mass density structure and the deviation from the reference structure, with a RMS value of $2 \times 10^{-6}$.

## 4.3. A case of moderate rotation

An interesting test concerns the transition from slow to moderate rotators. We performed a new run for $n = 1$ and $\bar{\epsilon}_s = 0.95$ (which is close to the axis ratio of Jupiter, for instance). The results are shown in Fig. 7 after reconstruction of the 2D map for the mass

<figure><img src="image_6.png" /><figcaption>Fig. 6. Same as Fig. 3, but for $\bar{\epsilon}_s = 0.99$ and $n = 1$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 7. Same as Fig. 3, but for $\bar{\epsilon}_s = 0.95$ and $n = 1$.</figcaption></figure>

<footnote>$^4$As Clairaut’s equation is first order accurate in $\epsilon^2$, only the zeroth-order on $\hat{\rho}$ is needed to have $\epsilon^2$ in the limit of slow rotations.</footnote>
<figure><img src="image_6.png" /><figcaption>Fig. 5. Mass density profile and squared eccentricity profile (left panels) obtained by solving Eq. (4) with the SCF algorithm, and comparison with the exact nonrotating first-order solution from Clairaut’s fundamental equation (right panels).</figcaption></figure>

## 4.2. The case of slow rotation

In the limit of slow rotation, various approximations can be found (see, e.g., Chandrasekhar 1933). Of particular interest is Clairaut’s theory, which is first-order accurate in the square of the eccentricity, $\epsilon^2 \ll 1$. It happens that there are a few closed-form solutions to Clairaut’s second-order differential equation compatible with physically realistic BCs. Among them, Tisserand (1891) and Marchenko (2000) have a Legendre-Laplace solution for $n = 1$: the solution is the same as for the nonrotating case (i.e., a sine cardinal for the mass density)$^4$. For Eq. (19), the eccentricity profile reads

$$ \bar{\epsilon}^2(\omega) = 1 - (1 - \bar{\epsilon}_s^2) \frac{(\pi^2 \omega^2 - 3) \sin(\pi\omega) + 3\pi\omega \cos(\pi\omega)}{3\omega^2 [\pi\omega \cos(\pi\omega) - \sin(\pi\omega)]} \quad (20) $$

We ran the code for $n = 1$ and $\bar{\epsilon}_s = 0.99$. The results are displayed in Fig. 5. We see that the mass density profile of the rotating case departs only slightly from a sine cardinal (the deviation is of order $10^{-3}$). More importantly, the actual method yields an eccentricity profile that is very close to Eq. (20), within a few $10^{-5}$ in absolute, which is highly satisfying. This confirms the efficiency of the SSB approximation at slow rotation. Figure 6 displays the mass density structure and the deviation from the reference structure, with a RMS value of $2 \times 10^{-6}$.

## 4.3. A case of moderate rotation

An interesting test concerns the transition from slow to moderate rotators. We performed a new run for $n = 1$ and $\bar{\epsilon}_s = 0.95$ (which is close to the axis ratio of Jupiter, for instance). The results are shown in Fig. 7 after reconstruction of the 2D map for the mass

<figure><img src="image_7.png" /><figcaption>Fig. 6. Same as Fig. 3, but for $\bar{\epsilon}_s = 0.99$ and $n = 1$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Same as Fig. 3, but for $\bar{\epsilon}_s = 0.95$ and $n = 1$.</figcaption></figure>

<footnote>$^4$ As Clairaut’s equation is first order accurate in $\epsilon^2$, only the zeroth-order on $\hat{\rho}$ is needed to have $\epsilon^2$ in the limit of slow rotations.</footnote>
<figure><img src="image_8.png" /><figcaption>**Fig. 8.** The $(j^2, \omega^2)$ diagram for several polytropic indices $n$ labeled on the curves. The solid lines correspond to the solutions obtained from DROP code. The empty circles are the solutions reported in the tables of Hachisu (1986a). The full circles are the solutions obtained with the SCF method with dimension reduction. The color-coding indicates the decimal logarithm of the relative virial parameter $|VP/W|$.</figcaption></figure>

density. We see that the deviation is again maximum at the surface, midway between the pole and the equator. The highest value of the RMS for the mass density is about $3 \times 10^{-4}$, while the mean value is about $7 \times 10^{-5}$. At the center, near the pole and the equator, the precision is excellent, with more than six correct digits.

### 4.4. Reproduction of Hachisu's $(j^2, \omega^2)$ sequences of equilibrium

In the 1980s Hachisu and collaborators perfomed a wide exploration of figures of equilibrium through a series of fundamental papers (Hachisu 1986a,b). They computed sequences of axially symmetrical equilibria by varying the surface axis ratio $\tilde{\epsilon}_s$, both for stars and rings (out of range here). The configurations are gathered in the $(j^2, \omega^2)$-plane, where $\omega$ and $j$ are the normalized rotation rate $\Omega$ and angular momentum $J$, respectively, defined as

$$ \begin{cases} j^2 = \frac{J^2}{4\pi G M^3 V^{1/3}} \\ \omega^2 = \frac{\Omega^2}{4\pi G M V^{-1}} \end{cases}. \qquad (21) $$

Hachisu (1986a) showed that, in contrast to the Maclaurin uncompressible sequence, all sequences are open, depending on the polytropic index $n > 0$ of the gas: the larger the value of $n$, the larger the gap. Figure 8 shows a few points of the sequences obtained with the SSB approximation for $n = \{0.1, 0.5, 1.5, 3, 4\}$, together with the data extracted from Hachisu's paper (same axis ratio). Values obtained in full 2D with the DROP code are also plotted for comparison. We see that the agreement is excellent, except for points close to critical rotations. This is obviously due to the deviation of the external surface to a spheroid: the volume is overestimated in this case, which causes a double shift in the diagram, toward higher $\omega$-values and lower $j$-values. The best results are obtained for weakly compressible gas. The virial parameters, on the other hand, are better and better as $n$ increases (see Sect. 3.6).

## 5. Cases with mass density jumps

### 5.1. Changes

Bodies made of different inhomogeneous domains separated by mass density jumps are of immense interest. Such cases are studied in Paper IV. We now consider $\mathcal{K}$ domains where the mass

density $\rho_k(\omega)$ varies (domain no. 1 is for the innermost domain and $k = \mathcal{K}$ is for the outermost one). Then we can write the full mass density profile from the center to the surface as

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

where $\omega_k$ is the position, along the equatorial axis, of the mass density jump between domains $k$ and $k + 1$, and $\mathcal{H}$ is Heaviside's distribution. For convenience we set $\rho_{\mathcal{K}+1} = 0$ to keep a single sum in Eq. (22). The radial derivative is then given by the expression

$$ \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), \qquad (23) $$

where $\delta$ is Dirac's distribution. So, for any continuous function $f(\omega', \omega)$ in the interval $[0, 1]^2$, we have

$$ \int_{\hat{\rho}(0)}^{\hat{\rho}(1)} d\hat{\rho}(\omega') f(\omega', \omega) = \sum_{k=1}^{\mathcal{K}} \int_{\hat{\rho}_k(\omega_{k-1})}^{\hat{\rho}_k(\omega_k)} d\hat{\rho}_k(\omega') f(\omega', \omega) - \sum_{k=1}^{\mathcal{K}} \frac{\alpha_k - 1}{\alpha_k} \hat{\rho}_k(\omega_k) f(\omega_k, \omega), \qquad (24) $$

where $\omega_0 = 0$ again for convenience, and the mass density jump $\alpha_k$ is defined by

$$ \alpha_k = \frac{\rho_k(\omega_k)}{\rho_{k+1}(\omega_k)}. \qquad (25) $$

It follows that the structure of a rigidly rotating body made of several inhomogeneous domains can be treated with dimension reduction by solving Eqs. (4a) and (4b), where all the integrals are estimated from Eq. (24).

### 5.2. Presence of an ambient pressure

As the present formalism relies mainly on a barotropic EOS, we assume the external pressure to be constant along the outermost surface of the object: $P_{\text{amb}}$. Thus, this value corresponds to a
<figure><img src="image_9.png" /><figcaption>Fig. 8. The ($j^2$, $\omega^2$) diagram for several polytropic indices $n$ labeled on the curves. The solid lines correspond to the solutions obtained from DROP code. The empty circles are the solutions reported in the tables of Hachisu (1986a). The full circles are the solutions obtained with the SCF method with dimension reduction. The color-coding indicates the decimal logarithm of the relative virial parameter $|VP/W|$.</figcaption></figure>

density. We see that the deviation is again maximum at the surface, midway between the pole and the equator. The highest value of the RMS for the mass density is about $3 \times 10^{-4}$, while the mean value is about $7 \times 10^{-5}$. At the center, near the pole and the equator, the precision is excellent, with more than six correct digits.

### 4.4. Reproduction of Hachisu's ($j^2$, $\omega^2$) sequences of equilibrium

In the 1980s Hachisu and collaborators perfomed a wide exploration of figures of equilibrium through a series of fundamental papers (Hachisu 1986a,b). They computed sequences of axially symmetrical equilibria by varying the surface axis ratio $\bar{\epsilon}_s$, both for stars and rings (out of range here). The configurations are gathered in the ($j^2$, $\omega^2$)-plane, where $\omega$ and $j$ are the normalized rotation rate $\Omega$ and angular momentum $J$, respectively, defined as

$$ \begin{cases} j^2 = \frac{J^2}{4\pi G M^3 V^{1/3}} \\ \omega^2 = \frac{\Omega^2}{4\pi G M V^{-1}} \end{cases}. \qquad (21) $$

Hachisu (1986a) showed that, in contrast to the Maclaurin uncompressible sequence, all sequences are open, depending on the polytropic index $n > 0$ of the gas: the larger the value of $n$, the larger the gap. Figure 8 shows a few points of the sequences obtained with the SSB approximation for $n = \{0.1, 0.5, 1.5, 3, 4\}$, together with the data extracted from Hachisu's paper (same axis ratio). Values obtained in full 2D with the DROP code are also plotted for comparison. We see that the agreement is excellent, except for points close to critical rotations. This is obviously due to the deviation of the external surface to a spheroid: the volume is overestimated in this case, which causes a double shift in the diagram, toward higher $\omega$-values and lower $j$-values. The best results are obtained for weakly compressible gas. The virial parameters, on the other hand, are better and better as $n$ increases (see Sect. 3.6).

## 5. Cases with mass density jumps

### 5.1. Changes

Bodies made of different inhomogeneous domains separated by mass density jumps are of immense interest. Such cases are studied in Paper IV. We now consider $\mathcal{K}$ domains where the mass

density $\rho_k(\bar{\omega})$ varies (domain no. 1 is for the innermost domain and $k = \mathcal{K}$ is for the outermost one). Then we can write the full mass density profile from the center to the surface as

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

where $\bar{\omega}_k$ is the position, along the equatorial axis, of the mass density jump between domains $k$ and $k + 1$, and $\mathcal{H}$ is Heaviside's distribution. For convenience we set $\rho_{\mathcal{K}+1} = 0$ to keep a single sum in Eq. (22). The radial derivative is then given by the expression

$$ \frac{d\hat{\rho}}{d\bar{\omega}} = \sum_{k=1}^{\mathcal{K}} \left[ \frac{d\hat{\rho}_k}{d\bar{\omega}} - \frac{d\hat{\rho}_{k+1}}{d\bar{\omega}} \right] \mathcal{H}(\bar{\omega}_k - \bar{\omega}) - \sum_{k=1}^{\mathcal{K}} [\hat{\rho}_k(\bar{\omega}) - \hat{\rho}_{k+1}(\bar{\omega})] \delta(\bar{\omega}_k - \bar{\omega}), \qquad (23) $$

where $\delta$ is Dirac's distribution. So, for any continuous function $f(\bar{\omega}', \bar{\omega})$ in the interval $[0, 1]^2$, we have

$$ \int_{\hat{\rho}(0)}^{\hat{\rho}(1)} d\hat{\rho}(\bar{\omega}') f(\bar{\omega}', \bar{\omega}) = \sum_{k=1}^{\mathcal{K}} \int_{\hat{\rho}_k(\bar{\omega}_{k-1})}^{\hat{\rho}_k(\bar{\omega}_k)} d\hat{\rho}_k(\bar{\omega}') f(\bar{\omega}', \bar{\omega}) - \sum_{k=1}^{\mathcal{K}} \frac{\alpha_k - 1}{\alpha_k} \hat{\rho}_k(\bar{\omega}_k) f(\bar{\omega}_k, \bar{\omega}), \qquad (24) $$

where $\bar{\omega}_0 = 0$ again for convenience, and the mass density jump $\alpha_k$ is defined by

$$ \alpha_k = \frac{\rho_k(\bar{\omega}_k)}{\rho_{k+1}(\bar{\omega}_k)}. \qquad (25) $$

It follows that the structure of a rigidly rotating body made of several inhomogeneous domains can be treated with dimension reduction by solving Eqs. (4a) and (4b), where all the integrals are estimated from Eq. (24).

### 5.2. Presence of an ambient pressure

As the present formalism relies mainly on a barotropic EOS, we assume the external pressure to be constant along the outermost surface of the object: $P_{\text{amb}}$. Thus, this value corresponds to a
<figure><img src="image_9.png" /><figcaption>**Fig. 9.** Mass density profile obtained by solving Eq. (4b) with the SCF algorithm and $N = 256$, with configuration $\bar{\epsilon}_s = 1$, $n = 5$, and $\hat{\rho}_{amb} = 0.04$ (configuration C), compared to the exact solution of Schuster (left panel), and the absolute deviation between these two solutions for several values of $N$ (right panel).</figcaption></figure>

<table><caption>**Table 4.** Data associated with Fig. 9 (static pressurized case).</caption><thead><tr><th></th><th>Exact</th><th>This work</th></tr></thead><tbody><tr><td>$M/\rho_c R_e^3$</td><td>0.60718</td><td>0.60720</td></tr><tr><td>$I/\rho_c R_e^5$</td><td>0.15262</td><td>0.15262</td></tr><tr><td>Steps for convergence</td><td>0</td><td>21</td></tr><tr><td>Virial parameter</td><td>0</td><td>$2 \times 10^{-5}$</td></tr></tbody></table>

cutoff for the mass density at $\hat{\rho}_{amb}$, and then to a mass density jump at $\bar{\omega} = 1$. In the case of a single domain object ($K = 1$ and $\bar{\omega}_1 = 1$ here), we have

$$ \hat{\rho}(\bar{\omega}) = \hat{\rho}_1(\bar{\omega})H(1 - \bar{\omega}), \quad (26) $$

where $\hat{\rho}_1(1) = \hat{\rho}_{amb}$. An interesting test case for ambient pressure is the $n = 5$ polytrope. The analytical solution in the nonrotating case is due to Schuster (e.g., Horedt 2004; see note 3). In the present context, the solution must be truncated at the right isobar $P(1) = P_{amb}$. The analytical solution becomes

$$ \hat{\rho}(\bar{\omega}) = \frac{1}{[1 + \bar{\omega}^2(\hat{\rho}_{amb}^{-2/5} - 1)]^{5/2}}, \quad (27) $$

which verifies $\hat{\rho}(0) = 1$ and $\hat{\rho}(1) = \hat{\rho}_{amb}$. We computed the solution from the SSB approximation. We note that the method is not appropriate for the case $\hat{\rho}_{amb} = 0$ because the mass and especially the radius are infinite for $P_{amb} = 0$. The reconstructed mass density map is displayed in Fig. 9, and the main data are listed in Table 4. As for the $n = 1$ static case, we see that the agreement is excellent and the deviation from the exact solution depends just on the resolution, as expected. We note that the virial parameter was adapted to the context of an ambient pressure. We have

$$ VP = W + 2T + U - U_{amb}, \quad (28) $$

where $U_{amb} = 3P_{amb}V$ (see, e.g., Cox & Giuli 1968).

In the rotating case, Schuster's solution is not adapted anymore, and we switch back to the DROP code for the reference. For this illustration we consider the configuration given in Table 5. It is an intermediate rotator with the same index and ambient pressure as before. A comparison between the full 2D solution and the reconstructed density field, as well as the difference between

<figure><img src="image_10.png" /><figcaption>**Fig. 10.** Same as Fig. 3, but for $\bar{\epsilon}_s = 0.90$, $n = 5.0$, and $\hat{\rho}_{amb} = 0.04$ (configuration C); see also Table 5.</figcaption></figure>

<table><caption>**Table 5.** Same caption and same conditions as for Table 2, but for the rotating case with ambient pressure.</caption><thead><tr><th>Configuration C</th><th>DROP<sup>(†)</sup> code</th><th>This work</th></tr></thead><tbody><tr><td>$\bar{\epsilon}_s$</td><td>$\leftarrow 0.900$</td><td>$\leftarrow 0.900$</td></tr><tr><td>$n$</td><td>$\leftarrow 5$</td><td>$\leftarrow 5$</td></tr><tr><td>$\hat{\rho}_{amb}$</td><td>$\leftarrow 0.04$</td><td>$\leftarrow 0.04$</td></tr><tr><td>$M/(\rho_c R_e^3)$</td><td>$5.2781 \times 10^{-1}$</td><td>$5.2871 \times 10^{-1}$</td></tr><tr><td>$I/(\rho_c R_e^5)$</td><td>$1.3129 \times 10^{-1}$</td><td>$1.3169 \times 10^{-2}$</td></tr><tr><td>$V/R_e^3$</td><td>$3.7596$</td><td>$3.7699$</td></tr><tr><td>$\hat{\Omega}^2$</td><td>$7.7911 \times 10^{-2}$</td><td>$7.8114 \times 10^{-2}$<br>($\star$)$7.8109 \times 10^{-2}$</td></tr><tr><td>$J/(G\rho_c^3 R_e^{10})^{1/2}$</td><td>$3.6648 \times 10^{-2}$</td><td>$3.6805 \times 10^{-2}$</td></tr><tr><td>$-W/(G\rho_c^2 R_e^5)$</td><td>$2.2649 \times 10^{-1}$</td><td>$2.2709 \times 10^{-1}$</td></tr><tr><td>$T/(G\rho_c^2 R_e^5)$</td><td>$5.1147 \times 10^{-3}$</td><td>$5.1434 \times 10^{-3}$</td></tr><tr><td>$U/(G\rho_c^5 R_e^5)$</td><td>$2.7139 \times 10^{-1}$</td><td>$2.7190 \times 10^{-1}$</td></tr><tr><td>$U_{amb}/(G\rho_c^2 R_e^5)$</td><td>$5.5118 \times 10^{-2}$</td><td>$5.5283 \times 10^{-2}$</td></tr><tr><td>$[VP/W]$</td><td>$5 \times 10^{-5}$</td><td>$8 \times 10^{-4}$</td></tr></tbody></table>

the two maps, is shown in Fig. 10. We see that the RMS error<fnref n="5" /> is on the order of $10^{-4}$. The main global quantities are given in Table 5. The agreement is once again very good, with local deviations on the order of $10^{-4}$ on the density profile. For global properties, the relative deviations are on the order of $10^{-3}$.

<footnote n="5">⁵ As the true boundary of the fluid is slightly sub-elliptical, some nodes are out of the fluid in the 2D problem, while $\hat{\rho} = \hat{\rho}_{amb}$. To avoid an excessive overestimation of the RMS error, we subtracted the ambient mass density from both maps to obtain this value.</footnote>
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Mass density profile obtained by solving Eq. (4b) with the SCF algorithm and $N = 256$, with configuration $\bar{\epsilon}_s = 1$, $n = 5$, and $\hat{\rho}_{amb} = 0.04$ (configuration C), compared to the exact solution of Schuster (left panel), and the absolute deviation between these two solutions for several values of $N$ (right panel).</figcaption></figure>

<table><caption>**Table 4.** Data associated with Fig. 9 (static pressurized case).</caption><thead><tr><th></th><th>Exact</th><th>This work</th></tr></thead><tbody><tr><td>$M/\rho_c R_e^3$</td><td>0.60718</td><td>0.60720</td></tr><tr><td>$I/\rho_c R_e^5$</td><td>0.15262</td><td>0.15262</td></tr><tr><td>Steps for convergence</td><td>0</td><td>21</td></tr><tr><td>Virial parameter</td><td>0</td><td>2 &times; 10<sup>-5</sup></td></tr></tbody></table>

cutoff for the mass density at $\hat{\rho}_{amb}$, and then to a mass density jump at $\bar{\omega} = 1$. In the case of a single domain object ($K = 1$ and $\bar{\omega}_1 = 1$ here), we have

$$ \hat{\rho}(\bar{\omega}) = \hat{\rho}_1(\bar{\omega})H(1 - \bar{\omega}), \quad (26) $$

where $\hat{\rho}_1(1) = \hat{\rho}_{amb}$. An interesting test case for ambient pressure is the $n = 5$ polytrope. The analytical solution in the nonrotating case is due to Schuster (e.g., Horedt 2004; see note 3). In the present context, the solution must be truncated at the right isobar $P(1) = P_{amb}$. The analytical solution becomes

$$ \hat{\rho}(\bar{\omega}) = \frac{1}{[1 + \bar{\omega}^2(\hat{\rho}_{amb}^{-2/5} - 1)]^{5/2}}, \quad (27) $$

which verifies $\hat{\rho}(0) = 1$ and $\hat{\rho}(1) = \hat{\rho}_{amb}$. We computed the solution from the SSB approximation. We note that the method is not appropriate for the case $\hat{\rho}_{amb} = 0$ because the mass and especially the radius are infinite for $P_{amb} = 0$. The reconstructed mass density map is displayed in Fig. 9, and the main data are listed in Table 4. As for the $n = 1$ static case, we see that the agreement is excellent and the deviation from the exact solution depends just on the resolution, as expected. We note that the virial parameter was adapted to the context of an ambient pressure. We have

$$ VP = W + 2T + U - U_{amb}, \quad (28) $$

where $U_{amb} = 3P_{amb}V$ (see, e.g., Cox & Giuli 1968).

In the rotating case, Schuster's solution is not adapted anymore, and we switch back to the DROP code for the reference. For this illustration we consider the configuration given in Table 5. It is an intermediate rotator with the same index and ambient pressure as before. A comparison between the full 2D solution and the reconstructed density field, as well as the difference between

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Same as Fig. 3, but for $\bar{\epsilon}_s = 0.90$, $n = 5.0$, and $\hat{\rho}_{amb} = 0.04$ (configuration C); see also Table 5.</figcaption></figure>

<table><caption>**Table 5.** Same caption and same conditions as for Table 2, but for the rotating case with ambient pressure.</caption><thead><tr><th>Configuration C</th><th>DROP<sup>(†)</sup> code</th><th>This work</th></tr></thead><tbody><tr><td>$\bar{\epsilon}_s$</td><td>$\leftarrow$ 0.900</td><td>$\leftarrow$ 0.900</td></tr><tr><td>$n$</td><td>$\leftarrow$ 5</td><td>$\leftarrow$ 5</td></tr><tr><td>$\hat{\rho}_{amb}$</td><td>$\leftarrow$ 0.04</td><td>$\leftarrow$ 0.04</td></tr><tr><td>$M/(\rho_c R_e^3)$</td><td>$5.2781 \times 10^{-1}$</td><td>$5.2871 \times 10^{-1}$</td></tr><tr><td>$I/(\rho_c R_e^5)$</td><td>$1.3129 \times 10^{-1}$</td><td>$1.3169 \times 10^{-2}$</td></tr><tr><td>$V/R_e^3$</td><td>$3.7596$</td><td>$3.7699$</td></tr><tr><td>$\hat{\Omega}^2$</td><td>$7.7911 \times 10^{-2}$</td><td>$7.8114 \times 10^{-2}$<br>($\star$)7.8109 $\times 10^{-2}$</td></tr><tr><td>$J/(G\rho_c^3 R_e^{10})^{1/2}$</td><td>$3.6648 \times 10^{-2}$</td><td>$3.6805 \times 10^{-2}$</td></tr><tr><td>$-W/(G\rho_c^2 R_e^5)$</td><td>$2.2649 \times 10^{-1}$</td><td>$2.2709 \times 10^{-1}$</td></tr><tr><td>$T/(G\rho_c^2 R_e^5)$</td><td>$5.1147 \times 10^{-3}$</td><td>$5.1434 \times 10^{-3}$</td></tr><tr><td>$U/(G\rho_c^5 R_e^5)$</td><td>$2.7139 \times 10^{-1}$</td><td>$2.7190 \times 10^{-1}$</td></tr><tr><td>$U_{amb}/(G\rho_c^2 R_e^5)$</td><td>$5.5118 \times 10^{-2}$</td><td>$5.5283 \times 10^{-2}$</td></tr><tr><td>$[VP/W]$</td><td>$5 \times 10^{-5}$</td><td>$8 \times 10^{-4}$</td></tr></tbody></table>

the two maps, is shown in Fig. 10. We see that the RMS error<fnref n="5" /> is on the order of $10^{-4}$. The main global quantities are given in Table 5. The agreement is once again very good, with local deviations on the order of $10^{-4}$ on the density profile. For global properties, the relative deviations are on the order of $10^{-3}$.

<footnote n="5">⁵ As the true boundary of the fluid is slightly sub-elliptical, some nodes are out of the fluid in the 2D problem, while $\hat{\rho} = \hat{\rho}_{amb}$. To avoid an excessive overestimation of the RMS error, we subtracted the ambient mass density from both maps to obtain this value.</footnote>
<figure><img src="image_11.png" /><figcaption>**Fig. 11.** Main outputs for the Earth using the PREM. From left to right: normalized mass density of the PREM (Dziewonski & Anderson 1981) which is used as an input in this work, and the gradient of the eccentricity squared, squared eccentricity, and normalized pressure determined from the SSB approximation.</figcaption></figure>

### 5.3. The Earth as an example

Dziewonski & Anderson (1981) obtained a seismologic model for the internal structure of the Earth, named the Preliminary Reference Earth Model (PREM). They derived profiles for the mass density, pressure, and gravitational field of the Earth under the hypothesis of a purely spherical planet. In this section we use the mass density profile of the PREM as an input to model a rotating Earth with the SSB approximation; it is plotted in Fig. 11a. It is clear that nonspherical effects have an impact on the mass density distribution. However, the Earth is a slow rotor, so the deviation is expected to be on the order of $1 - \bar{\epsilon}_s \sim 3 \times 10^{-3}$ in relative. The PREM is a model with $K = 10$ domains, and as many mass density jumps$^6$. As $\hat{\rho}(\omega)$ is an input, the SSB approximation developed in this work returns to a single fixed-point problem where only step (i) of the cycle is needed (i.e., we iterate only on Eq. (4b)). Only 15 iterations are required by the SCF method (with $N = 1024$ per domain). We show the results for $d\epsilon^2/d\omega$ and $\epsilon^2(\omega)$ in Figs. 11b and c, respectively. We also report the eccentricity profile computed by Ragazzo (2020), who solved Clairaut's equation with the same piece-wise prescription for the mass density. We see that both solutions are in excellent agreement. The pressure profile $P(\omega)$ at equilibrium is visible in Fig. 11d (see Appendix E for this calculation). It compares very well to the PREM values (the central dimensionless value is $\sim 0.782$ compared to $0.784$ for the PREM). The relative deviations do not exceed $\sim 3 \times 10^{-3}$, which is on the order of the flattening of the outermost surface. This departure is not due to the spheroidal approximation made in this work (as $|VP/W| \sim 7 \times 10^{-7}$), but to the fact that the observed mass density profile has been spherically averaged by the PREM and that centrifugal effects have not been taken into account in their calculations of the pressure profile.

The global properties deduced from SSB approximation are reported in Table 6. However, the data in this table must be rescaled prior to any comparison with observational data. For this purpose we chose to use the mass and equatorial radius given by Chambat et al. (2010), and reported in Table 7. We then deduced $\rho_c = 13\,083.8 \, \text{kg m}^{-3}$, whose value remains close to the PREM value of $1.3088 \times 10^4 \, \text{kg m}^{-3}$. Table 7 compares the momentum of inertia, rotation rate, and first two gravitational moments. We see that the normalized moment of inertia and the rotation rates are in good agreement with the observational data, again within $\sim 3 \times 10^{-3}$ in relative. Regarding the gravitational

<figure><img src="image_12.png" /><figcaption>**Table 6.** Output dimensionless quantities obtained for the rotating Earth from the SSB approximation, for $N = 10240$ in total.</figcaption></figure>

<table><caption>**Table 7.** Physical properties of the Earth rescaled; $R_v = R_e \bar{\epsilon}_s^{1/3}$ is the mean volumetric radius ($\approx 6.371 \times 10^6$ m).</caption><thead><tr><th></th><th>This work</th></tr></thead><tbody><tr><td>$\bar{\epsilon}_s$</td><td>$\leftarrow 0.99665$</td></tr><tr><td>$M/(\rho_c R_e^3)$</td><td>$1.7592$</td></tr><tr><td>$I/(\rho_c R_e^5)$</td><td>$5.8188 \times 10^{-1}$</td></tr><tr><td>$\hat{\Omega}^2$</td><td>$6.1199 \times 10^{-3}$</td></tr><tr><td>$\hat{P}_c$</td><td>$7.8200 \times 10^{-1}$</td></tr><tr><td>$J/(G\rho_c^3 R_e^{10})^{1/2}$</td><td>$4.5521 \times 10^{-2}$</td></tr><tr><td>$-W/(G\rho_c^2 R_e^5)$</td><td>$2.0631$</td></tr><tr><td>$T/(G\rho_c^2 R_e^5)$</td><td>$1.7805 \times 10^{-3}$</td></tr><tr><td>$U/(G\rho_c^2 R_e^5)$</td><td>$2.0596$</td></tr><tr><td>$|VP/W|$</td><td>$7 \times 10^{-7}$</td></tr></tbody></table>

<table><thead><tr><th>Data</th><th>Observed<sup>(†)</sup></th><th>Second-order<br>Clairaut<sup>(†)</sup></th><th>This work</th></tr></thead><tbody><tr><td>$M$ [kg]</td><td>$5.97218 \times 10^{24}$</td><td></td><td></td></tr><tr><td>$R_e$ [m]</td><td>$6.378137 \times 10^6$</td><td></td><td></td></tr><tr><td>$I/MR_v^2$</td><td>$3.3069 \times 10^{-1}$</td><td></td><td>$3.3151 \times 10^{-1}$</td></tr><tr><td>$\Omega$ [s<sup>-1</sup>]</td><td>$7.2921 \times 10^{-5}$</td><td></td><td>$7.3104 \times 10^{-5}$</td></tr><tr><td>$J_2$</td><td>$1.0826 \times 10^{-3}$</td><td>$1.0712 \times 10^{-3}$</td><td>$1.0771 \times 10^{-3}$</td></tr><tr><td>$-J_4$</td><td>$1.620 \times 10^{-6}$</td><td>$2.96 \times 10^{-6}$</td><td>$2.8233 \times 10^{-6}$</td></tr></tbody></table>

Notes. $^{(†)}$See Chambat et al. (2010) and references therein.

moments, we see that $J_2$ is close to the observed value, roughly within 1%. The result is worse for $J_4$, a deviation reaching 75%. A similar discrepancy is found by Chambat et al. (2010), who solved the second-order Clairaut's equation with a different (yet similar) mass density profile for the Earth. As quoted by these authors, this deviation from the observed gravitational moments is due to nonhydrostatic characteristics of the Earth, which are not taken into account in our equation set.

## 6. Summary

In this article we showed that the 2D structure of a rotating self-gravitating fluid can be determined with good precision from

<footnote>$^6$As the atmosphere is not taken into account in the PREM, a mass density jump is present at the outermost surface: from liquid water to the outer space.</footnote>
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Main outputs for the Earth using the PREM. From left to right: normalized mass density of the PREM (Dziewonski & Anderson 1981) which is used as an input in this work, and the gradient of the eccentricity squared, squared eccentricity, and normalized pressure determined from the SSB approximation.</figcaption></figure>

### 5.3. The Earth as an example

Dziewonski & Anderson (1981) obtained a seismologic model for the internal structure of the Earth, named the Preliminary Reference Earth Model (PREM). They derived profiles for the mass density, pressure, and gravitational field of the Earth under the hypothesis of a purely spherical planet. In this section we use the mass density profile of the PREM as an input to model a rotating Earth with the SSB approximation; it is plotted in Fig. 11a. It is clear that nonspherical effects have an impact on the mass density distribution. However, the Earth is a slow rotor, so the deviation is expected to be on the order of $1 - \bar{\epsilon}_s \sim 3 \times 10^{-3}$ in relative. The PREM is a model with $K = 10$ domains, and as many mass density jumps$^6$. As $\hat{\rho}(\omega)$ is an input, the SSB approximation developed in this work returns to a single fixed-point problem where only step (i) of the cycle is needed (i.e., we iterate only on Eq. (4b)). Only 15 iterations are required by the SCF method (with $N = 1024$ per domain). We show the results for $d\epsilon^2/d\omega$ and $\epsilon^2(\omega)$ in Figs. 11b and c, respectively. We also report the eccentricity profile computed by Ragazzo (2020), who solved Clairaut's equation with the same piece-wise prescription for the mass density. We see that both solutions are in excellent agreement. The pressure profile $P(\omega)$ at equilibrium is visible in Fig. 11d (see Appendix E for this calculation). It compares very well to the PREM values (the central dimensionless value is ~0.782 compared to 0.784 for the PREM). The relative deviations do not exceed ~3 $\times 10^{-3}$, which is on the order of the flattening of the outermost surface. This departure is not due to the spheroidal approximation made in this work (as $|VP/W| \sim 7 \times 10^{-7}$), but to the fact that the observed mass density profile has been spherically averaged by the PREM and that centrifugal effects have not been taken into account in their calculations of the pressure profile.

The global properties deduced from SSB approximation are reported in Table 6. However, the data in this table must be rescaled prior to any comparison with observational data. For this purpose we chose to use the mass and equatorial radius given by Chambat et al. (2010), and reported in Table 7. We then deduced $\rho_c = 13\,083.8 \, \text{kg m}^{-3}$, whose value remains close to the PREM value of $1.3088 \times 10^4 \, \text{kg m}^{-3}$. Table 7 compares the momentum of inertia, rotation rate, and first two gravitational moments. We see that the normalized moment of inertia and the rotation rates are in good agreement with the observational data, again within ~3 $\times 10^{-3}$ in relative. Regarding the gravitational

<figure><img src="image_13.png" /><figcaption>**Table 6.** Output dimensionless quantities obtained for the rotating Earth from the SSB approximation, for $N = 10240$ in total.</figcaption></figure>

<table><caption>**Table 7.** Physical properties of the Earth rescaled; $R_v = R_e \bar{\epsilon}_s^{1/3}$ is the mean volumetric radius ($\approx 6.371 \times 10^6$ m).</caption><thead><tr><th></th><th>This work</th></tr></thead><tbody><tr><td>$\bar{\epsilon}_s$</td><td>$\leftarrow 0.99665$</td></tr><tr><td>$M/(\rho_c R_e^3)$</td><td>$1.7592$</td></tr><tr><td>$I/(\rho_c R_e^5)$</td><td>$5.8188 \times 10^{-1}$</td></tr><tr><td>$\hat{\Omega}^2$</td><td>$6.1199 \times 10^{-3}$</td></tr><tr><td>$\hat{P}_c$</td><td>$7.8200 \times 10^{-1}$</td></tr><tr><td>$J/(G\rho_c^3 R_e^{10})^{1/2}$</td><td>$4.5521 \times 10^{-2}$</td></tr><tr><td>$-W/(G\rho_c^2 R_e^5)$</td><td>$2.0631$</td></tr><tr><td>$T/(G\rho_c^2 R_e^5)$</td><td>$1.7805 \times 10^{-3}$</td></tr><tr><td>$U/(G\rho_c^2 R_e^5)$</td><td>$2.0596$</td></tr><tr><td>$|VP/W|$</td><td>$7 \times 10^{-7}$</td></tr></tbody></table>

<table><thead><tr><th>Data</th><th>Observed<sup>(†)</sup></th><th>Second-order<br>Clairaut<sup>(†)</sup></th><th>This work</th></tr></thead><tbody><tr><td>$M$ [kg]</td><td>$5.97218 \times 10^{24}$</td><td></td><td></td></tr><tr><td>$R_e$ [m]</td><td>$6.378137 \times 10^6$</td><td></td><td></td></tr><tr><td>$I/MR_v^2$</td><td>$3.3069 \times 10^{-1}$</td><td></td><td>$3.3151 \times 10^{-1}$</td></tr><tr><td>$\Omega$ [s<sup>-1</sup>]</td><td>$7.2921 \times 10^{-5}$</td><td></td><td>$7.3104 \times 10^{-5}$</td></tr><tr><td>$J_2$</td><td>$1.0826 \times 10^{-3}$</td><td>$1.0712 \times 10^{-3}$</td><td>$1.0771 \times 10^{-3}$</td></tr><tr><td>$-J_4$</td><td>$1.620 \times 10^{-6}$</td><td>$2.96 \times 10^{-6}$</td><td>$2.8233 \times 10^{-6}$</td></tr></tbody></table>

Notes. (†)See Chambat et al. (2010) and references therein.

moments, we see that $J_2$ is close to the observed value, roughly within 1%. The result is worse for $J_4$, a deviation reaching 75%. A similar discrepancy is found by Chambat et al. (2010), who solved the second-order Clairaut's equation with a different (yet similar) mass density profile for the Earth. As quoted by these authors, this deviation from the observed gravitational moments is due to nonhydrostatic characteristics of the Earth, which are not taken into account in our equation set.

## 6. Summary

In this article we showed that the 2D structure of a rotating self-gravitating fluid can be determined with good precision from

<footnote>$^6$As the atmosphere is not taken into account in the PREM, a mass density jump is present at the outermost surface: from liquid water to the outer space.</footnote>
the theory of Nested Spheroidal Figures of Equilibrium, which assumes that isopycnics are perfect spheroids. The method uses the general IDE established in Staelen & Huré (2024), connecting the local eccentricity $ε$ of isopycnics to the local mass density $ρ$. As this equation involves only the equatorial radius, it can be coupled to the centrifugal balance along the polar axis, providing a specific one-dimensional projection of the genuine bi-dimensional problem. As shown, the new problem consists in solving two fixed-points problems coupled together. From the equatorial solution {$ρ(a), ε(a)$}, which can be determined by a simple SCF method, we recover the full structure by propagating the mass density along the ellipse, from the equator to the pole, and for any equatorial radius $a \in [0, R_e]$. We provided a series of examples supporting the efficiency and versatility of this SSB approximation. In particular, the method is well suited to systems made of heterogeneous domains separated by mass density jumps. It can account for an ambient pressure. The method is also flexible in terms of barotropic EOS. Here we used a polytrope, but any kind of $P(ρ)$-relation is usable. The main results of this paper can be summarized as follows:

- The SSB approximation is exact for nonrotating gas, because the NSFoE is also exact in this case (all confocal parameters are null). The method is then equivalent to a standard Lane-Emden solver. The precision is then fully governed by the numerical schemes implemented.

- For slow rotators (i.e., small flattenings), the method has an excellent precision. Depending on the EOS, the solution can reproduce the full 2D mass density profile, typically within $10^{-5}$ in absolute (dimensionless profile), even at low to moderate numerical resolution.

- For fast rotators (i.e., large flattenings), the SSB approximation furnishes good and reliable results, whatever the EOS.

- For hard EOSs (typically $0 < n < 1$), the deviation of the true surface from a spheroid is very small, but there is a significant amount of matter close to the surface. The precision of SSB approximation is very good provided the system stays far from the critical rotation state. Near the sequence ends, the precision is acceptable, with typically 1% in the mass density (RMS value).

- For soft EOSs ($n > 1$) the precision is very good because the amount of matter located close to the surface has a negligible contribution to the total mass (and gravitational potential), although the deviation of the true surface from a spheroid is significant.

- Compared to the full 2D problem in which the surface and isopycnics are not constrained, the mass density in the vicinity of the center computed with the SSB approximation has an excellent precision. This is also true for values close to the pole and to the equator.

- The computing time is reduced by at least two orders of magnitude with respect to the full 2D problem. This is a direct consequence of dimension reduction. This enables us to reach very high numerical resolutions in a short time, which is particularly attractive for generating grids of models for instance.

**Acknowledgements.** We are grateful to F. Chambat for informed comments on the paper before submission. We thank the anonymous referee for suggestions to improve the paper and valuable references, in particular on the state-of-art regarding stellar models.

## References

Abramyan, M. G., & Kaplan, S. A. 1974, Astrophysics, 10, 358

Amedt, P., Lanza, A., & Abramowicz, M. A. 1989, ApJ, 343, 437

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

Binney, J., & Tremaine, S. 1987, Galactic dynamics (Princeton: Princeton University Press), 747

Chambat, F., Ricard, Y., & Valette, B. 2010, Geophys. J. Int., 183, 727

Chandrasekhar, S. 1933, MNRAS, 93, 390

Chandrasekhar, S. 1969, Ellipsoidal Figures of Equilibrium (London: Yale Univ. Press)

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

Cisneros Parra, J. U., Martínez Herrera, F. J., & Montalvo Castro, J. D. 2015, Rev. Mex. Astron. Astrofis., 51, 121

Cox, J. P., & Giuli, R. T. 1968, Principles of Stellar Structure Volume I: Physical Principles, 1st edn. (Philadelphia: Gordon and Breach)

Debras, F., & Chabrier, G. 2018, A&A, 609, A97

Domiciano de Souza, A., Kervella, P., Moser Faes, D., et al. 2014, A&A, 569, A10

Dziewonski, A. M., & Anderson, D. L. 1981, Phys. Earth Planet. Inter., 25, 297

Espinosa Lara, F., & Rieutord, M. 2013, A&A, 552, A35

Even, W., & Tohline, J. E. 2009, ApJ, 184, 248

Gradshteyn, I. S., & Ryzhik, I. M. 2014, in Table of Integrals, Series, and Products, 8th edn., eds. D. Zwillinger, & V. Moll (Boston: Academic Press), 776

Hachisu, I. 1986a, ApJS, 61, 479

Hachisu, I. 1986b, ApJS, 62, 461

Hadjifotinou, K. G. 2000, A&A, 354, 328

Horedt, G. P. 2004, in Polytropes - Applications in Astrophysics and Related Fields, Astrophys. Space Sci. Lib., 306

Houdayer, P. S., & Reese, D. R. 2023, A&A, 675, A181

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

Huré, J.-M., Hersant, F., & Nasello, G. 2018, MNRAS, 475, 63

James, R. A. 1964, ApJ, 140, 552

Kadam, K., Motl, P. M., Frank, J., Clayton, G. C., & Marcello, D. C. 2016, MNRAS, 462, 2237

Kiuchi, K., Nagakura, H., & Yamada, S. 2010, ApJ, 717, 666

Kluźniak, W., & Rosińska, D. 2013, MNRAS, 434, 2825

Kovetz, A. 1968, ApJ, 154, 999

Lander, S. K., & Jones, D. I. 2009, MNRAS, 395, 2162

Liu, F. K. 1996, MNRAS, 281, 1197

Mach, P. 2012, J. Math. Phys., 53, 062503

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

Mishra, B., & Vaidya, B. 2015, MNRAS, 447, 1154

Montalvo, D., Martinez, F. J., & Cisneros, J. 1983, Rev. Mex. Astron. Astrofis., 5, 293

Nettelmann, N. 2017, A&A, 606, A139

Odrzywołek, A. 2003, MNRAS, 345, 497

Ostriker, J. P., & Mark, J. W.-K. 1968, ApJ, 151, 1075

Ragazzo, C. 2020, São Paulo Journal of Mathematical Sciences, 14, 1

Rampalli, R., Smock, A., Newton, E. R., Daniel, K. J., & Curtis, J. L. 2023, ApJ, 958, 76

Rieutord, M. 2006, in SF2A-2006: Semaine de l'Astrophysique Francaise, eds. D. Barret, F. Casoli, G. Lagache, A. Lecavelier, & L. Pagani, 501

Rieutord, M., Espinosa Lara, F., & Putigny, B. 2016, J. Comput. Phys., 318, 277

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

Roxburgh, I. W., & Strittmatter, P. A. 1966, MNRAS, 133, 345

Seidov, Z. F., & Kuzakhmedov, R. K. 1978, Sov. Astron., 22, 711

Sharma, V. D. 1977, Phys. Lett. A, 60, 381

Srivastava, S. 1962, ApJ, 136, 680

Staelen, C., & Huré, J.-M. 2024, MNRAS, 527, 863 (Paper IV)

Tassoul, J.-L. 1978, Theory of Rotating Stars (Princeton: Princeton University Press)

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)

Tohline, J. E. 2021, a (MediaWiki-based) Vistrails.org publication, https://www.vistrails.org/index.php/User:Tohline

Tomimura, Y., & Eriguchi, Y. 2005, MNRAS, 359, 1117

Venditti, F. C. F., de Almeida Junior, A. K., & Prado, A. F. B. A. 2020, Planet. Space Sci., 192, 105063

Véronet, A. 1912, Journal de mathématiques pures et appliquées 6e série, 8, 331

Zharkov, V. N., & Trubitsyn, V. P. 1970, Sov. Astron., 13, 981
the theory of Nested Spheroidal Figures of Equilibrium, which assumes that isopycnics are perfect spheroids. The method uses the general IDE established in Staelen & Huré (2024), connecting the local eccentricity $\epsilon$ of isopycnics to the local mass density $\rho$. As this equation involves only the equatorial radius, it can be coupled to the centrifugal balance along the polar axis, providing a specific one-dimensional projection of the genuine bi-dimensional problem. As shown, the new problem consists in solving two fixed-points problems coupled together. From the equatorial solution {$\rho(a), \epsilon(a)$}, which can be determined by a simple SCF method, we recover the full structure by propagating the mass density along the ellipse, from the equator to the pole, and for any equatorial radius $a \in [0, R_e]$. We provided a series of examples supporting the efficiency and versatility of this SSB approximation. In particular, the method is well suited to systems made of heterogeneous domains separated by mass density jumps. It can account for an ambient pressure. The method is also flexible in terms of barotropic EOS. Here we used a polytrope, but any kind of $P(\rho)$-relation is usable. The main results of this paper can be summarized as follows:

- The SSB approximation is exact for nonrotating gas, because the NSFoE is also exact in this case (all confocal parameters are null). The method is then equivalent to a standard Lane-Emden solver. The precision is then fully governed by the numerical schemes implemented.

- For slow rotators (i.e., small flattenings), the method has an excellent precision. Depending on the EOS, the solution can reproduce the full 2D mass density profile, typically within $10^{-5}$ in absolute (dimensionless profile), even at low to moderate numerical resolution.

- For fast rotators (i.e., large flattenings), the SSB approximation furnishes good and reliable results, whatever the EOS.

- For hard EOSs (typically $0 < n < 1$), the deviation of the true surface from a spheroid is very small, but there is a significant amount of matter close to the surface. The precision of SSB approximation is very good provided the system stays far from the critical rotation state. Near the sequence ends, the precision is acceptable, with typically 1% in the mass density (RMS value).

- For soft EOSs ($n > 1$) the precision is very good because the amount of matter located close to the surface has a negligible contribution to the total mass (and gravitational potential), although the deviation of the true surface from a spheroid is significant.

- Compared to the full 2D problem in which the surface and isopycnics are not constrained, the mass density in the vicinity of the center computed with the SSB approximation has an excellent precision. This is also true for values close to the pole and to the equator.

- The computing time is reduced by at least two orders of magnitude with respect to the full 2D problem. This is a direct consequence of dimension reduction. This enables us to reach very high numerical resolutions in a short time, which is particularly attractive for generating grids of models for instance.

**Acknowledgements.** We are grateful to F. Chambat for informed comments on the paper before submission. We thank the anonymous referee for suggestions to improve the paper and valuable references, in particular on the state-of-art regarding stellar models.

## References

- Abramyan, M. G., & Kaplan, S. A. 1974, Astrophysics, 10, 358

- Amendt, P., Lanza, A., & Abramowicz, M. A. 1989, ApJ, 343, 437

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

- Binney, J., & Tremaine, S. 1987, Galactic dynamics (Princeton University Press), 747

- Chambat, F., Ricard, Y., & Valette, B. 2010, Geophys. J. Int., 183, 727

- Chandrasekhar, S. 1933, MNRAS, 93, 390

- Chandrasekhar, S. 1969, Ellipsoidal Figures of Equilibrium (London: Yale Univ. Press)

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

- Cisneros Parra, J. U., Martínez Herrera, F. J., & Montalvo Castro, J. D. 2015, Rev. Mex. Astron. Astrofis., 51, 121

- Cox, J. P., & Giuli, R. T. 1968, Principles of Stellar Structure Volume I: Physical Principles, 1st edn. (Philadelphia: Gordon and Breach)

- Debras, F., & Chabrier, G. 2018, A&A, 609, A97

- Domiciano de Souza, A., Kervella, P., Moser Faes, D., et al. 2014, A&A, 569, A10

- Dziewonski, A. M., & Anderson, D. L. 1981, Phys. Earth Planet. Inter., 25, 297

- Espinosa Lara, F., & Rieutord, M. 2013, A&A, 552, A35

- Even, W., & Tohline, J. E. 2009, ApJ, 184, 248

- Gradshteyn, I. S., & Ryzhik, I. M. 2014, in Table of Integrals, Series, and Products, 8th edn., eds. D. Zwillinger, & V. Moll (Boston: Academic Press), 776

- Hachisu, I. 1986a, ApJS, 61, 479

- Hachisu, I. 1986b, ApJS, 62, 461

- Hadjifotinou, K. G. 2000, A&A, 354, 328

- Horedt, G. P. 2004, in Polytropes - Applications in Astrophysics and Related Fields, Astrophys. Space Sci. Lib., 306

- Houdayer, P. S., & Reese, D. R. 2023, A&A, 675, A181

- 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

- Huré, J.-M., Hersant, F., & Nasello, G. 2018, MNRAS, 475, 63

- James, R. A. 1964, ApJ, 140, 552

- Kadam, K., Motl, P. M., Frank, J., Clayton, G. C., & Marcello, D. C. 2016, MNRAS, 462, 2237

- Kiuchi, K., Nagakura, H., & Yamada, S. 2010, ApJ, 717, 666

- Kluźniak, W., & Rosińska, D. 2013, MNRAS, 434, 2825

- Kovetz, A. 1968, ApJ, 154, 999

- Lander, S. K., & Jones, D. I. 2009, MNRAS, 395, 2162

- Liu, F. K. 1996, MNRAS, 281, 1197

- Mach, P. 2012, J. Math. Phys., 53, 062503

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

- Mishra, B., & Vaidya, B. 2015, MNRAS, 447, 1154

- Montalvo, D., Martinez, F. J., & Cisneros, J. 1983, Rev. Mex. Astron. Astrofis., 5, 293

- Nettelmann, N. 2017, A&A, 606, A139

- Odrzywołek, A. 2003, MNRAS, 345, 497

- Ostriker, J. P., & Mark, J. W.-K. 1968, ApJ, 151, 1075

- Ragazzo, C. 2020, São Paulo Journal of Mathematical Sciences, 14, 1

- Rampalli, R., Smock, A., Newton, E. R., Daniel, K. J., & Curtis, J. L. 2023, ApJ, 958, 76

- Rieutord, M. 2006, in SF2A-2006: Semaine de l'Astrophysique Francaise, eds. D. Barret, F. Casoli, G. Lagache, A. Lecavelier, & L. Pagani, 501

- Rieutord, M., Espinosa Lara, F., & Putigny, B. 2016, J. Comput. Phys., 318, 277

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

- Roxburgh, I. W., & Strittmatter, P. A. 1966, MNRAS, 133, 345

- Seidov, Z. F., & Kuzakhmedov, R. K. 1978, Sov. Astron., 22, 711

- Sharma, V. D. 1977, Phys. Lett. A, 60, 381

- Srivastava, S. 1962, ApJ, 136, 680

- Staelen, C., & Huré, J.-M. 2024, MNRAS, 527, 863 (Paper IV)

- Tassoul, J.-L. 1978, Theory of Rotating Stars (Princeton: Princeton University Press)

- 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)

- Tohline, J. E. 2021, a (MediaWiki-based) Vistrails.org publication, https://www.vistrails.org/index.php/User:Tohline

- Tomimura, Y., & Eriguchi, Y. 2005, MNRAS, 359, 1117

- Venditti, F. C. F., de Almeida Junior, A. K., & Prado, A. F. B. A. 2020, Planet. Space Sci., 192, 105063

- Véronet, A. 1912, Journal de mathématiques pures et appliquées 6e série, 8, 331

- Zharkov, V. N., & Trubitsyn, V. P. 1970, Sov. Astron., 13, 981
# Appendix A: The functions $\chi, \mu, \nu$, and $\eta$

We define for convenience

$$e^2(c) = \frac{\omega'^2 \varepsilon^2(\omega')}{\omega^2[1+c]}, \qquad (A.1)$$

and $\bar{e}(c) = \sqrt{1 - e^2(c)} \ge 0$, where the continuous confocal parameter $c \equiv c(\omega', \omega)$ is given by

$$c(\omega', \omega) = \frac{\omega'^2 \varepsilon^2(\omega')}{\omega^2} - \varepsilon^2(\omega). \qquad (A.2)$$

In these conditions, and given the function $\mathcal{A}(\bar{\varepsilon})$ defined by Eq. (B.1), $\chi$ is defined as

$$\chi(\omega', \omega; \varepsilon) = \begin{cases} \frac{\omega'^3 \bar{\varepsilon}(\omega')}{\omega^4} \left[ \mathcal{A}(\bar{\varepsilon}(0)) - \frac{\mathcal{A}(\bar{\varepsilon}(c))}{\sqrt{1+c}} \right], & \omega' < \omega \\ 0, & \omega' \ge \omega \end{cases} \qquad (A.3)$$

$\mu$ is defined as

$$\mu(\omega', \omega; \varepsilon) = \begin{cases} \frac{\omega' \bar{\varepsilon}(\omega') \mathcal{A}(\bar{\varepsilon}(c))}{\omega \varepsilon^2(\omega') \sqrt{1+c}} - \frac{\omega' \bar{\varepsilon}(\omega')}{\omega \varepsilon^2(\omega') \bar{\varepsilon}(\omega)}, & \omega' < \omega \\ \frac{\bar{\varepsilon}(\omega')}{\varepsilon^2(\omega')} \mathcal{A}(\bar{\varepsilon}(\omega')) - \frac{1}{\varepsilon^2(\omega')}, & \omega' \ge \omega \end{cases}, \qquad (A.4)$$

$\eta$ is defined as

$$\eta(\omega', \omega; \varepsilon) = \begin{cases} \omega' \omega \frac{\bar{\varepsilon}(\omega')}{\varepsilon^2(\omega')} \left[ \mathcal{A}(\bar{\varepsilon}(c)) \sqrt{1+c} - \bar{\varepsilon}(\omega) \right], & \omega' < \omega \\ \omega^2 [1+c] \frac{\bar{\varepsilon}(\omega')}{\varepsilon^2(\omega')} \mathcal{A}(\bar{\varepsilon}(\omega')) - \omega^2 \frac{\bar{\varepsilon}^2(\omega)}{\varepsilon^2(\omega')}, & \omega' \ge \omega \end{cases}, \qquad (A.5)$$

and $\kappa$ is defined as

$$\kappa(\omega', \omega; \varepsilon) = \begin{cases} \frac{\omega' \bar{\varepsilon}(\omega')}{\omega \varepsilon^2(\omega')} \left\{ [1 - 2e^2(0)] \mathcal{A}(\bar{\varepsilon}(0)) - 2\bar{\varepsilon}(\omega) - \bar{\varepsilon}(0) + 2\mathcal{A}(\bar{\varepsilon}(c)) \sqrt{1+c} \right\}, & \omega' < \omega \\ \frac{3 - 2\varepsilon^2(\omega)}{\varepsilon^2(\omega')} [\bar{\varepsilon}(\omega') \mathcal{A}(\bar{\varepsilon}(\omega')) - 1] + 1, & \omega' \ge \omega \end{cases}. \qquad (A.6)$$

# Appendix B: A single function for the prolate and oblate cases

As can be seen in Paper IV, which is devoted to oblate configurations, all functions $\kappa, \chi, \mu$, (and $\eta$) depend on $\varepsilon$, which appears in the argument of an arcsin function. During the numerical tests, we deal with transition states where the deep isopycnic surfaces are slightly prolate, resulting in a purely imaginary eccentricity. This naturally creates numerical issues, which can be circumvented by a simple prolongation since $\arcsin(z) = \arcsinh(i z)/i$, where $i$ is the imaginary unit and $z$ is a complex number (see, e.g., Gradshteyn & Ryzhik 2014). So, for the numerical computation we define

$$\mathcal{A}(\bar{\varepsilon}) = \begin{cases} \frac{\arcsin(\sqrt{1 - \bar{\varepsilon}^2})}{\sqrt{1 - \bar{\varepsilon}^2}}, & 0 < \bar{\varepsilon} < 1, \\ 1, & \bar{\varepsilon} = 1, \\ \frac{\arcsinh(\sqrt{\bar{\varepsilon}^2 - 1})}{\sqrt{\bar{\varepsilon}^2 - 1}}, & \bar{\varepsilon} > 1, \end{cases} \qquad (B.1)$$

which enables us to make the transition from oblate to prolate configurations. It can be verified (see Appendix A), by introducing this function $\mathcal{A}$ in all integrand kernels defined in Paper IV, that only real terms like $\bar{\varepsilon}$ or $\varepsilon^2$ finally survive, whatever the case, either oblate or prolate.
# Appendix A: The functions $\chi$, $\mu$, $\nu$, and $\eta$

We define for convenience

$$e^2(c) = \frac{\omega'^2 \varepsilon^2(\omega')}{\omega^2[1+c]}, \qquad (A.1)$$

and $\bar{e}(c) = \sqrt{1 - e^2(c)} \ge 0$, where the continuous confocal parameter $c \equiv c(\omega', \omega)$ is given by

$$c(\omega', \omega) = \frac{\omega'^2 \varepsilon^2(\omega')}{\omega^2} - \varepsilon^2(\omega). \qquad (A.2)$$

In these conditions, and given the function $\mathcal{A}(\bar{\varepsilon})$ defined by Eq. (B.1), $\chi$ is defined as

$$\chi(\omega', \omega; \varepsilon) = \begin{cases} \frac{\omega'^3 \bar{\varepsilon}(\omega')}{\omega^4} \left[ \mathcal{A}(\bar{\varepsilon}(0)) - \frac{\mathcal{A}(\bar{\varepsilon}(c))}{\sqrt{1+c}} \right], & \omega' < \omega \\ 0, & \omega' \ge \omega \end{cases} \qquad (A.3)$$

$\mu$ is defined as

$$\mu(\omega', \omega; \varepsilon) = \begin{cases} \frac{\omega' \bar{\varepsilon}(\omega') \mathcal{A}(\bar{\varepsilon}(c))}{\omega \varepsilon^2(\omega') \sqrt{1+c}} - \frac{\omega' \bar{\varepsilon}(\omega')}{\omega \varepsilon^2(\omega') \bar{\varepsilon}(\omega)}, & \omega' < \omega \\ \frac{\bar{\varepsilon}(\omega')}{\varepsilon^2(\omega')} \mathcal{A}(\bar{\varepsilon}(\omega')) - \frac{1}{\varepsilon^2(\omega')}, & \omega' \ge \omega \end{cases}, \qquad (A.4)$$

$\eta$ is defined as

$$\eta(\omega', \omega; \varepsilon) = \begin{cases} \omega' \omega \frac{\bar{\varepsilon}(\omega')}{\varepsilon^2(\omega')} \left[ \mathcal{A}(\bar{\varepsilon}(c)) \sqrt{1+c} - \bar{\varepsilon}(\omega) \right] & \omega' < \omega \\ \omega'^2 [1+c] \frac{\bar{\varepsilon}(\omega')}{\varepsilon^2(\omega')} \mathcal{A}(\bar{\varepsilon}(\omega')) - \omega'^2 \frac{\bar{\varepsilon}^2(\omega)}{\varepsilon^2(\omega')}, & \omega' \ge \omega \end{cases}, \qquad (A.5)$$

and $\kappa$ is defined as

$$\kappa(\omega', \omega; \varepsilon) = \begin{cases} \frac{\omega' \bar{\varepsilon}(\omega')}{\omega \varepsilon^2(\omega')} \left\{ [1 - 2e^2(0)] \mathcal{A}(\bar{\varepsilon}(0)) - 2\bar{\varepsilon}(\omega) - \bar{\varepsilon}(0) + 2\mathcal{A}(\bar{\varepsilon}(c)) \sqrt{1+c} \right\}, & \omega' < \omega \\ \frac{3 - 2\varepsilon^2(\omega)}{\varepsilon^2(\omega')} [\bar{\varepsilon}(\omega') \mathcal{A}(\bar{\varepsilon}(\omega')) - 1] + 1, & \omega' \ge \omega \end{cases}. \qquad (A.6)$$

# Appendix B: A single function for the prolate and oblate cases

As can be seen in Paper IV, which is devoted to oblate configurations, all functions $\kappa$, $\chi$, $\mu$, (and $\eta$) depend on $\varepsilon$, which appears in the argument of an arcsin function. During the numerical tests, we deal with transition states where the deep isopycnic surfaces are slightly prolate, resulting in a purely imaginary eccentricity. This naturally creates numerical issues, which can be circumvented by a simple prolongation since $\arcsin(z) = \arcsinh(i z)/i$, where $i$ is the imaginary unit and $z$ is a complex number (see, e.g., Gradshteyn & Ryzhik 2014). So, for the numerical computation we define

$$\mathcal{A}(\bar{\varepsilon}) = \begin{cases} \frac{\arcsin(\sqrt{1 - \bar{\varepsilon}^2})}{\sqrt{1 - \bar{\varepsilon}^2}}, & 0 < \bar{\varepsilon} < 1, \\ 1, & \bar{\varepsilon} = 1, \\ \frac{\arcsinh(\sqrt{\bar{\varepsilon}^2 - 1})}{\sqrt{\bar{\varepsilon}^2 - 1}}, & \bar{\varepsilon} > 1, \end{cases} \qquad (B.1)$$

which enables us to make the transition from oblate to prolate configurations. It can be verified (see Appendix A), by introducing this function $\mathcal{A}$ in all integrand kernels defined in Paper IV, that only real terms like $\bar{\varepsilon}$ or $\varepsilon^2$ finally survive, whatever the case, either oblate or prolate.
## Appendix C: The enthalpy profile

Along the polar axis $R = 0$, the Bernoulli equation Eq. (1a) reads

$$H(0, Z) + \Psi(0, Z) = C^{\text{te}}. \qquad (C.1)$$

In order to estimate $\Psi(0, Z)$, we need to go back to Paper II. The gravitational potential of a heterogeneous system made of L homogeneous spheroidal layers is given in that article by Eq. (14). At point $A_j$, at the intersection between the polar axis and the interface between layers $j$ and $j+1$, it reads

$$\frac{\Psi|_{A_j}}{-2\pi G a_j^2} = \sum_{i=1}^{j} (\rho_i - \rho_{i+1}) \frac{\bar{\varepsilon}_i}{\varepsilon_i^3} \left[ (1 + c_{i,j}) \arcsin\left(\frac{q_{i,j} \varepsilon_i}{\sqrt{1 + c_{i,j}}}\right) - q_{i,j} \varepsilon_i \bar{\varepsilon}_j \right] + \sum_{i=1}^{j} (\rho_i - \rho_{i+1}) \frac{\bar{\varepsilon}_i}{\varepsilon_i^3} \left[ (1 + c_{i,j}) \arcsin(\varepsilon_i) (q_{i,j}^2 \varepsilon_i^2 + \bar{\varepsilon}_j^2) - \frac{\varepsilon_i \bar{\varepsilon}_j^2}{\bar{\varepsilon}_i^2} \right], \qquad (C.2)$$

where $a_j$ is the equatorial radius of layer $j$, $\varepsilon_j$ is its eccentricity, $\rho_j$ is its mass density, $\bar{\varepsilon}_j = \sqrt{1 - \varepsilon_j^2}$, $q_{i,j} = a_i/a_j$, and $c_{i,j} = q_{i,j}^2 \varepsilon_i^2 - \varepsilon_j^2$. Oblate spheroids were assumed in Paper II, but it can be generalized to any spheroid by introducing the function $\mathcal{A}$ defined in Appendix B. In the continuous limit $L \to \infty$, we have

$$\frac{\Psi(R=0, Z)}{G \rho_c R_c^2} = 2\pi \int_0^1 d\hat{\rho}(\sigma')\eta(\sigma', \sigma; \varepsilon). \qquad (C.3)$$

As $H(0, Z)$ is a function of $\sigma$ only, we define $\hat{H}(\sigma) = H(0, Z)/(G\rho_c R_c^2)$. By Eq. (C.1), we have

$$\hat{H}(\sigma) + 2\pi \int_{0}^{1} d\hat{\rho}(\sigma')\eta(\sigma', \sigma; \varepsilon) = \hat{H}(1) + 2\pi \int_{0}^{1} d\hat{\rho}(\sigma')\eta(\sigma', 1; \varepsilon), \qquad (C.4)$$

where the constant $C^{\text{te}}$ was evaluated at the pole of the whole object (i.e., at $\sigma = 1$). We then see that Eq. (C.4) thus yields Eq. (4a).

## Appendix D: The spherical limit

For $\varepsilon^2 \ll 1$, we have

$$\mathcal{A}(\bar{\varepsilon}) \approx 1 + \frac{\varepsilon^2}{6} \qquad (D.1)$$

for both $\bar{\varepsilon} < 1$ and $\bar{\varepsilon} > 1$. So we can expand functions $\kappa$, $\chi$, $\mu$, and $\eta$ at zeroth order to obtain these functions at $\bar{\varepsilon} = 1$ (i.e., in the nonrotating spherical case). We have $\kappa(\sigma', \sigma) = 0$, which naturally imposes $\hat{\Omega}^2 = 0$ (see Eq. (10)), which means that the body is nonrotating, as expected. Regarding the functions $\chi$ and $\mu$ appearing in the expression for $d\varepsilon^2/d\sigma$, we have $\chi(\sigma', \sigma) = 0$ and

$$\mu(\sigma', \sigma) = \begin{cases} -\frac{\sigma'^3}{3\sigma^3}, & \sigma' < \sigma \\ -\frac{1}{3}, & \sigma' \ge \sigma \end{cases}. \qquad (D.2)$$

This thus imposes $d\varepsilon^2/d\sigma = 0$ with BC $\varepsilon^2(1) = 0$. We then recover $\varepsilon^2(\sigma) = 0$ (i.e., all isopycnic surfaces are spherical, as expected). Finally, the function $\eta$ needed to compute the enthalpy field becomes

$$\eta(\sigma', \sigma) = \begin{cases} \frac{2\sigma'^3}{3\sigma}, & \sigma' < \sigma \\ \sigma'^2 - \frac{\sigma^2}{3}, & \sigma' \ge \sigma \end{cases}. \qquad (D.3)$$

## Appendix E: Pressure from the enthalpy field

For any isentropic barotrope, the enthalpy and the pressure are related by

$$d\hat{H} = d\hat{P}/\hat{\rho}, \qquad (E.1)$$

where $\hat{P} = P/(G\rho_c^2 R_c^2)$ is the dimensionless pressure. So, the pressure profile is given by

$$\hat{P}(\sigma) = \hat{P}(1) - \int_{\sigma}^{1} d\sigma' \hat{\rho}(\sigma') \frac{d\hat{H}}{d\sigma'}. \qquad (E.2)$$

We can compute the enthalpy gradient from Eq. (4a). We can show that

$$\frac{d\hat{H}}{d\sigma} = -2\pi \left[ 2\sigma\bar{\varepsilon}^2(\sigma) - \bar{\varepsilon}^2 \frac{d\varepsilon^2}{d\sigma^2} \right] \int_{\hat{\rho}(0)}^{\hat{\rho}(1)} d\hat{\rho}(\sigma')\mu(\sigma', \sigma). \qquad (E.3)$$
## Appendix C: The enthalpy profile

Along the polar axis $R = 0$, the Bernoulli equation Eq. (1a) reads

$$H(0, Z) + \Psi(0, Z) = C^{\text{te}}. \qquad (C.1)$$

In order to estimate $\Psi(0, Z)$, we need to go back to Paper II. The gravitational potential of a heterogeneous system made of L homogeneous spheroidal layers is given in that article by Eq. (14). At point $A_j$, at the intersection between the polar axis and the interface between layers $j$ and $j + 1$, it reads

$$\frac{\Psi|_{A_j}}{-2\pi G a_j^2} = \sum_{i=1}^{j} (\rho_i - \rho_{i+1}) \frac{\bar{\varepsilon}_i}{\varepsilon_i^3} \left[ (1 + c_{i,j}) \arcsin\left(\frac{q_{i,j}\varepsilon_i}{\sqrt{1+c_{i,j}}}\right) - q_{i,j}\varepsilon_i\bar{\varepsilon}_j \right] + \sum_{i=1}^{j} (\rho_i - \rho_{i+1}) \frac{\bar{\varepsilon}_i}{\varepsilon_i^3} \left[ (1 + c_{i,j}) \arcsin(\varepsilon_i) (q_{i,j}^2\varepsilon_i^2 + \bar{\varepsilon}_j^2) - \frac{\varepsilon_i\bar{\varepsilon}_j^2}{\bar{\varepsilon}_i^2} \right], \qquad (C.2)$$

where $a_j$ is the equatorial radius of layer $j$, $\varepsilon_j$ is its eccentricity, $\rho_j$ is its mass density, $\bar{\varepsilon}_j = \sqrt{1 - \varepsilon_j^2}$, $q_{i,j} = a_i/a_j$, and $c_{i,j} = q_{i,j}^2\varepsilon_i^2 - \varepsilon_j^2$. Oblate spheroids were assumed in Paper II, but it can be generalized to any spheroid by introducing the function $\mathcal{A}$ defined in Appendix B. In the continuous limit $L \to \infty$, we have

$$\frac{\Psi(R=0, Z)}{G\rho_c R_c^2} = 2\pi \int_0^1 d\hat{\rho}(\omega')\eta(\omega', \omega; \varepsilon). \qquad (C.3)$$

As $H(0, Z)$ is a function of $\omega$ only, we define $\hat{H}(\omega) = H(0, Z)/(G\rho_c R_c^2)$. By Eq. (C.1), we have

$$\hat{H}(\omega) + 2\pi \int_{0}^{1} d\hat{\rho}(\omega')\eta(\omega', \omega; \varepsilon) = \hat{H}(1) + 2\pi \int_{0}^{1} d\hat{\rho}(\omega')\eta(\omega', 1; \varepsilon), \qquad (C.4)$$

where the constant $C^{\text{te}}$ was evaluated at the pole of the whole object (i.e., at $\omega = 1$). We then see that Eq. (C.4) thus yields Eq. (4a).

## Appendix D: The spherical limit

For $\varepsilon^2 \ll 1$, we have

$$\mathcal{A}(\bar{\varepsilon}) \approx 1 + \frac{\varepsilon^2}{6} \qquad (D.1)$$

for both $\bar{\varepsilon} < 1$ and $\bar{\varepsilon} > 1$. So we can expand functions $\kappa$, $\chi$, $\mu$, and $\eta$ at zeroth order to obtain these functions at $\bar{\varepsilon} = 1$ (i.e., in the nonrotating spherical case). We have $\kappa(\omega', \omega) = 0$, which naturally imposes $\hat{\Omega}^2 = 0$ (see Eq. (10)), which means that the body is nonrotating, as expected. Regarding the functions $\chi$ and $\mu$ appearing in the expression for $d\varepsilon^2/d\omega$, we have $\chi(\omega', \omega) = 0$ and

$$\mu(\omega', \omega) = \begin{cases} -\frac{\omega'^3}{3\omega^3}, & \omega' < \omega \\ -\frac{1}{3}, & \omega' \ge \omega \end{cases}. \qquad (D.2)$$

This thus imposes $d\varepsilon^2/d\omega = 0$ with BC $\varepsilon^2(1) = 0$. We then recover $\varepsilon^2(\omega) = 0$ (i.e., all isopycnic surfaces are spherical, as expected). Finally, the function $\eta$ needed to compute the enthalpy field becomes

$$\eta(\omega', \omega) = \begin{cases} \frac{2\omega'^3}{3\omega}, & \omega' < \omega \\ \omega'^2 - \frac{\omega^2}{3}, & \omega' \ge \omega \end{cases}. \qquad (D.3)$$

## Appendix E: Pressure from the enthalpy field

For any isentropic barotrope, the enthalpy and the pressure are related by

$$d\hat{H} = d\hat{P}/\hat{\rho}, \qquad (E.1)$$

where $\hat{P} = P/(G\rho_c^2 R_c^2)$ is the dimensionless pressure. So, the pressure profile is given by

$$\hat{P}(\omega) = \hat{P}(1) - \int_{\omega}^{1} d\omega' \hat{\rho}(\omega') \frac{d\hat{H}}{d\omega'}. \qquad (E.2)$$

We can compute the enthalpy gradient from Eq. (4a). We can show that

$$\frac{d\hat{H}}{d\omega} = -2\pi \left[ 2\omega\bar{\varepsilon}^2(\omega) - \omega^2 \frac{d\varepsilon^2}{d\omega^2} \right] \int_{\hat{\rho}(0)}^{\hat{\rho}(1)} d\hat{\rho}(\omega')\mu(\omega', \omega). \qquad (E.3)$$
Appendix F: Varying the surface axis ratio and the polytropic index: Computing versus precision

<table><caption>**Table F.1.** Data obtained from dimension reduction compared to values obtained for the full 2D problem with the DROP code, obtained for several surface axis ratios $\bar{\epsilon}_s$ and polytropic indices $n$. The computing times spent to achieve convergence are purely indicative, and may vary a little with the computer load.</caption><thead><tr><th></th><th>N</th><th>n</th><th>ε̄<sub>s</sub></th><th>iterations</th><th>comp. time (s)</th><th>M/(ρ<sub>c</sub>R<sub>e</sub><sup>3</sup>)</th><th>Ω<sup>2</sup>/(Gρ<sub>c</sub>)</th><th>|VP/W|</th><th>RMS</th></tr></thead><tbody><tr><td rowspan="12">this work<br>DROP code</td><td rowspan="3">256</td><td rowspan="3">1.5</td><td rowspan="3">0.950</td><td>32</td><td>0.199</td><td>6.492 × 10<sup>-1</sup></td><td>5.314 × 10<sup>-2</sup></td><td>1 × 10<sup>-4</sup></td><td rowspan="3">4 × 10<sup>-5</sup></td></tr><tr><td>32</td><td>37.0</td><td>6.490 × 10<sup>-1</sup></td><td>5.310 × 10<sup>-2</sup></td><td>1 × 10<sup>-5</sup></td></tr><tr><td>33</td><td>0.217</td><td>5.979 × 10<sup>-1</sup></td><td>1.034 × 10<sup>-1</sup></td><td>5 × 10<sup>-4</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">1.5</td><td rowspan="3">0.900</td><td>44</td><td>48.4</td><td>5.973 × 10<sup>-1</sup></td><td>1.032 × 10<sup>-1</sup></td><td>1 × 10<sup>-5</sup></td><td rowspan="3">2 × 10<sup>-4</sup></td></tr><tr><td>35</td><td>0.230</td><td>5.452 × 10<sup>-1</sup></td><td>1.501 × 10<sup>-1</sup></td><td>1 × 10<sup>-3</sup></td></tr><tr><td>42</td><td>46.6</td><td>5.437 × 10<sup>-1</sup></td><td>1.497 × 10<sup>-1</sup></td><td>1 × 10<sup>-5</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">1.5</td><td rowspan="3">0.800</td><td>37</td><td>0.249</td><td>4.906 × 10<sup>-1</sup></td><td>1.921 × 10<sup>-1</sup></td><td>2 × 10<sup>-3</sup></td><td rowspan="3">7 × 10<sup>-4</sup></td></tr><tr><td>42</td><td>46.1</td><td>4.881 × 10<sup>-1</sup></td><td>1.912 × 10<sup>-1</sup></td><td>1 × 10<sup>-5</sup></td></tr><tr><td>39</td><td>0.261</td><td>4.339 × 10<sup>-1</sup></td><td>2.280 × 10<sup>-1</sup></td><td>3 × 10<sup>-3</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">1.5</td><td rowspan="3">0.750</td><td>44</td><td>48.0</td><td>4.302 × 10<sup>-1</sup></td><td>2.266 × 10<sup>-1</sup></td><td>1 × 10<sup>-5</sup></td><td rowspan="3">1 × 10<sup>-3</sup></td></tr><tr><td>43</td><td>0.285</td><td>3.747 × 10<sup>-1</sup></td><td>2.561 × 10<sup>-1</sup></td><td>5 × 10<sup>-3</sup></td></tr><tr><td>56</td><td>59.7</td><td>3.698 × 10<sup>-1</sup></td><td>2.538 × 10<sup>-1</sup></td><td>2 × 10<sup>-5</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">1.5</td><td rowspan="3">0.650</td><td>50</td><td>0.327</td><td>3.126 × 10<sup>-1</sup></td><td>2.736 × 10<sup>-1</sup></td><td>6 × 10<sup>-3</sup></td><td rowspan="3">2 × 10<sup>-3</sup></td></tr><tr><td>56</td><td>58.2</td><td>3.071 × 10<sup>-1</sup></td><td>2.704 × 10<sup>-1</sup></td><td>2 × 10<sup>-5</sup></td></tr><tr><td>50</td><td>0.324</td><td>2.701 × 10<sup>-1</sup></td><td>2.777 × 10<sup>-1</sup></td><td>6 × 10<sup>-3</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">1.5</td><td rowspan="3">0.617</td><td>74</td><td>73.3</td><td>2.648 × 10<sup>-1</sup></td><td>2.740 × 10<sup>-1</sup></td><td>2 × 10<sup>-5</sup></td><td>2 × 10<sup>-3</sup></td></tr><tr><td rowspan="12">this work<br>DROP code</td><td rowspan="3">256</td><td rowspan="3">0.5</td><td rowspan="3">0.950</td><td>16</td><td>0.102</td><td>2.162</td><td>1.189 × 10<sup>-1</sup></td><td>5 × 10<sup>-4</sup></td><td rowspan="3">7 × 10<sup>-4</sup></td></tr><tr><td>23</td><td>28.5</td><td>2.162</td><td>1.189 × 10<sup>-1</sup></td><td>3 × 10<sup>-4</sup></td></tr><tr><td>24</td><td>0.152</td><td>1.197</td><td>8.259 × 10<sup>-2</sup></td><td>1 × 10<sup>-4</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">1.0</td><td rowspan="3">0.950</td><td>51</td><td>56.4</td><td>1.197</td><td>8.253 × 10<sup>-1</sup></td><td>2 × 10<sup>-5</sup></td><td rowspan="3">7 × 10<sup>-5</sup></td></tr><tr><td>40</td><td>0.244</td><td>3.359 × 10<sup>-1</sup></td><td>3.089 × 10<sup>-2</sup></td><td>1 × 10<sup>-4</sup></td></tr><tr><td>44</td><td>47.3</td><td>3.358 × 10<sup>-1</sup></td><td>3.087 × 10<sup>-2</sup></td><td>1 × 10<sup>-5</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">2.5</td><td rowspan="3">0.950</td><td>51</td><td>0.306</td><td>1.610 × 10<sup>-1</sup></td><td>1.589 × 10<sup>-2</sup></td><td>8 × 10<sup>-5</sup></td><td rowspan="3">4 × 10<sup>-5</sup></td></tr><tr><td>50</td><td>53.5</td><td>1.610 × 10<sup>-1</sup></td><td>1.589 × 10<sup>-2</sup></td><td>3 × 10<sup>-5</sup></td></tr><tr><td>70</td><td>0.425</td><td>6.859 × 10<sup>-2</sup></td><td>7.032 × 10<sup>-3</sup></td><td>8 × 10<sup>-5</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">3.0</td><td rowspan="3">0.950</td><td>65</td><td>68.1</td><td>6.853 × 10<sup>-2</sup></td><td>7.023 × 10<sup>-3</sup></td><td>4 × 10<sup>-5</sup></td><td rowspan="3">5 × 10<sup>-5</sup></td></tr><tr><td>88</td><td>0.516</td><td>2.402 × 10<sup>-2</sup></td><td>2.507 × 10<sup>-3</sup></td><td>1 × 10<sup>-4</sup></td></tr><tr><td>88</td><td>91.0</td><td>2.398 × 10<sup>-2</sup></td><td>2.501 × 10<sup>-3</sup></td><td>7 × 10<sup>-5</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">4.0</td><td rowspan="3">0.950</td><td>127</td><td>0.750</td><td>5.861 × 10<sup>-3</sup></td><td>6.163 × 10<sup>-4</sup></td><td>4 × 10<sup>-4</sup></td><td rowspan="3">2 × 10<sup>-4</sup></td></tr><tr><td>128</td><td>130</td><td>5.818 × 10<sup>-3</sup></td><td>6.112 × 10<sup>-4</sup></td><td>2 × 10<sup>-4</sup></td></tr><tr><td>31</td><td>0.004</td><td>6.505 × 10<sup>-1</sup></td><td>5.336 × 10<sup>-2</sup></td><td>5 × 10<sup>-4</sup></td></tr><tr><td rowspan="12">this work<br>DROP code</td><td rowspan="3">32</td><td rowspan="3">1.5</td><td rowspan="3">0.950</td><td>34</td><td>0.067</td><td>6.483 × 10<sup>-1</sup></td><td>5.307 × 10<sup>-2</sup></td><td>8 × 10<sup>-4</sup></td><td rowspan="3">5 × 10<sup>-4</sup></td></tr><tr><td>31</td><td>0.014</td><td>6.496 × 10<sup>-1</sup></td><td>5.320 × 10<sup>-2</sup></td><td>2 × 10<sup>-4</sup></td></tr><tr><td>34</td><td>0.497</td><td>6.489 × 10<sup>-1</sup></td><td>5.308 × 10<sup>-2</sup></td><td>2 × 10<sup>-4</sup></td></tr><tr><td rowspan="3">64</td><td rowspan="3">1.5</td><td rowspan="3">0.950</td><td>31</td><td>0.050</td><td>6.493 × 10<sup>-1</sup></td><td>5.315 × 10<sup>-2</sup></td><td>1 × 10<sup>-4</sup></td><td rowspan="3">1 × 10<sup>-4</sup></td></tr><tr><td>35</td><td>5.10</td><td>6.490 × 10<sup>-1</sup></td><td>5.310 × 10<sup>-2</sup></td><td>5 × 10<sup>-5</sup></td></tr><tr><td>31</td><td>0.783</td><td>6.493 × 10<sup>-1</sup></td><td>5.314 × 10<sup>-2</sup></td><td>1 × 10<sup>-4</sup></td></tr><tr><td rowspan="3">128</td><td rowspan="3">1.5</td><td rowspan="3">0.950</td><td>31</td><td>33</td><td>6.490 × 10<sup>-1</sup></td><td>5.311 × 10<sup>-2</sup></td><td>3 × 10<sup>-6</sup></td><td rowspan="3">4 × 10<sup>-5</sup></td></tr><tr><td>31</td><td>346</td><td>6.492 × 10<sup>-1</sup></td><td>5.314 × 10<sup>-2</sup></td><td>1 × 10<sup>-4</sup></td></tr><tr><td>31</td><td>2.97</td><td>6.492 × 10<sup>-1</sup></td><td>5.314 × 10<sup>-2</sup></td><td>1 × 10<sup>-4</sup></td></tr><tr><td>this work</td><td>1024</td><td>1.5</td><td>0.950</td><td>31</td><td>11.9</td><td>6.492 × 10<sup>-1</sup></td><td>5.314 × 10<sup>-2</sup></td><td>1 × 10<sup>-4</sup></td><td></td></tr><tr><td>this work</td><td>2048</td><td>1.5</td><td>0.950</td><td>31</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>
## Appendix F: Varying the surface axis ratio and the polytropic index: Computing versus precision

<table><caption>**Table F.1.** Data obtained from dimension reduction compared to values obtained for the full 2D problem with the DROP code, obtained for several surface axis ratios $\bar{\epsilon}_s$ and polytropic indices $n$. The computing times spent to achieve convergence are purely indicative, and may vary a little with the computer load.</caption><thead><tr><th></th><th>N</th><th>n</th><th>ε̄<sub>s</sub></th><th>iterations</th><th>comp. time (s)</th><th>M/(ρ<sub>c</sub>R<sub>e</sub><sup>3</sup>)</th><th>Ω<sup>2</sup>/(Gρ<sub>c</sub>)</th><th>|VP/W|</th><th>RMS</th></tr></thead><tbody><tr><td rowspan="12">this work<br>DROP code</td><td rowspan="3">256</td><td rowspan="3">1.5</td><td rowspan="3">0.950</td><td>32</td><td>0.199</td><td>6.492 × 10<sup>-1</sup></td><td>5.314 × 10<sup>-2</sup></td><td>1 × 10<sup>-4</sup></td><td rowspan="3">4 × 10<sup>-5</sup></td></tr><tr><td>32</td><td>37.0</td><td>6.490 × 10<sup>-1</sup></td><td>5.310 × 10<sup>-2</sup></td><td>1 × 10<sup>-5</sup></td></tr><tr><td>33</td><td>0.217</td><td>5.979 × 10<sup>-1</sup></td><td>1.034 × 10<sup>-1</sup></td><td>5 × 10<sup>-4</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">1.5</td><td rowspan="3">0.900</td><td>44</td><td>48.4</td><td>5.973 × 10<sup>-1</sup></td><td>1.032 × 10<sup>-1</sup></td><td>1 × 10<sup>-5</sup></td><td rowspan="3">2 × 10<sup>-4</sup></td></tr><tr><td>35</td><td>0.230</td><td>5.452 × 10<sup>-1</sup></td><td>1.501 × 10<sup>-1</sup></td><td>1 × 10<sup>-3</sup></td></tr><tr><td>42</td><td>46.6</td><td>5.437 × 10<sup>-1</sup></td><td>1.497 × 10<sup>-1</sup></td><td>1 × 10<sup>-5</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">1.5</td><td rowspan="3">0.800</td><td>37</td><td>0.249</td><td>4.906 × 10<sup>-1</sup></td><td>1.921 × 10<sup>-1</sup></td><td>2 × 10<sup>-3</sup></td><td rowspan="3">7 × 10<sup>-4</sup></td></tr><tr><td>42</td><td>46.1</td><td>4.881 × 10<sup>-1</sup></td><td>1.912 × 10<sup>-1</sup></td><td>1 × 10<sup>-5</sup></td></tr><tr><td>39</td><td>0.261</td><td>4.339 × 10<sup>-1</sup></td><td>2.280 × 10<sup>-1</sup></td><td>3 × 10<sup>-3</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">1.5</td><td rowspan="3">0.750</td><td>44</td><td>48.0</td><td>4.302 × 10<sup>-1</sup></td><td>2.266 × 10<sup>-1</sup></td><td>1 × 10<sup>-5</sup></td><td rowspan="3">1 × 10<sup>-3</sup></td></tr><tr><td>43</td><td>0.285</td><td>3.747 × 10<sup>-1</sup></td><td>2.561 × 10<sup>-1</sup></td><td>5 × 10<sup>-3</sup></td></tr><tr><td>56</td><td>59.7</td><td>3.698 × 10<sup>-1</sup></td><td>2.538 × 10<sup>-1</sup></td><td>2 × 10<sup>-5</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">1.5</td><td rowspan="3">0.650</td><td>50</td><td>0.327</td><td>3.126 × 10<sup>-1</sup></td><td>2.736 × 10<sup>-1</sup></td><td>6 × 10<sup>-3</sup></td><td rowspan="3">2 × 10<sup>-3</sup></td></tr><tr><td>56</td><td>58.2</td><td>3.071 × 10<sup>-1</sup></td><td>2.704 × 10<sup>-1</sup></td><td>2 × 10<sup>-5</sup></td></tr><tr><td>50</td><td>0.324</td><td>2.701 × 10<sup>-1</sup></td><td>2.777 × 10<sup>-1</sup></td><td>6 × 10<sup>-3</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">1.5</td><td rowspan="3">0.617</td><td>74</td><td>73.3</td><td>2.648 × 10<sup>-1</sup></td><td>2.740 × 10<sup>-1</sup></td><td>2 × 10<sup>-5</sup></td><td>2 × 10<sup>-3</sup></td></tr><tr><td rowspan="12">this work<br>DROP code</td><td rowspan="3">256</td><td rowspan="3">0.5</td><td rowspan="3">0.950</td><td>16</td><td>0.102</td><td>2.162</td><td>1.189 × 10<sup>-1</sup></td><td>5 × 10<sup>-4</sup></td><td rowspan="3">7 × 10<sup>-4</sup></td></tr><tr><td>23</td><td>28.5</td><td>2.162</td><td>1.189 × 10<sup>-1</sup></td><td>3 × 10<sup>-4</sup></td></tr><tr><td>24</td><td>0.152</td><td>1.197</td><td>8.259 × 10<sup>-2</sup></td><td>1 × 10<sup>-4</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">1.0</td><td rowspan="3">0.950</td><td>51</td><td>56.4</td><td>1.197</td><td>8.253 × 10<sup>-1</sup></td><td>2 × 10<sup>-5</sup></td><td rowspan="3">7 × 10<sup>-5</sup></td></tr><tr><td>40</td><td>0.244</td><td>3.359 × 10<sup>-1</sup></td><td>3.089 × 10<sup>-2</sup></td><td>1 × 10<sup>-4</sup></td></tr><tr><td>44</td><td>47.3</td><td>3.358 × 10<sup>-1</sup></td><td>3.087 × 10<sup>-2</sup></td><td>1 × 10<sup>-5</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">2.0</td><td rowspan="3">0.950</td><td>51</td><td>0.306</td><td>1.610 × 10<sup>-1</sup></td><td>1.589 × 10<sup>-2</sup></td><td>8 × 10<sup>-5</sup></td><td rowspan="3">3 × 10<sup>-5</sup></td></tr><tr><td>50</td><td>53.5</td><td>1.610 × 10<sup>-1</sup></td><td>1.589 × 10<sup>-2</sup></td><td>3 × 10<sup>-5</sup></td></tr><tr><td>70</td><td>0.425</td><td>6.859 × 10<sup>-2</sup></td><td>7.032 × 10<sup>-3</sup></td><td>8 × 10<sup>-5</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">2.5</td><td rowspan="3">0.950</td><td>65</td><td>68.1</td><td>6.853 × 10<sup>-2</sup></td><td>7.023 × 10<sup>-3</sup></td><td>4 × 10<sup>-5</sup></td><td rowspan="3">5 × 10<sup>-5</sup></td></tr><tr><td>88</td><td>0.516</td><td>2.402 × 10<sup>-2</sup></td><td>2.507 × 10<sup>-3</sup></td><td>1 × 10<sup>-4</sup></td></tr><tr><td>88</td><td>91.0</td><td>2.398 × 10<sup>-2</sup></td><td>2.501 × 10<sup>-3</sup></td><td>7 × 10<sup>-5</sup></td></tr><tr><td rowspan="3">256</td><td rowspan="3">3.0</td><td rowspan="3">0.950</td><td>127</td><td>0.750</td><td>5.861 × 10<sup>-3</sup></td><td>6.163 × 10<sup>-4</sup></td><td>4 × 10<sup>-4</sup></td><td rowspan="3">2 × 10<sup>-4</sup></td></tr><tr><td>128</td><td>130</td><td>5.818 × 10<sup>-3</sup></td><td>6.112 × 10<sup>-4</sup></td><td>2 × 10<sup>-4</sup></td></tr><tr><td>127</td><td>0.004</td><td>6.505 × 10<sup>-1</sup></td><td>5.336 × 10<sup>-2</sup></td><td>5 × 10<sup>-4</sup></td></tr><tr><td rowspan="12">this work<br>DROP code</td><td rowspan="3">32</td><td rowspan="3">1.5</td><td rowspan="3">0.950</td><td>34</td><td>0.067</td><td>6.483 × 10<sup>-1</sup></td><td>5.307 × 10<sup>-2</sup></td><td>8 × 10<sup>-4</sup></td><td rowspan="3">5 × 10<sup>-4</sup></td></tr><tr><td>31</td><td>0.014</td><td>6.496 × 10<sup>-1</sup></td><td>5.320 × 10<sup>-2</sup></td><td>2 × 10<sup>-4</sup></td></tr><tr><td>34</td><td>0.497</td><td>6.489 × 10<sup>-1</sup></td><td>5.308 × 10<sup>-2</sup></td><td>2 × 10<sup>-4</sup></td></tr><tr><td rowspan="3">64</td><td rowspan="3">1.5</td><td rowspan="3">0.950</td><td>31</td><td>0.050</td><td>6.493 × 10<sup>-1</sup></td><td>5.315 × 10<sup>-2</sup></td><td>1 × 10<sup>-4</sup></td><td rowspan="3">1 × 10<sup>-4</sup></td></tr><tr><td>31</td><td>5.10</td><td>6.490 × 10<sup>-1</sup></td><td>5.310 × 10<sup>-2</sup></td><td>5 × 10<sup>-5</sup></td></tr><tr><td>35</td><td>0.783</td><td>6.493 × 10<sup>-1</sup></td><td>5.314 × 10<sup>-2</sup></td><td>1 × 10<sup>-4</sup></td></tr><tr><td rowspan="3">128</td><td rowspan="3">1.5</td><td rowspan="3">0.950</td><td>31</td><td>33</td><td>6.490 × 10<sup>-1</sup></td><td>5.311 × 10<sup>-2</sup></td><td>3 × 10<sup>-6</sup></td><td rowspan="3">4 × 10<sup>-5</sup></td></tr><tr><td>31</td><td>346</td><td>6.492 × 10<sup>-1</sup></td><td>5.314 × 10<sup>-2</sup></td><td>1 × 10<sup>-4</sup></td></tr><tr><td>31</td><td>2.97</td><td>6.492 × 10<sup>-1</sup></td><td>5.314 × 10<sup>-2</sup></td><td>1 × 10<sup>-4</sup></td></tr><tr><td>this work</td><td>1024</td><td>1.5</td><td>0.950</td><td>31</td><td>11.9</td><td>6.492 × 10<sup>-1</sup></td><td>5.314 × 10<sup>-2</sup></td><td>1 × 10<sup>-4</sup></td><td></td></tr><tr><td>this work</td><td>2048</td><td>1.5</td><td>0.950</td><td>31</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>