Two algorithms to construct a consistent ﬁrst order theory of equilibrium ﬁgures of close binary systems

One of the main problems in celestial mechanics is the study of the shape adopted by extended deformable celestial bodies in its equilibrium conﬁgura-tion. In this paper, a new point of view about classical theories on equilibrium ﬁgures in close binary systems is oﬀered. Classical methods are based on the evaluation of the self-gravitational, centrifugal and tidal potentials. The most common technique used by classical methods shows convergence problems. To solve this problem up to ﬁrst order in amplitudes two algorithms has been developed, the ﬁrst one based in the Laplace method to develop the inverse of the distance and the second one based on the asymptotic properties of the numerical quadrature formulas.


Introduction
The main aim of this work is to develop a consistent theory to determine the self-gravitational potential at an arbitrary point of each component of a detached close binary system in a state of slow rotation, when its equilibrium configuration has been reached.This means that the binary system is observed from a coordinate system rotating in association with the studied component, once the hydrostatic equilibrium is reached.The configuration of hydrostatic equilibrium is equivalent to a state of rigid rotation, which corresponds to the minimum potential [4], [8], [7].In accordance, among others [4], [9], [7], this state implies the identification of equipotential, isobaric and isopycnic surfaces.
The mechanisms by which a close binary system reaches its equilibrium state have been studied among others by [2], [3].This process implies circularizing the orbit and synchronizing the rotation and translation movements.Let M 1 and M 2 be two masses and R the distance between the center of mass of its components.The equilibrium state implies Kepler's third law The state of the close binary system can be modeled by the equations: where P is the pressure, ρ is the density, Ψ is the total potential, △ is the Laplace operator, G is the constant of universal gravitation and − → ω is the angular velocity of the system.To integrate these equations in a general case of mass distribution, a state equation, P = P (ρ), relating pressure with density is required.
To assess the full potential is necessary to calculate the self-gravitational potential Ω, the centrifugal potential V c and the tidal potential V t .
The physical effects produced by the disturbing potentials (centrifugal and tidal) are deformations of the equipotential surfaces of the components with respect to the spherical symmetry that stars at rest would possess.
To study the potential of a generic point of the close binary system, it is convenient to define the coordinate system OXY Z. On it, O is the center of mass of the primary component, OX is the axis rotating in the direction of the secondary component, OZ is an axis in the direction of the angular velocity of the system and OY is defined in a way that OXY Z may form a direct trihedral.The value of V c and V t in this reference system is described by Finlay [4], Kopal [8], [9] and López [10], [12].
The solution of the problem can be addressed on the basis of the well-known equations determining the structure of a self-gravitatory fluid, initially isolated, which leads to a solution with spherical symmetry.
From this initial solution, the general problem includes centrifugal forces due to rotation and tidal forces generated by the secondary component.This problem can be analyzed from the solution of an unperturbed problem using the perturbation theory, successive approximations of the density, pressure and temperature data, and the definition of the shape of the component as the lower surface containing it.
An alternative treatment, due to Clairaut [7], is to consider that density, pressure and temperature of each component are only function of the total potential, so that the complete solution of the problem is equivalent to determine the equipotential surfaces.An advantage of this procedure is that pressure does not appear explicitly in the basic equations which determine the equipotential surfaces.The equations depend on the density, remaining invariant with the form it takes.
In the Clairaut method, the figure of each component is defined by the lower equipotential surface containing the component.The study of the close binary systems has been addressed using the Clairaut method by several authors [9], [10], extending the technique used by Finlay and Kopal, among others, to the problem of equilibrium figures for rotating isolated bodies.Kopal [9] studies the amplitudes of the deformations with respect to the sphere by adding the terms of rotation, tide and the rotation-tide interaction.On the other hand, Lopez [10] studies the problem globally, up to the first order, using as a basis for the developments of equipotential surfaces the spherical functions.
Unfortunately, the method used by these authors [9], [10] offers some problems regarding the convergence of the used developments and its solution depends on the desideratum of Laplace, which has not been yet demonstrated.In case of rotating figures, these problems have been studied by López [12].The study shows inconsistencies in Kopal's approach.To overcome these disadvantages, López [12] introduces two new algorithms which do not require compliance with the desideratum of Laplace.López shows that, up to first order, the inner and outer terms of the self-gravitational potential given by the Kopal theory are wrong.However, the full potential up to the first order coincides with the one obtained by Kopal.This demonstrates that the results of Kopal theory are consistent because they do not depend on Laplace desideratum.
This paper focuses on extending the results obtained for rotating deformable bodies showing rotational symmetry to the most general case of the close binary systems.In this case we show that the terms of the self-gravitational potential calculated using classical methods are not correct.Despite this, it is shown that the total potential, up to first order in the small amount ω 2 , coincides with that obtained by classical methods.In the case we are treating, spherical functions and their products are involved, which makes it far more complex than the problem considering only rotation, because in this case only the Legendre polynomials and their products must be considered.
In this section, the general backgrounds of the problem are explained.On them, the inconsistency of the classical theory due to its reliance on the unproven hypothesis known as the desideratum of Laplace is shown.It is also shown the need to build new methods to solve the first order theory for the close binary systems which are independent of the Laplace desideratum.
Section 2 describes further the classical method proposed by Finlay and Kopal, based on the method of Clairaut as well as its extension for close binary systems.In this section, difficulties offered by the classical theory due to the dependence of this theory with the unproven desideratum of Laplace are also shown.
In section 3, a new algorithm based on an extension of the method used by López [12] for close binary systems is shown.This method is based on a modified Laplace method for calculating the inverse of the distance between two bodies moving in elliptical orbits.This method obviates the need of depending on Laplace desideratum to evaluate, up to first order, the self-gravitational potential in a point of a component of a close binary system.
Section 4 introduces a second algorithm that extends the technique developed by López [12], to close binary systems.This method, consisting on a generalization of the Wavre method, permits to obtain the rigorous calculation of the different terms of the self-gravitational potential of each of the components.These terms do not coincide with those obtained by the classical theory but, nevertheless, it is also shown that, up to first order in ω 2 , the self-gravitational potential coincides with the one obtained by classical theory.
In section 5, the main conclusions of this work and its comparison with the results obtained by the classical theory are discussed.

