Nested spheroidal figures of equilibrium IV. On heterogeneous configurations

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


INTRODUCTION
Unveiling the internal structure of celestial bodies is a longstanding and fundamental challenge in astrophysics.Theories have emerged three centuries ago, with a principal interest in the Earth's interior.In the limit of slow rotation, Clairaut (1743) showed the isopycnic surfaces are spheroids, i.e. ellipsoids of revolution.Using the theory of Maclaurin (1742) for homogeneous spheroids, he obtained a second-order, ordinary differential equation linking the flattening of isopycnics to the mass-density profile.This equation has been more recently extended by Lanzano (1962Lanzano ( , 1974)), the shape of the external surface being expanded over Legendre polynomials  2 up to the -th order.Unfortunately, Clairaut's equation admits essentially no analytical solutions (with some exceptions, see Tisserand 1891;Marchenko 2000).Slow rotators are accessible from the "modified" Lane-Emden equation in the form of series (Chandrasekhar 1933;Kovetz 1968).Besides, Clairaut's equation is limited to small flattenings (i.e. to low rotation rates), while many systems do not belong to the category of slow rotators.This is the case of the giant planets in the Solar System.For Jupiter and Ceres, the flattening parameter  ≈ 0.07, and this is even larger for Saturn (Tricarico 2014;Rambaux et al. 2015).Achernar represents an extreme configuration (Carciofi et al. 2008).New developments remain therefore necessary to model the structure of spinning objects, especially for moderate to fast rotation rates (e.g.Lanzano 1962;Ragazzo 2020).
The determination of the gravitational potential of rotating bodies has always demanded a high analytical effort or substantial computational resources, or both (this exceeds the present context).The spheroidal shape is appealing, as its gravitational potential is known with a closed form (see e.g.Chandrasekhar 1969).Kong et al. (2015) ★ E-mail:clement.staelen@u-bordeaux.fr have investigated the validity of the hypothesis of spheroidal isopycnics by comparing of the "true" shape obtained by numerical means to "perfect" spheroids.They showed that discrepancies are small in amplitude.In fact, this remains true at moderate/large rotation, but unsurprisingly fails close to the mass-shedding limit (Hachisu 1986).Using the gravitational potential of a heterogeneous spheroid, Roberts (1963) used the tensor Virial theorem to derive equations valid for fast rotators, but no self-consistent solutions was produced.
In Huré (2022a,b, hereafter, Paper I and II, respectively), we have investigated the conditions of equilibrium of a piece-wise, heterogeneous system made of L homogeneous layers bounded by spheroidal surfaces.The theory of Nested Spheroidal Figures of Equilibrium (hereafter, NSFoE) assumes that these surfaces stay close to confocality (in the sense of oblate spheroidal coordinates; see Sect. 2) and that each layer can rotate at its own rate (Véronet 1912;Bizyaev et al. 2015).A wide range of configurations is then reachable, from quasispheres to very flat, disk-like objects.It must be pointed out that such solutions remain approximate, although the Virial parameter relative to the gravitational energy is very small (Staelen 2022).This is a consequence of Poincaré-Hamy theorems: a rigidly rotating body with a spheroidal stratification is not an exact figure of equilibrium 1 .In this article, we investigate the solutions in the case where the number L of layers is infinite, which corresponds to a fully heterogeneous body, and for a global rigid rotation.It is therefore a natural continuation of Paper II.Another motiviation of the article is the case of moderate/fast rotators, characterized by a significant oblateness or flatenning (larger than a percent typically).It is therefore interesting to see to what extent classical theories, which are often limited to slow rotation, remain valid or fail.In this purpose, it is necessary to compare any analytical result with numerical solutions.In the present case, this is achieved by using the DROP code which solves the problem for a polytopic equation-of-state (EoS), various flattenings and rotation profiles (Huré & Hersant 2017;Basillais & Huré 2021).
It is obvious that the present approach is not supposed to surpass sophisticated models for stars and planets, which are dynamically and thermodynamically more complex that what the hypothesis made here allow.Stars are widely prone to mixing, transport and circulation.Planets, closer to rigid rotation, have a more simple layered structure (except at the very surface) and isopynics surfaces are believed to be very close to spheroids, as suggested by the inversion of gravitational moments (e.g.Hubbard 2013;Nettelmann et al. 2021).After a brief summary on the theory NSFoE, we show in Sec. 2 how the discrete set of equations can be converted into a differential equation.For rigid, global rotation, this is equivalent to an Integro-Differential Equation (IDE) for the ellipticity of isopycnics.In a first example, we feed this IDE with the numerical solutions obtained from the Self-Consistent-Field (SCF) method (e.g.Hachisu 1986), and show that this approach is not only coherent but quite accurate.In Sec. 3, we study the behavior of the equation in the limit of small flattenings, which happens at slow rotation.In particular, we show that the formalism is fully compatible with classical theories, namely the fundamental second-order differential equation established by Clairaut (1743), the solutions obtained by Chandrasekhar (1933) from the "modified" Lane-Emden equation.We also make a comparison with the equation of Roberts (1963).The question of internal jumps is adressed in Sec. 4, where we derive a modified IDE and test it.In the concluding section, we propose a criterion characterizing the transition from slow to fast rotators, and give a few perspectives.

