An improved algorithm of second order to construct consistent theories of equilibrium figures of rotating celestial bodies

One of the main problems in celestial mechanics is the study of the figure adopted by large deformable bodies in slow rotation around an axis with a constant angular velocity ~ω when they reach their equilibrium configuration. This figure corresponds to the lowest equipotential surface containing the entire mass and, in order to determine it, to calculate its potential at an arbitrary point is required. Classical methods address this problem generally by performing a series development of the inverse-distance by using Clairaut’s coordinates. These methods show convergence problems, already in first order in ω, so that to avoid them they must assume no demonstrated hypotheses. The authors of this paper warned and proved this fact in a previous work, for which they used two methods: 1. Taking into account the asymptotic properties of numeric quadrature formulas. 2. By a process similar to that used by Laplace to develop the inverse of the distance between two planets. Thus, the authors demonstrated that, although the intermediate formulas obtained by the classical methods are wrong, the self-gravitational potential, up to first order in ω, obtained with them, were coincident with the obtained by other methods. Preprint submitted to Journal of Computational and Applied Mathematics July 2, 2018 Manuscript Click here to view linked References

1. Taking into account the asymptotic properties of numeric quadrature formulas. 2. By a process similar to that used by Laplace to develop the inverse of the distance between two planets.
Thus, the authors demonstrated that, although the intermediate formulas obtained by the classical methods are wrong, the self-gravitational potential, up to first order in ω 2 , obtained with them, were coincident with the obtained by other methods.

Introduction
In the following paper the study of the equilibrium configurations of a celestial body formed by a fluid mass endowed with a constant angular velocity ω around a fixed axis is discussed. That axis, which passes through the celestial body centre of masses O is noted as OZ.
As it is well known [4], [7], [8] the stable equilibrium configuration coincides with the minimum potential, which is equivalent to a hydrostatic equilibrium state and this implies a rigid rotation. Consequently, the equipotential, isobaric and isopycnic surfaces coincide.
Let OX, OY be two perpendicular axes connected to the body in rotation and forming a direct trihedral OXY Z.
Let Ψ, P and ρ be the total potential at one point, the pressure and the density at one point in the body, respectively. Then, the problem can be modeled by the Poisson equation and the hydrostatic equilibrium equation given by: The integration of these equations requires a state equation P = P (ρ).
To solve the problem, either direct integration methods or the Clairaut method [6] can be used. The last method uses either the potential or other biunivocally variable determined by it as an independent variable.
In this work as Clairaut coordinates of one point of the body will be used those determined by (a, θ, λ). Variable a is defined as the radius of the O centered sphere containing the same mass as the equipotential surface containing the point and the variables θ and λ are those corresponding to the spherical coordinates of the point.
The theories of deformable equilibrium figures in uniform rotation developed by Finlay [4] and Kopal [7], [8] were based on the use of Clairaut coordinates. However, these theories, as demonstrated in López [10], require for their validity the unproven Laplace desideratum [13].

José A. López Ortí, Manuel Forner Gumbau, Miguel Barreda Rochera
To solve this problem, two algorithms were developed in [10] to approximate the problem solution to the first order in ω 2 . There it was proved that several of the expressions obtained in the works of Finlay [4] and Kopal [7], [8] are not true yet in the first order in ω 2 . In the same work it was also shown that, in spite of this detected incorrectness, the value of the obtained potential in the mentioned works is correct until first order in ω 2 . López [11] extends this technique to the case of binary closed systems obtaining correct results until first order in ω 2 .
In this work we will try to extend up to second order in ω 2 the first algorithm previously mentioned [10]. This obviously entails much more complex and extensive developments than those of first order theory.
The present section offers the general problem and the research objectives. In section 2 the classical theory corresponding to the equilibrium figures of celestial bodies in uniform rotation is presented. In Section 3, the numerical quadrature algorithm previously developed by the authors in first order [10] are extended to second order in ω 2 .
In section 4 a numerical example is presented. Although it is not directly related to the main object of this work, it can be of use to show the efficiency of the Clairaut method to deal with this type of problems.
Section 5, finally, explains the main results and conclusions of this work.