Classical theory of close binary systems
In this section, the classical theory about equilibrium figures in close binary systems due to [4], [8], [9], [10] is presented.These theories are based in the use of the Clairaut coordinates system.Let (r, θ, λ) be the spherical coordinates of a point of the primary component.These coordinates are related to Cartesian coordinates (x, y, z) by The Clairaut coordinates system (a, θ, λ) is related to the spherical coordinates system by r = r(a, θ, λ).The total potential Ψ depends only on the a coordinate and then, for each constant a, the surface r = r(a, θ, λ) is an equipotential surface.
The value of a is defined by the following condition: the mass of the component enclosed by the equipotential surface r(a, θ, λ) coincides with the mass contained in the sphere of center O and radius a. Obviously, in the case of equilibrium, the density ρ depends only on the coordinate.To assess the total potential at a point contained in the primary component, it is necessary to obtain the self-gravitational potential Ω, the tidal potential V t and the centrifugal potential V c .
Let P be a point in spherical coordinates (r, θ, λ) of the system primary component.Let r 1 be the smallest sphere centered at O containing the primary component and dm ′ the mass element of a point P ′ with coordinates (r ′ , θ ′ , λ ′ ) of the component.
The self-gravitational potential is given by where ∆ is the distance between the points of coordinates (r, θ, λ) and (r ′ , θ ′ , λ ′ ), and The inverse of the distance can be obtained through the classical development where γ is the angle between the vectors − − → OP and − − → OP ′ .By the theorem of addition of spherical harmonics in real form, we know The self-gravitational potential Ω is evaluated by where where U is the potential due to the outside part of the sphere of radius r and V the potential due to the inside part of this sphere.
To evaluate these integrals it is convenient to replace 1 ∆ by (6), and so we have where Taking into account that the Jacobian of transformation between spherical and Clairaut coordinates (a ′ , θ ′ , λ ′ ) is ∂r ′ ∂a ′ , we have and following [4], [9] we can write the last equations as: where log denotes Naperian logarithms, and a(θ ′ , λ ′ ) is the value of the parameter a corresponding to the equipotential surface containing the point of spherical coordinates (r, θ ′ , λ ′ ), and a 1 the value of Clairaut coordinate a on the equipotential surface determitated by the shape of the component.Classical theory assumes that developments (11) and ( 13) are equivalent.However, this assumption is not true because ( 11) is calculated on the sphere and ( 13) over the equipotential surfaces.
To evaluate the last integrals, we consider [7] , where ξ(a, θ, λ) is a small quantity.In the other hand, we assume that the development for r converges: The functions f n,m (a) are called amplitudes and the functions Y n,m (θ, λ) are the spherical functions.The spherical functions are an orthogonal and complete system on the sphere.To manage the products of spherical functions we van see [1], [6], [5].The absolute value amplitude functions |f n,m (a)| increases with a [7], [13], [8].
In the case of an isolated body in rotation [8], [9], this development can be reduced as r = a(1 f 2n (a)P 2n (cos θ)).In the case of a close binary system by symmetry reasons, we have f n,m (a) = 0 if m < 0 or n + m is odd.
To evaluate the integrals (13), it is convenient to approach r ′k and log r ′ up to first order in amplitudes, and so: Replacing (15) in ( 13) and taking into account (10), we obtain the development [10] where, for n = 0, and On the other hand, the function V0 G represents the mass M (a) contained in the sphere of radius a which coincides with the enclosed in the equipotential surface r(a, θ, λ).In consequence Developing now r 3 up to second order and replacing the product of spherical harmonics by their developments, it is obtained From this, we can deduce that the amplitude f 0,0 (a) is of second order for amplitudes.
Moreover, if we assume the condition that the center of coordinates is in the center of mass of the main component, then From the first equation of (23), we have Moreover, the centrifugal potential and the tidal potential in first order [11] are expressed respectively by Notice that up to first order in amplitudes the disturbing potential V c + V t due to the third Kepler law, only contains terms in Y 0,0 , Y 2,0 , and Y 2,2 .The total potential Ψ = Ω + V c + V t is given by: And so Ψ(a) = Ψ 0,0 (a) and Ψ n,m (a) = 0 if n + m = 0. Following a similar way as Kopal [8], [9] only the amplitudes f 2,0 (a) and f 2,2 (a) are of first order in ω 2 .
Unfortunately the classical method involves several inconsistencies.The values of U n and V n given by ( 11) and ( 13) are not the same as will be later demonstrated, because in the first form the domain of integration they are limited by sphere centred in O containing the point P and in the second form this one is the region limited by the equipotential surface containing P .
An alternative method to solve this problem used by Liapunov [7] can be describe as: let us define V a as the region limited by the equipotential surface containing the point P with Clairaut coordinates (a, θ, λ) and be V e the region of the component of the equipotential surface of P .The self-gravitational potential Ω can be decomposed as and replacing we have where At this point it is important to mention the equivalence of equations ( 13) and (34).
Let's remember that in order to use (33), the assumption of the non-demonstrate Laplace desideratum (can be found in [7] and [13]) is necessary to be able to accept the convergence of the development (33).These hypotheses assume that by replacing the first part of (32) in the first term of (31) and the second part of (32) in the second term of (31), convergence is maintained.