Equation set for the theory of NSFoE
We adopt the same theoretical background and same notations as in Paper I and Paper II, which can be summarized as follows.We consider L oblate, non-intersecting spheroidal surfaces E  ,  ∈ [[1, L]] with semi-major axis   , semi-minor axis 0 <   ≤   and eccentricity as depicted in Fig. 1.These surfaces define L layers.We note with index  ≥ 2 the layer bounded by E  −1 and E  (we then have   −1 <   and   −1 <   ).Index 1 corresponds to the deepest layer, bounded by surface E 1 only.Furthermore, we assume that each layer  is homogeneous, with mass density   , and rotates rigidly around the -axis at a rate Ω  .A key point in the theory of NFSoE is the possibility of asynchroneous motion of layers, i.e.Ω  −1 ≠ Ω  .In this article, however, we will consider a subclass of configurations characterized by synchroneous rotations.
A fundamental parameter that controls the applicability of the theory is the "confocal parameter"  ,  , defined for each pair (E  , E  ) by This parameter is positive if a surface E  , interior to a surface E  is, in terms of oblate spheroidal coordinates, more oblate than layer .
As quoted in the introduction, only systems with confocal spheroidal surfaces (i.e. ,  = 0 for all pairs) correspond to an exact equilibrium (Poincaré 1888;Hamy 1890).Then, the equilibrium of any layered systems in rigid rotation with non-zero confocal parameters is necessarily approximate.As shown in Paper I and Paper II, equilibria with | ,  | ≲ 0.3 typically are found to be very close to numerical simulations obtained with the DROP-code (Huré & Hersant 2017;Basillais & Huré 2021).Besides, the -parameter is generally found to be slightly negative, meaning that isopycnic surfaces tend to be more spherical with depth in the system, or equivalently, that the ellipticity of isopycnics increases from the center to the surface.However, note that if models of stars and planets mainly agree on such a "standard" stratification, there is no argument or observational proof that definitively rules out a reversal, for some objects.This may depend on physical mechanisms at work and on the formation process.Prolate shapes can be induced by circulations or magnetic fields (e.g.Fujisawa & Eriguchi 2014, and references therein).
The starting point of the present work is (27) of Paper II, which links the properties of all layers together.This is a set of coupled, L − 1 algebraic equations, which read, for  < L, where and is Maclaurin's function defined by (see Paper I), and is the dimensionless rotation rate normalised to the mass density of the uppermost layer ( is the gravitational constant).For  = L (the upper layer), we have 2.2 From a discrete set of layers to a continuum We now seek for equilibrium configurations where  is continuous and derivable from the center to the surface.We first consider configurations without any mass-density jumps (mass-density jumps are considered in Sec. 5).As each layer  has its specific massdensity, and specific set of confocal parameters  ,  , the theory of NSFoE is expected to be capable of such a prolongation, provided these confocal parameters are all "small" enough.When L drastically increases, the extension of layer  in the equatorial plane is   −   −1 ≡ Δ  → 0. In a similar way, at the polar axis, we have   −   −1 ≡ Δ  → 0. Furthermore, the difference in the mass-density between two consecutive layers is In these conditions, (3) can be rewritten as In this form, we see that (9) has the convenient form for the continuous case, as in the limit Δ ρ → 0 the sums over  tend to integrals.To express these integrals, the equatorial radius of layer  is rewritten in the form of the dimensionless, continuous variable Equivalentally, this is the semi-major axis of the isopycnic surface E() ≡ E  , normalised to the equatorial radius of the body.In a similar manner, we associate with the equatorial radius of layer .As long as  1 ≠ 1, we have  1 → 0, otherwise a minimal radius is required.Yet, numerical solutions obtained from the DROP-code (see  show that the deeper the layer, the smaller its flattening, namely ∇ > 0, in agreement with classical theories.Thus, in this work, we will freely take  1 → 0, so that ( ′ , ) ∈ [0, 1] 2 .
In the perspective of a continuous mass-density profile, we must consider situations where the mass density vanishes onto the external surface.In general, the adimensionning used for the discrete theory is not appropriate and must be reconsidered.This is easily corrected.In this purpose, we choose the central mass density  c =  1 as the new reference, instead of  L .The main parameters of layer  are then where  () is the eccentricity of the isopycnic surface E() and ε () is its axis ratio.In a similar way, the confocal parameter is now a continuous variable, namely, from ( 2) With these definitions, the left-hand side (LHS) of ( 9) becomes while the right-hand side (RHS) of ( 9), more complex, can be written in compact form as where the two functions  in and  out are explicitely given in Appendix A; see (A1) and (A2).Formally, these depend on ,  ′ , and , i.e.  ≡ ( ′ , ; ).As  depends on  ′ or , then there are only two variables on input.Despite apparences, these functions behave very well over the integration range.Among interesting properties, we have  in (, ) =  out (, ).This is particularly important and attractive for numerical applications.We also see that  in (, 0) = 0 for  ≠ 0 ( in is not defined for this value).The typical shape of the these functions is visible in Fig. 2b, where we have plotted  in and  out as functions of  ′ for fives values of .For this exemple, we have prescribed a parabolic profile for the eccentricity (see Fig. 2a), as observed in many numerical experiments (see, e.g.configuration A discussed below).It follows from ( 14) and ( 15) that (9) reads, in the continuous limit This equation is the main equation of the present problem.It links the eccentricity of the isopycnic surfaces and their mass density (in fact, its derivative) to the variations of the rotation rate.It enables to reach configurations where both the mass density and rotation rate vary smoothly with the equatorial radius.We can apply the same transformation to (7), and we obtain at the surface We see that the RHS of ( 16) and ( 17) coincide for  = 1, within a factor (d ρ/d)| =1 .So, the two LHS must also coincide for  = 1, which imposes the condition d Ω2 d Thus, the squared rotation has an extremum at the surface.

