Approaching the structure of rotating bodies from dimension reduction

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.


Introduction
The structure of static polytropic stars, in the classical sense, is traditionnally 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 ⋆ A minimum driver program is available at https://github.com/clstaelen/ssba 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 Ω i of each layer i ∈ [1, L], the parameter set of the spheroids E i , and the massdensities ρ 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 → ∞, we showed that, for global rigid rotation, the eccentricity ε = [1 − b 2 /a 2 ] 1/2 of isopycnic surfaces (with semimajor axis a and semiminor axis b) and the mass density profile ρ 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 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 2006;Espinosa Lara & Rieutord 2013;Rieutord et al. 2016;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 selfconsistent field (SCF) method.In the present case, we have to solve two coupled fixed-point problems, one for the eccentricity profile ε(a) and one for the mass density profile ρ(a), in the range a ∈ [0, R e ].The full 2D structure is reconstructed by unfolding the equatorial solution {ρ(a), ε(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.

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 where Ψ is the gravitational potential, Φ = − dR Ω 2 R is the centrifugal potential, Ω is the local rotation rate, P is the pressure, ρ is mass density of the fluid, and H = dP/ρ is the enthalpy (the self-gravitating flow is isentropic).The link between H and ρ in Eq.( 1c) is usually done via the EOS, namely P(ρ).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 Φ 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 (L in total, ρ 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 where 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 {ρ i ,

Dimension reduction (equatorial projection)
In Paper IV we showed that, in the limit where the number 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 ε(a) and the mass density profile ρ(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 ϖ = a/R e ∈ [0, 1], and ρ = ρ/ρ c is the dimensionless mass density (ρ c is the central mass density).In these conditions, the IDE reads where χ and µ are continuous smooth functions in the whole domain ϖ ∈ [0, 1].These are defined in Paper IV (see also Appendix A, respectively Eq.(A.3) and Eq.(A.4)).As explicitly stated, χ and µ depend on the radius ϖ and on the eccentricity profile ε(ϖ).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 where Ĥ = H/Gρ c R 2 e , and η is reported in Appendix A (see Eq.(A.5)).It follows that Eq.( 1), in the framework of the theory of NSFoE, becomes 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.

Unfolding and global quantities
Once the 1D solution { ρ(ϖ), ε(ϖ)} 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, θ, φ), we actually have where the two coordinates are easily determined by using the parametric formula with θ ∈ [0, π/2] and β ∈ [0, π/2], and where is the axis ratio of the isopycnic E(a, b) considered.The conversion to cylindrical coordinates is straightforward: R = R e ϖ cos(β) and Z = R e ϖ ε(ϖ) sin(β).We can subsequently deduce all global quantities (e.g., mass, volume).In (ϖ, β, φ) coordinates, the volume element writes For most global properties of the system, the integrations over β are trivial.For instance, the mass M and the moment of inertia I are reduced to a single integral over ϖ and we have According to Paper IV, the rotation rate Ω is given by where Ω = Ω/ √ Gρ c and κ 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 Ω2 (ϖ).

Comments
As can be seen in Paper IV, the three functions χ, µ, and η have small amplitude, but all take real, negative values in the range of ϖ 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 ρ ≥ 0 is a firm condition.A similar remark holds for Eq. ( 10).For Eq. (4b), things are less restrictive.Actually, negative values of ε 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 ε = 0 + → 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 ε(ϖ) ≶ 1 instead of ε (see below, and Appendix B).

Cycle and convergence criterion
It is well known that Eqs.( 1a) and (1b) correspond to a fixedpoint problem: ρ = f (ρ) 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 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 ρ(0) (ϖ) and ε(0) (ϖ), and let the profiles evolve until stabilization.This can be accomplished following the iterative scheme, for t ≥ 1: 1. ε(t) is obtained from Eq. ( 7), after integrating Eq. (4b); 2. Ĥ(t) is computed according to Eq. (4a); 3. ρ(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, ε in Eq.( 11b) can be computed using the ρ-profile updated from Eq. (11a), or ρ in Eq. ( 11a) can be computed from the ε-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 take1,2 δ (t) ≲ 10 −14 , where where we define still for t ≥ 1 and ϖ ∈ [0, 1] and use a similar definition for ∆ (t) ε.

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 ϖ = 1 where εs = R p /R e is the axis ratio of the outermost surface.At the center of the body, we have Unfortunately, ε(0) cannot be easily deduced.
Here we decided to recast Eq. (4b) in integral form, i.e., 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: ϖ i ∈ {0, 1 N , . . ., ϖ i = i N , . . ., 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 as the starting guess.Initially, we thus have Ĥ(0) (0) = Ĥc at the center, and ε(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.

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)ρ 1/n , where K is a positive constant and the polytropic index n > 0 is finite.In this case the relationship between Ĥ and ρ is where R 2 e (this relation assumes ρ ≥ 0).The mass density ρ(t) inside the cycle, including the seed, is deduced from the EOS through Eq. (4c).

Example of cycle convergence
The first example is for ε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 ρ(ϖ) and ε 2 (ϖ) during the cycle, as well as the deviations ∆ ρ and ∆ε 2 from one step to the other.Figure 2f gives the convergence parameter δ (t) defined by Eq. ( 12) from  10) and its mean value (orange), and the convergence parameter δ computed from Eq. ( 12) as a function of step t (panel f).
the beginning to convergence.We see that δ (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 δ (t) ∼ 10 −5 , which is on the order of the accuracy of the numerical scheme (i.e., 1/N 2 ∼ 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 κ involved when computing Ω 2 changes sign for most pairs (ϖ ′ , ϖ), 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 reolution.Actually, a test with N = 1024 shows a mean value ⟨ Ω2 ⟩ ≈ 1.3014 × 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 ki- netic energy T , the internal energiy 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 (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.

Example of two-dimensional reconstruction (deprojection)
Once a solution in the form of a pair of profiles ρ(ϖ) and ε(ϖ) 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 εs = 0.75 and n = 1.5, which could correspond to a fast-rotating fully convective star.We note that εs ≈ 0.74 for Achernar (Domiciano de Souza et al. 2014).We can compare the ρ-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 × 10 −3 in absolute (this value is on the order of the virial 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.

Varying the surface axis ratio and the polytropic index: Computing vs. precision
We performed similar comparisons by varying the surface axis ratio εs and the polytropic index n, again with a moderate numerical resolution corresponding to N = 256.The results are gathered in Table F.1 in Appendix F, 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 εs ≈ 0.62.
In the second series the surface eccentricity and the resolution are fixed ( εs = 0.95, N = 256) and n varies from 0.5 to 4. In the last series, the configuration is fixed ( ε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 × 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 ( ρ ≲ 10 −3 at ϖ ≳ 0.7).
The computing time reported in the table is obtained on a standard laptop, without any specification optimization.For n = 1.5 and εs = 0.95, the number of iterations is about 30 and is not sensitive to N. We find ∼ (N/581.295) 1.930 seconds for convergence with the SSB approximation, to be compared to the full 2D problem ∼ (N/77.824) 3.089 seconds.

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 Absolute deviation (log.) Fig. 4. Mass density profile obtained by solving Eq. ( 4) with the SCF algorithm (N = 256) in the static case ( ε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.
density jumps).The convergence criterion and the numerical resolution are (unless stated otherwise) the same as in Sect.3.5.

Static case with index unity
We consider here a static gas with polytropic index n = 1 (see the Introduction).All isobars are spherical and ε = 0 for any ϖ.In this case the solution of the Lane-Emden equation is analytical, 3 namely 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 ε(ϖ) = 0 in the four functions κ, χ, µ, and η, an expansion is required in the limit ε → 0. We find that ε = 0 is in fact a regular singular point, and it follows that κ = χ = 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 quadrature scheme in this case.

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, ε 2 ≪ 1.It happens that there are a few closedform 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). 4For Eq.( 19), the eccentricity profile reads We ran the code for n = 1 and ε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 × 10 −6 .

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 ε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 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 × 10 −4 , while the mean value is about 7 × 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 , ω 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 εs , both for stars and rings (out of range here).The configurations are gathered in the ( j 2 , ω 2 )-plane, where ω and j are the normalized rotation rate Ω and angular momentum J, respectively, defined as (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 ω-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 section 3.6).

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 K domains where the mass density ρ k (ϖ) varies (domain no. 1 is for the innermost domain and k = K is for the outermost one).Then we can write the full mass density profile from the center to the surface as where ϖ k is the position, along the equatorial axis, of the mass density jump between domains k and k + 1, and H is Heaviside's distribution.For convenience we set ρ K+1 = 0 to keep a single sum in Eq. ( 22).The radial derivative is then given by the expression where δ is Dirac's distribution.So, for any continuous function f (ϖ ′ , ϖ) in the interval [0, 1] 2 , we have where ϖ 0 = 0 again for convenience, and the mass density jump α k is defined by 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).

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 amb .Thus, this value corresponds to a cutoff for the mass density at ρamb , and then to a mass density jump at ϖ = 1.In the case of a single domain object (K = 1 and ϖ 1 = 1 here), we have where ρ1 (1) = ρamb .An interesting test case for ambient pressure is the n = 5 polytrope.The analytical solution in the nonrotating case is due to Schuster (Horedt e.g., 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 which verifies ρ(0) = 1 and ρ(1) = ρamb .We computed the solution from the SSB approximation.We note that the method is not appropriate for the case ρ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 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 the two maps, is shown in Fig. 10.We see that the RMS error5 5.
Table 5. Same caption and same conditions as for Table 2, but for the rotating case with ambient pressure.
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 .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 rotator, so the deviation is expected to be on the order of 1− εs ∼ 3×10 −3 in relative.The PREM is a model with K = 10 domains, and as many mass density jumps. 6As ρ(ϖ) is an input, the SSB approximation developed in this work returns to a single fixedpoint 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ε 2 /dϖ and ε 2 (ϖ) in Fig. 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(ϖ) 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 comparted to 0.784 for the PREM).The relative deviations do not exceed ∼ 3 × 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| ∼ 7 × 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 ρ c = 13083.8kg • m −3 , whose value remains close to the PREM value of 1.3088 × 10 4 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 × 10 −3 in relative.Regarding the gravitational 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 1.0826 × 10 −3 1.0712 × 10 −3 1.0771 × 10 −3 −J 4 1.620 × 10 −6 2.96 × 10 −6 2.8233 × 10 −6 † see Chambat et al. (2010) and references therein.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.

Summary
In this article we showed that the 2D structure of a rotating selfgravitating fluid can be determined with good precision from 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 isopynics 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 ∈ [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 endings, 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 isoycnics 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.

Fig. 1 .
Fig. 1.Schematic of the unfolding process.The one-dimensional equation set is first solved for the mass density ρ and eccentricity ε 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).

Fig. 2 .
Fig.2.Radial profiles for ρ (panel a) and ε (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 δ computed from Eq. (12) as a function of step t (panel f).

Fig. 3 .
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 εs = 0.75 and n = 1.5, with 257 nodes per direction (see also Table2).

Fig. 5 .
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).

Fig. 9 .
Fig. 9. Mass density profile obtained by solving Eq.(4b) with the SCF algorithm and N = 256, with configuration εs = 1, n = 5, and ρ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).Table 4. Data associated with Fig. 9 (static pressurized case).exact this work M/ρ c R 3 e 0.60718 0.60720 I/ρ c R 5 e 0.15262 0.15262 Steps for convergence 21 Virial parameter 0 2 × 10 −5

Fig. 11 .
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.

Table 1 .
Global quantities computed for configuration A discussed in Sect.3.3 with the SSB approximation for N = 256 (last column).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.
⋆ value obtained for N = 1024

Table 2 .
Same caption and same conditions as for Table1, but for a fast rotator.

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

Table 6 .
Output dimensionless quantities obtained for the rotating Earth from the SSB approximation, for N = 10240 in total.

Table 7 .
Physical properties of the Earth rescaled; R v = R e ε1/3