First order analytical algorithm
In the classical method it is necessary to evaluate the potential through (33), but unfortunately, the series defined by (6) do not converge in the layer defined by r ∈ [r min (a), r max (a)] with To solve this problem a method based on the Laplace development of the inverse of the distance between two planets [13] will be used.López [12] apply this technique to obtain the inverse of the distance between two points in the rotational case, and we propose to extend this technique to the case of two points of a component of a close binary system.To do that, we define The inverse of the distance ∆ −1 = (r 2 − r ′2 − 2r r ′ cos γ) −1/2 between dm and dm ′ , ∆ −1 = D(r, r ′ ), can be obtained to first order in amplitudes by means of a Taylor development up to first order in amplitudes: Y n,m (θ, λ), up to first order in amplitudes we have: From ( 36) and ( 37) we obtain where The mass element dm ′ is expressed by (12).From (37) From ( 12) and (39), the mass element dm ′ expressed in first order terms is obtained Outside, we have: The self-gravitational potential Ω is expressed as the sum of the external potential and the internal potential where Now, we are going to develop the internal potential W up to first order.To do that we will use the development in the interior of D(a, a ′ ): From ( 38) and ( 46) From ( 40), ( 45) and ( 47) Then, from (48) and ( 49) we obtain Consequently, from (50) and ( 51) In the same way that W has been calculated, K development is obtained and so Developing (64) up to the first order in f p,q (a) and f p,q (a ′ ) On the other hand Approximating (66) up to the first order in f p,q (a), it stays From ( 65) and ( 67) On one hand, from the addition theorem of spherical harmonics in real form, and on the other, according to the orthonormality of spherical harmonics, where δ i,j is Kronecker's delta.
From ( 63), ( 68) and (69), and considering, in this work, q ≥ 0, m ≥ 0, it stays ) In a similar way From (34) and taking into account (70), ( 71), (72) we obtain and for n = 0 In short, we have come to show that both, the development of the outer potential (Un) and the development of the inner potential (V n ), are formed by adding up the development of the outer potential of classical theory (K n ) plus others terms (75) and by adding up the term of the inner potential of classical theory (W n ) plus other terms (75), respectively.
Given a From ( 76), ( 77) and (78), and approximating f n,m (a) up to the first order in amplitudes, we obtain which is the expression of the self-gravitational potential in its classical form.