The case of global rigid rotation: the general Integro-Differential Equation (IDE) for the eccentricity
We see that ( 16) is capable of modeling a wide range of situations, from rigid to differential rotation, and independently, from homogeneous to heterogeneous mass-density profiles.In this work, we focus on rigidly rotating bodies, so we have d Ω = 0.It means that ( 18) is naturally satisfied.Therefore, (16) becomes This equation yields the rotation rate of the body once the configuration is known through ρ() and  ().Alternatively, it can be used to constrain the solutions if the rotation law is prescribed in advance.
We can take the derivative of ( 19) with respect to .In this purpose, we use Leibniz's integral rule, namely, for a given derivable function , d d where  0 is a constant.In the present case, it leads to where we have used the property that  is continuous at  ′ =  (see above).In fact, the partial derivatives can be put in the form and where ,  and  are defined in the Appendix A; see (A3), ( A4) and (A5), respectively.Like , these functions depend on 3 quantities, ,  ′ and , but implicitly only on the 2 space variables  and  ′ .An illustration is given in Fig. 3 for the parabolic eccentricity profile considered previously.An important property is that (, ) = 0 and (, ) = (), which means that the derivatives of the functions are equal at the connection, i.e.

𝜕𝜅 in 𝜕𝜛
This is visible in Fig. 2b.These functions have also a relatively small amplitude, which, again, is very practial for any numerical treatment.From ( 22) and ( 23), ( 21) becomes 2 which then links directly the eccentricity of the isopycnic surfaces to their mass density.
A consequence of (25) comes from the case  = 0 (i.e. the center of the body).Indeed, at this point, we have As ( ′ ) < 0 ∀ ′ (this is seen from its definition; see (A5) in the Appendix) and d ρ < 0 from stability consideration, (26) yields whatever the mass density profile.