Classical Theory
Let M be a deformable isolated mass with uniform rotation around its mass center, endowed with angular velocity ω and whose mass distribution is given by ρ(x, y, z). Let assume that the mass distribution function is differentiable inside of the mass.
Let OXY Z be the coordinate system associated with the body, defined as follows: • O is the mass center of the rotating body.
• OZ is the axis determined by the straight line passing through O and parallel to the angular velocity vector ( ω) of the rotating body.
• OX is a line passing through O and contained in the OXY plane so that the trihedron OXY Z is direct.
Then, the potential at an inner point whose radius vector has r = (x, y, z) as components will be determined by the following expression: Where Ω is the self-gravitational potential, V c is the centrifugal potential, G is the universal gravitational constant, dm ′ is the mass element of an arbitrary inner coordinate point (x ′ , y ′ , z ′ ) and ∆ is the distance between the coordinate points (x, y, z) and (x ′ , y ′ , z ′ ).
An improved algorithm of second order to construct consistent theories...
An inner point of the body in coordinate rotation (x, y, z) is expressed in spherical coordinates as follows: x = r cos θ cos λ, y = r cos θ sin λ, z = r sin θ.
The inverse of the distance is developed as follows where γ is the angle between the radii vectors of r = (r, θ, λ) and r ′ = (r ′ , θ ′ , λ ′ ) of the points (x, y, z) and (x ′ , y ′ , z ′ ), respectively, expressed in spherical coordinates and P n are the Legendre polynomials. Then, the self-gravitational potential can be expressed as where In the Clairaut coordinate system (a, θ, λ) a parameter is constant on each equipotential surface. In this work we have chosen the parameter a so that it is the radius of the sphere centered in O such that it has the same mass than the equipotential surface ( Figure 1). The Clairaut coordinate system is related to the spherical coordinate system by the relation r = r(a, θ, λ). In consequence, x = r(a, θ, λ) cos θ cos λ y = r(a, θ, λ) cos θ sin λ z = r(a, θ, λ) sin θ. In the Clairaut coordinate system the radius vector r of an equipotential surface is developed ( [4], [6], [7], [8]) as follows: where f n,m (a) are the functions of amplitude and Y n,m (θ, λ) the spherical functions in real form [12]. For symmetry reasons, the radius vector r can be developed as follows: where P n are the Legendre polynomials.
When ω is small (in a slow rotation case), the amplitudes f 2m (a) are small quantities with respect to the unit.
In order to calculate integrals (6) and (7), Finlay [4] and Kopal [8] assume that and where a 1 is the first root of the equation ρ(a) = 0.
To evaluate these integrals up to second order in amplitudes, the following approximations are used: To develop as a linear combination of Legendre polynomials the powers and cross products of Legendre polynomials that are obtained from the different powers of Σ in (14) and (15) the Adams-Newman formula [1] is used On the other hand, taking into account the spherical harmonics addition theorem in real form from its orthonormality and from the fact that where δ n,s is Kronecker's delta. Consequently, from (5)-(19) and after various algebraic manipulations the self-gravitational potential can be expressed as follows: Approaching up to first order in the amplitudes we have In addition, the centrifugal potential V c approximated to the first order in amplitudes is resulting that total potential (2) depends only on parameter a Ψ(a) = Ψ 0 (a), Ψ n (a) = 0 ; n = 1, 2, . . .
Then, from (2), (20) and (23) it is obtained It should be noted that in equations (25), (26) and (27) only f 2 (a) is of the first order with respect to the amplitudes.
Approaching now up to second order in the amplitudes and taking into account (28) and for n = 2, 4 36 35 An improved algorithm of second order to construct consistent theories...
On the one hand, it should be noted that the assumption made by Finlay [4] and Kopal [8] about the fact that (6) and (7) are equivalent to (11), (12) and (13) respectively, is based on The Laplace desideratum which, as indicated in [13], is not a proven fact but a conjecture.
On the other hand, in [10], without making use of the Laplace desideratum, the correct development up to first order in amplitudes of U n (6) and V n (7) is obtained. It is worth noting that, in first order in amplitudes (6) and (7) are not equivalent to (11), (12) and (13), respectively. However, the self-gravitational potential (Ω), developed up to first order in amplitudes by Kopal [8] and by López [10] coincides. Then, although developments up to first order in amplitudes, of the external (U n ) and inner potentials V n carried out by Kopal [8] are incorrect, [10] demonstrates that the development, up to first order in amplitudes, of the self-gravitational potential (Ω) obtained by Kopal [8] and by López [10], is identical.
Kopal, in [8] pp 28-36, explicitly obtains the developments of f 2 (a), f 4 (a) and f 6 (a). One of the consequences derived from these developments is the proof that f 2 (a) is of first order in ω 2 .
Consequently, in accordance with the above considerations, the results obtained by Kopal [8] involving the self-gravitational potential (Ω) and, particularly, the fact that the amplitude f 2 (a) is of first order in ω 2 , are corroborated up to first order.