Concluding Remarks
This research proves that Finlay and Kopal method used by the classical theory reveals inconsistencies already in the first order.These inconsistencies are derived from the assumption, relying on Laplace desideratum, that the developments ( 11) and ( 13) are equivalent.
In order to avoid confusions, in this work we have redefined the name of the developments displayed in (13), so that U n and V n are renamed as K n and W n , respectively (34).
By using a method that does not presuppose the Laplace desideratum and from the developments displayed in (11), U n is proved to be different from K n and V n is proved to be different from W n , as is brought to light in (75).
Moreover, the method used by Liapunov to solve this problem depends on the desideratum of Laplace, which has not been demonstrated.
López [12] obtained a solution in first order in ω 2 for the case of rotating bodies.This solution was achieved by using two independent algorithms: one based on the development of Laplace for calculating the inverse of the distance and another based on the calculation of the inner and outer terms of self-gravitational potential by approximation up to the first order, using the asymptotic properties of a numerical quadrature formula.
To extend these algorithms to the more complex case of close binary systems, a package of codes allowing the formal manipulation of the product of spherical harmonics to obtain the result as a linear combination of them has been developed.The source codes of this package are accessible under GPL version 4 conditions.
In section 3, the self-gravitational potential at an arbitrary point P of the primary component is calculated as the sum of the potentials created by the mass of the portion of the inner and outer component to the equipotential surface containing that point.This calculation is made by approximation, up to the first order, of the inverse of the distance between the point P and two points, one inside and one outside the equipotential surface passing through P .This approach up to the first order, which coincides with the one developed by the classical theory, is achieved using the package codes created to obtain the product of spherical harmonics as a linear combination of them and making use of the orthogonality properties of the product of spherical harmonics.It must be emphasized that the method used for calculating the self-gravitational potential is independent of Laplace's desideratum.
In section 4, the values of the terms U n and V n of the self-gravitational potential are obtained by an approximation, up to the first order, of the sphere containing the point P in Clairaut coordinates.These values have been obtained as the sum of K n , W n and the value of the integral in the small region between the sphere and the equipotential surface passing through the point P .The above integral has been obtained by exact numerical quadrature up to the first order.All of these elements permit showing that the values of U n and V n given by the classical theory are incorrect.However, the approximation up to the first order of the self-gravitational potential Ω coincides with that obtained in classical theory.
Consequently, in this work the results, up to the first order, of the classical theory about close binary systems are consistently demonstrated, that is to say, without turning to Laplace desideratum.
Finally, we propose to extend the obtained results to higher orders than the first one, both in the case of only rotation as in the case of the close binary systems.

As |f 1 , 1
(a)| is increasing, then f 1,1 (a) and d f 1,1 (a) d a have the same sign.Therefore, d a 4 f 1,1 (a) d a has the same sign as f 1,1 (a).Furthermore, as ρ(a) > 0 ∀ a < a 1 , we can deduce that, in first order, f 1,1 (a) = d f 1,1 (a) d a = 0.That is to say, order of f 1,1 (a) is of an order higher than the first.