A note on the condition of immersion
As early quoted in this section, an important hypothesis of the theory of NSFoE is the non-intersection of the interfaces between layers.In the continuous limit, this means that the isopycnic surfaces must not cross each other.Two ellipses intersect if and only if the one with the largest major axis has also the smallest minor axis.So, if the polar radius , i.e. the minor axis, is given by It is easily shown from the definitions of  and ε that So, (28) can be written as which then imposes an upper limit for the eccentricity gradient.
Note that requiring (d . the body would be infinitely flat.

The particular case of homogeneity: Maclaurin formula recovered
A first check of the ( 25) is performed by considering the Maclaurin spheroid.In this case d ρ = 0, and the mass density profile is where H is Heaviside's step function.In the sense of distributions, the derivative of the mass density is where  is Dirac distribution.We can now use (25) to obtain , which is the only unknown of the problem.However, for a body where the mass density is a constant, the notion of isopycnic surfaces appears as a non-sense.Yet, the Poincaré-Wavre theorem implies that for a body where the rotation rate is constant on cylinders, isopycnic and isobaric surfaces must coincide2 .Rigid rotation has a rate which is obviously constant on cylinders, so we can consider  as the eccentricity of isobaric surfaces.It can be shown that ( 16) and ( 25) become and respectively, where we used the properties of the Dirac distribution.From (A5), we see that  never vanishes, so the only solution to ( 34) is which means that isobaric surfaces are similar spheroids.So, by expliciting  out , (33) becomes where  s ≡  (1) and εs ≡ ε (1) are values at the surface.As expected, we fully recover the results from Maclaurin's theory.

Checking the IDE from a numerical reference
Unfortunately, without any prior knowledge on the mass-density profile or the eccentricity, (25) can not solely be used to determine any internal structure.However, we can test the reliability of the above approach.In this purpose, we find more practical to rewrite (25) in the form The test then consists in computing both sides of this expression by solutions obtained numerically from a Self-Consistent Field (SCF)method.We use the DROP code (Huré & Hersant 2017;Basillais & Huré 2021) as the numerical reference.This code, which has been extensively used, solves the full structure of rotating, self-gravitating fluids for a wide range of flattenings, equation of states and rotation profiles3 .A fundamental ingredient is the closure relationship between pressure  and mass density .In the paper throughtout, we use a polytropic EoS, namely where  and  (the polytropic index) are positive constants.Once the SCF-cycle has converged, the mass-density (, ) is known.Then, we have to determine isopycnic surfaces, denoted S  , from center to surface.Clearly, isopycnics are not exact spheroids (in general, these are sligthly depressed in the middle), but any equilibrium surface S  crosses the polar axis and equatorial axis respectively at points A  (0,   ) and B  (  , 0) (see Fig. 1).From these two points, we can calculate a "pseudo-eccentricity", basically from (1).This pseudo-eccentricity is then of the form  ().We can then use this output, together with the mass-density () along the equatorial plane to compute , ,  and the derivative of the pseudo-eccentricity d 2 /d, and then check (37).We will also compare our results to Clairaut's integral equation, i.e. ( 43), which will be discussed in Sec.
3, and to Roberts' result,which is given in our notations in Appendix B ; see (B1).
There are four main sources of errors in this kind of numerical test.First, DROP releases numerical solutions whose accuracy depends on the resolution.Second, the determination of equilibrium surface S  (and then, points A  and B  ) is also not perfect.Next, the integrals in the RHS of (37) are also sensitive to the quadrature scheme, as well as the scheme for the derivative of the pseudo-eccentricities (here, we use 2nd-order schemes).Obviously, we do not expect (37) to be exactly satisfied.In turn, if both sides of this equation are very close for a broad variety of configuration, then it proves the reliability of the IDE.

An example
As a first illustration, we consider a rotating polytrope with surface axis-ratio εs = 0.75 and polytropic index  = 1.5, hereafter Configuration A. It corresponds to a fast rotator (for comparison, Achernar has an axis ratio in surface around 0.74; see Domiciano de Souza et al. 2014), which is also one of the structures given in the tables of Hachisu (1986).The mass-density, pseudo-eccentricity and the deviation of the outermost surface to an exact spheroid are displayed in Fig. 4a to c (left panels).The RHS and LHS of (37) are plotted versus  in 4 d.We see that the absolute deviation between these two estimates (panel e) is much less than 1% for most radii, and even mess than 0.1% in the outerpart of the body.This agreement is already remarkable as the conforcal parameters (center and surface values) are marginally acceptable (i.e.(, 1) ∈ [−0.4375, 0]).The figure also shows that the approximation is also valid for d 2 /d, as the discrepancy with DROP is also of order ∼ 10 −3 in this case.
From panel f, we see that the rotation rate Ω deduced from ( 19) is not strictly a constant, as would be expected.But, we see that it varies weakly and compares greatly with the rotation rate yielded by DROP, with an error below a percent.We see that Roberts' equation compare greatly with the numerical reference, except at short (where a divergence is seen) and large radii (with an error of ∼ 3 %).
We have calculated the main global properties of the polytrope, namely the mass , the volume  and the angular momentum , the gravitational, kinetic and internal energies, ,  and  respec-  tively (see the Appendix C) and compared with the tables of Hachisu (1986).The results are reported in Tab. 1.We see that the values obtained are slightly overestimated with the present approximation.This is due to the boundary of the fluid, which is below the corresponding spheroidal surface, as seen from Fig. 4c.Thus, the volume of the fluid, and all volume integrals following, are clearly greater than the outputs of the numerical reference.Furthermore, the value of the Virial parameter, i.e. |VP/ | ≈ 8 • 10 −4 ≪ 1, also validates the approximation in this case.

On critical rotations
We can go further in the comparison by looking at an extreme configuration, i.e. a configuration near the so-called "critical-rotations" (Hachisu 1986), where matter at the surface is bearly bounded to the system.Such objects deviate largely from spheroids and we expect the approximation to fail at this point.We first consider configuration B, with a "soft" EoS ( = 3).The configuration and its global properties are reported in Tab. 2 and the results are plotted in Fig. 5. Surprisingly, the agreement between the spheroidal approximation and the numerical reference is very good; see Fig. 5d and Tab. 2. For  > 0.2, we see that the discrepancy is ≲ 10 −3 in relative.For shorter radii, the gap is wider, due to the numerical precision of the derivatives, as the values themselves are "small" ; so any discrepancy is amplified.The approximation seems to stay valid at the surface, even though the deviation from a spheroid is large (see panel c).This can be explained by the mass density curve, namely panel a.Indeed, we see that, for  > 0.4, we have () ≪  c , so the contribution of this part to the gravitational potential (and thus, to the rotation rate and (25)) is negligible.So, as long as the isopycnics for  < 0.4 are close enough to spheroids, the approximation is still valid.We also have plotted in panel d of Fig. 5 the upper limit of the immersion criterion, i.e. (30).Interestingly, the squared eccentricity gradient seems to tend to this limit for  = 1, i.e. at the surface.This would imply that at the critical rotation, we have (d/d)| = e → 0, where we used the physical radii, namely the matter at the pole is crushed.
Another example of critical rotation is displayed in configuration C, where the EoS is "hard" ( = 0.5).The configuration and its global properties are reported in Tab. 3 and the results are plotted in Fig. 6.Here, the agreement between the spheroidal approximation reported here and the numerical reference is not good at all, with a relative error of at least 10% on d 2 /d and the global properties.Only the averaged rotation rate is correct, but we see from Fig. 6f that the rate itself is not a constant anymore (with an amplitude of, again, ∼ 10% of the mean value).This disagreement is explained by the large deviation of the external surface to a spheroid, which is not cancelled by the mass density profile, i.e. ρ() ≪ 1 only very close  to the surface ( = 1).So, the deviation from a spheroid has here a real impact, as the gravitational potential arising from this mass distribution is significatively different from the one produced by a spheroidally stratified object.However, we observe once again that the immersion criterion joins with the d 2 /d-curve computed from DROP at  = 1, reinforcing our conclusion of the previous example.

The IDE at first order
The case of slowly rotating structures is of great importance in the context of planetary and stellar interiors (e.g.Chandrasekhar & Roberts 1963;Zharkov & Trubitsyn 1970).Such situations suppose that the deviation to sphericity is small, i.e.  2 () ≪ 1.While the Earth or the Sun can probably be considered as slow rotators, this does not seem to be the case of Jupiter and Saturn.The functions defined by (A3), ( A4) and (A5) can then be expanded at first order in  2 .So, we obtain and respectively.Thus, at first order in  2 , (25) becomes Note that d 2 /d is already first order in  2 , so the first order terms arising from  and  can be neglected.

Clairaut's equation recovered
Except in some particular cases, the mass density vanishes continuously at the surface.By integrating (42) by parts, we obtain is the classicaly called the mean density (e.g.Tisserand 1891; Ragazzo 2020), evaluated from the center to the running radius.
In this form, ( 43) is suitable to eliminate the integral by differentiation.So, we derivate a second time with respect to the physical radius This result clearly recalls the fundamental equation derived by Clairaut (1743), namely where  = 1 − √ 1 −  2 is the flattening of the isopycnic surface,  is its polar radius and Let us show that ( 45) and ( 46) are fully compatible.At first order in  2 , we have 2  ≈  2 and  ≈ (1 −  2 /2), so Now, as the derivatives and the function  itself are already of first order in  2 , only the "zeroth" order in ⟨ ρ⟩ is needed.At this order, we have  ≈  and thus ⟨ ρ⟩ ≈ ρm .Hence, we conclude that ( 25) is equivalent to Clairaut's differential equation in the limit of small flattenings, at first order in  2 .Note that some authors (e.g.Ragazzo 2020) use the mean radius ( 2 ) 1/3 instead of  or .We can show by the same reasoning that the equations would still agree at first order.

An example. Comparison with Chandrasekhar's pertubative approach
To illustrate the compatibility between Clairaut's equation and ( 25), let us consider the numerical solution computed from DROP for a self-gravitating polytrope with εs = 0.99 and  = 1, hereafter configuration D ; see Tab. 4 for the details of the configuration and the associated global quantities.We have  2 s = 0.0199, which is expected to be "small enough" for the expansions made in the previous paragraph to be valid.We can therefore check our expansions as well as Clairaut's equation.The results are presented in Fig. 7 (same panels as for configuration A).We notice that the -profile is close to a quadratic.We see, again, the excellent agreement between the present approach and Clairaut's equation.Also, we see that the global quantities obtained with the IDE are close to the one obtained with DROP, with between four to six digits shared on the values.As quoted in the introduction, this is not a surprise, as Clairaut (1743) showed that for small deviations from the sphere, i.e. small flattenings, the isopycnic surfaces are ellipses in any meridian plane.
Moreover, slowly rotating polytropes have been studied by many authors, in particular by Chandrasekhar (1933).His approach is based on the Lane-Emden equation, supplemented by a small amplitude, rotational field.The equilibrium is solved in the form of series.Configurations with  = 1 (like configuration B) are interesting because the results arising from this theory are purely analytical and offer a interesting opportunity for comparisons.As the dimensionless rotation rate Ω2 is an input in Chandrasekhar's work (while the axis ratio ε is an output), the comparison is performed by injecting the rotation rate provided by DROP into Chandrasekhar's equations.The results are reported in Tab. 4 (column 2).We see that the comparison is satisfactory, the agreement being much better than 1%.Furthermore, the Virial quantities yielded by the spheroidal approximation are in excellent agreement with the numerical reference.

INTRODUCTION OF MASS-DENSITY JUMPS: THE MODIFIED IDE
Mass-density jumps are usually associated with a sudden change in the equation of state or in the mechanism transporting matter or energy.It is therefore interesting to render the present method as flexible as possible, and to account for such discontinuities.As often, we consider jumps as zero-thickness transitions, while, in real systems, these have always have certain spatial extension.Inspired by Sec.2.5, we can easily introduce mass-density jumps in the present  formalism by decomposing the mass-density profile as where K is the number of heterogeneous domains and ρ () =   ()/  is the mass density inside domain number  (we still normalise mass-densities to the central value   ).As for the discrete case, we have set ρK+1 () = 0 to keep a single sum, which means that the outer space is the very last domain, with index K + 1 and null mass density.There are therefore K jumps, located at Note that (49) allows for configurations with a surface discontinuity, i.e. at  K = 1.The derivative of this profile writes We can thus make use of the properties of the Heaviside and Dirac distributions to generalise (25).

Piece-wise rotation and discontinuity in the ellipticity
Let us consider that each domain  ∈ [[1, K]] rotates rigidly at its own rate Ω .So, for a given domain  0 , we have  ∈ ]  0 −1 ,   0 [, and ( 16) becomes4 where   = ρ (  )/ ρ+1 (  ) is the mass density jump at each interface .A major question concerns the behavior of this equation when applied to two adjacent domains.To answer this point, we write (51) at  − =   0 − Δ (inside layer  0 ) and at  + =   0 + Δ (inside layer  0 + 1), with Δ > 0. In the limit where Δ ≪ 1, the difference in the rotation rates between  − and  + satisfies at first order in Δ.If asynchroneous motion is possible, then the RHS of this expression must remain finite when Δ → 0. We see from (52) that this is possible only if the eccentricity undergoes a discontinuity at   0 , namely d 2 d where   0 is the eccentricity profile in the domain  0 .Note that (52) can not be used to quantify this jump, as we assumed a continuous eccentricity to arrive at this point.Indeed, if these jumps are considered from the beginning, they would cause discontinuities in the -functions, ,  and , which makes the calculations far more complex.
This "eccentricity jump" only states that the interfaces between layers are not isopycnic surfaces.The isopycnic in the inner layer (the "core") intersect the interface and is prolongated by another isopycnic in the outer layer (the "envelope") whose eccentricity has no reason to be the same.This statement has two interesting consequences: i) the "eccentricity jump" occurs not at a single value of  but on a whole range close to any interface ; ii) the potential of an incomplete Maclaurin spheroid being unknown analytically, the continuous version of the NSFoE cannot describe systems with rotational discontinuities.