A consistent second order theory
Since López's [10] first-order theory in amplitudes coincides with the classical Kopal's theory [8], we can assume that the f 2 (a) amplitude and its derivative d f2(a) da are of first order in ω 2 and the rest of the amplitudes f 2k (a) and its derivatives are of orders higher than the first.
As a consequence of López's work [10] and in order to prevent confusions, the expressions (11), (12) and (13) will be denote by K n , K 2 and W n , respectively. That is to say, where a 1 is the first root of the equation ρ(a) = 0.
To deal with the development of the self-gravitational potential Ω up to second order in the amplitudes, we will proceed to extend up to the second order the algorithm based on the numerical quadrature in [10].
The expressions for the radii vectors r and r ′ of the equipotential surface containing dm and dm ′ are given by r = a(1+Σ) and r ′ = a ′ (1+Σ ′ ), respectively. Let be (a, θ, λ) the Clairaut coordinates of this equipotential surface in the direction (θ ′ , λ ′ ) given by a(1 Expressing U n (6) and V n (7) in Clairaut's coordinates are obtained (42) To evaluate the integrals of the second member of 40, 41 and 42, we will proceed as follows: Equipotential surface S containing P is defined. Let be (a ′ , θ ′ , λ ′ ) the Clairaut coordinates of a point P ′ of the sphere containing P.
Then, in second order in amplitudes, we have On the other hand, approaching up to second order in amplitudes, it is obtained An improved algorithm of second order to construct consistent theories...
We define S 1 = Σ − Σ ′ y S 2 = Σ ′2 − ΣΣ ′ . Operating up to second order in the amplitudes, it is obtained (51) To evaluate (40) we will proceed as follows: Approximating up to second order in amplitudes the second integral of the second member of (53), by means of the numerical integration method of by means of the numerical integration of the trapezoid first order method, we obtain Since S 1 is of the first order in amplitudes, it is developed ∂r ′2−n ∂a ′ and ρ(a + aS 1 ) up to first order, taking into account that, in first order, f 2 (a ′ ) = f 2 (a) and (57) Substituting 57 in 54 and operating To evaluate the third integral of (53) to second order in the amplitudes we will use the rectangle method where ξ ∈ [a + aS 1 , a + aS 1 + aS 2 ]. Thus, approaching up to second order, we obtain a+aS1 a+aS1+aS2 ρ(a ′ ) ∂r ′2−n ∂a ′ da ′ = −(2 − n)a 2−n ρ(a)S 2 , for n = 2.
Taking into account (45) y (61)-(68); and then approaching up to second order in amplitudes the terms of Ω (5), it is obtained: U n r n + V n r −1−n = K n r n + W n r −1−n for n > 4.
Consequently, from (5) and (69)-(72) That is, even if U n = K n , U 2 = K 2 and V n = W n , as was just verified in (64), (62) and (68), it is proved that, up to second order in the amplitudes, the development of the self-gravitational potential (Ω) obtained by classical methods [8], coincides with the development calculated in this work (73).
Therefore, the method used in this section has permitted us to prove that the classical self-gravitational potential theory is correct up to second order.

Numerical example
In order to illustrate the efficiency of the Clairaut method, our proposal is to study, from a numerical point of view, the behavior of a radiative Sun-type star in slow rotation. To do this, we have first chosen as a state equation the polytrope model of order n = 3, a very appropriate model in this case [2], [3], [9]. A very wide-ranging modern study of polytropic models can be found in [5].
A polytrope of order n is characterized by the equation of state where P represents the pressure in an inner point, ρ is the density and K is a constant depending on the polytrope mass. If central density is represented by ρ c , total potential by Ψ, total central potential Ψ c and the next variable changes are done, then and For a Sun-type star with ν = 5.0 · 10 −3 we have K = 3.84116 · 10 14 , ρ c = 7.6309278 · 10 4 Kg/m 3 . For the specific case of an isolated mass at rest and therefore with spherical symmetry, the classic Lane-Emden equation is satisfied.

José A. López Ortí, Manuel Forner Gumbau, Miguel Barreda Rochera
In this work, the study of the equilibrium configurations of celestial bodies in uniform rotation is discussed, regardless of the two previous considerations. The method followed in this study consists of extending up to second order in ω 2 the algorithm of numerical quadrature referenced in [10].
This algorithm is based on the second-order approximation in ω 2 of the expressions U n and K n for n = 0, 2, 4 by an appropriate numerical quadrature technique.
In spite of the previous inequalities, by means of this algorithm has been proved that the development of the self-gravitational potential obtained in the classical theory up to second order in ω 2 , coincides with the one obtained in this work. This is so because the terms of the inner potential V n developments and those of the outer potential U n that differ from those of W n and K n , respectively, cancel one another when they are added to obtain the self-gravitational potential.
As has been shown in this work, as in [10] and [11], moving from a first order theory to a second order theory supposes a considerable increase in algorithmic complexity. Consequently, the number of calculations required to evaluate potentials beyond the second order is so large that it is extremely difficult to address them. All this reasons lead us to conjecture that the method exposed in this work are extendable to orders above the second.
Finally, to show the efficiency of Clairaut's method, we have proceeded to implement a numerical example using for that a radiative star of one solar mass with slow rotation. As is common in these cases, an index 3 polytrope has been used to model the problem.
The integration of the model has shown two difficulties being the first one due to the existence of a singular point at the origin and the second one caused the boundary conditions are canceled at the origin.
The first difficulty has been overcome by using analytical series developments, of greater degree than the numerical method used, around the origin.
The second difficulty has been solved by applying the shooting method.
In short, results shown in table 1 are consistent with the example used in this work.