Global, rigid rotation
By requiring Ω = Ω, ∀ ∈ [[1, K]], both sides in (52) are null in the limit Δ → 0, meaning no eccentricity jump occurs for systems in global rotation, so that the interfaces between layers are isopycnic surfaces.As the RHS of ( 51) is constant, we can, as in the singlelayer case, take its derivative with respect to  inside layer  0 .We find This expression is the IDE modified by the presence of jumps.Note that it can be recast in the form of (37).As for the single-layer case, (54) can not be solved alone as we have a single integro-differential equation for K + 1 unknown functions, namely the mass-density profiles ρ () and the eccentricity  ().A solution requires K equations of state and K Bernoulli's equations.

An example
Once again, we check the self-consistency of (54) by comparison with a numerical solution from DROP; see Subsec.2.6.We see that (54) can be written in the form of (37), i.e. we can obtain an equation   of the form d 2 /d = (, , ).So, as before, we use DROP outputs to compute both sides of (54) and we then compare the results.Configuration A' is a rotating body with surface axis ratio of 0.75, a core with polytropic index  1 = 1.5 and semi-polar axis  1 ε ( 1 ) = 0.35 and an envelope with polytropic index  2 = 3.This system could correspond to a highly object with a large convective core (whose mean radius is ∼ 40% of the star's radius) and a big radiative envelope; it may thus be considered as a very simple model for a fast-rotating high-mass star ( ≳ 1.2 M ⊙ ); see e.g.Maeder (2009).The global quantities are given in Tab 5 and the results are plotted in Fig. 8. Again, the agreement between the spheroidal approximation reported here and the numerical reference is remarkable, within a few tenths of a percent (except for the volume).The relative Virial parameter is also really good (10 −3 ≪ 1), which validates more the approach.We see that both squared excentricity gradients compare really well to each other (see panel d), the discrepancy being around ∼ 10 −3 in most of the object and around ∼ 10 −2 in the neighboring of the mass density jump, which is due to the numerical resolution in this region.Indeed, for each cylindrical radius, the interface is described by two or three points, which may not be enough to reach a good accuracy on the dynamics of the eccentricity in this region.This peak is also seen in the Ω2 curve (panel f), where the gap to the value yielded by DROP is also about a few tenth of a percent.

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

Open questions and perspectives
1. Rotational discontinuities.When rotational discontinuities are present, an eccentricity jump is mandatory, meaning the interfaces between layers were not isopycnic surfaces.The approximation of spheroidal isopycnic then fails in this case.However, if the rotational discontinuities (or equivalently, the eccentricity jumps) are small enough, it should be possible to derive an IDE for this as the range where the jump occurs becomes negligible.This point would merit an additional work.
2. From slow to fast rotator: a criterion.In the limit of small flattenings, our approach compares really well with the one developped by Chandrasekhar (1933) and we were able to recover Clairaut's equation at first order in  2 .This adresses the question of the limit between slow rotators (well described by Clairaut's theory) and fast rotators, which can be roughly answered as follows.Let us develop  at second order in  2 (for convienience, we use  instead of  or  as it is a function of a single variable).From (A5), we directly obtain Now, let  be the ratio of the fourth order term to the second order term.We have So, roughly, the error in the quantity d 2 /d made by using Clairaut's equation, i.e. (45), is of the order of .The corresponding axis ratio at the surface is then εs We see that for configuration D (Fig. 7d), which has an axis ratio of 0.99 at the surface, the maximum error is of order 10 −2 (we do not take into account the part  < 0.2 which is dominated by the errors of the finite-difference scheme).This then corresponds to the criterion (57).
3. Can we expand the IDE at higher orders ?As shown, expanding the IDE at first-order in  2 leads to Clairaut's equation.It would then be interesting to derive a second-order Clairaut equation, basically by exanding the IDE a second-order expansion in  2 .This would be another approach to the expansion of Clairaut's equation than Lanzano (1962Lanzano ( , 1974) ) who has performed a multipolar expansion of the shape of the object.However, preliminary calculations indicate that the problem might not be any easier than the equation set reported here.This point is still under investigation.

Do exact solutions to the IDE exist ?
As it is well known, analytical solutions are always powerful tools for making models and diagnosis tools, regarding observations.The existence of analytical solutions to the IDE in the form () would be very interesting, and it already represents an exciting perspective.Clairaut's equation is known to have a few analytical solution (e.g.Tisserand 1891;Marchenko 2000).Given the complexity of the IDE, we expect any analytical solution to be only approximate.Solutions via a series expansion or linearisation for example would be interesting to seek for.
5. Towards 2D-structures ?As quoted, (25) is not sufficient in itself to derive models for interiors of rotating bodies; it is the case of Clairaut's equation as well.The IDE has to be combined with an EOS and to Bernoulli's equation.However, the IDE enables to reduce the number of dimensions of the problem, from two to one, through the relationship ().The computation of the gravitational potential is skipped in this process (in fact, it is already incorporated in the IDE).This is very attractive, in particular in terms of computing time if a large number of structures have to be computed (see below).We are currently preparing an article dealing with the structure of spheroidal stars and planets from a SCF-method (Hachisu 1986) through this dimension reduction.
6. Inverse problems.Planets like Jupiter and Saturn do probably not belong to the category of slow rotators.The IDE could therefore be of great help in generating fast internal 2D-structures (with appropriate EOS), under the conditions of the hypothesis of the NSFoE.Next, it would be easy to compute the gravitational moments and to isolate solutions that match the values "measured" by space probes.Yet, as pictured by e.g.Miguel & Vazan (2023), high-order gravitational moments mostly describe the outer layers of the object, which is the most poorly described zone by the theory reported here; see also Basillais & Huré (2023) (Paper III) and references therein.As such, we expect only the first two moments to be accurate enough.Furthermore, as quoted by Nettelmann et al. (2021), the Concentric Maclaurin Spheroid (CMS) method by (Hubbard 2013) has high computational needs, meaning that a scan of a given parameter space is tedious.With a very fast algorithm, it could be possible to identify places in the parameter space compatible with the measured  2 , which could be further studied with more sophisticated algorithms (e.g. the CMS-method).Obviously, in the case of gaseous planets, the presence of complex winds at the very surface is not strictily compatible with the NSFoE (the 3D-structure of a gaseous planet with zonal winds has been studied by Kong et al. 2016).This is worst in stars where meridional circulations are present (see e.g.Zahn 1992). .

(C2)
The Jacobian matrix  of the transformation from cartesian coordinates to an isopycnic coordinate system (, , ), where  is the spherical azimutal angle, reads All the integrals of (C1) can now be computed numerically (via a trapezoidal rule for instance).
For the rotation rate, as it is not exactly a constant due to the spheroidal approximation, we can obtain an mean value by integrating over the mass, namely

Figure 1 .
Figure 1.Typical configuration for a heterogeneous body with finite number of homogeneous layers (L in total) bounded by spheroidal surfaces E  .

Figure 4 .
Figure 4. Results for configuration A ( εs = 0.75,  = 1.5); see Tab. 1 for global quantities.Left-hand side panels: outputs from the DROP code, i.e.(a) eccentricity of the isopycnic layers (the dashed thin black line is the quadratic profile used for Figs. 2 and 3); (b) radial mass density; (c) the deviation of the external surface from a spheroid.Right-hand side panels: Comparison between this work and the output from DROP, i.e.(d) d 2 /d as a function of ; (e) decimal logarithm of the gap between the analytical methods and the numerical reference for d 2 /d; (f) decimal logarithm of the gap between DROP and this work for the rotation rate.

Figure 8 .
Figure8.Same legend as for Fig.4, but for configuration A', which is a two-domain body (i.e.K = 2) with a mass-density jump at  =  1 ≈ 0.36 (marked with a vertical red dshed line; see Tab. 5).

Table 1 .
Hachisu (1986)A and corresponding global quantities.Results from the tables ofHachisu (1986)are reported in the first column.

Table 2 .
Same legend as Tab. 1, but for configuration B.

Table 3 .
Same legend as Tab. 1, but for configuration C.

Table 5 .
Same legend as for Tab. 2 but for configuration A'.