Periodic orbits around the collinear equilibrium points for binary Sirius, Procyon, Luhman 16, α-Centuari and Luyten 726-8 systems: the spatial case

An investigation of three-dimensional periodic orbits and their stability emanating from the collinear equilibrium points of the restricted three-body problem with oblate and radiating primaries is presented. A simulation is done by using five binary systems: Sirius, Procyon, Luhman 16, α-Centuari and Luyten 726-8. Firstly, based on the topological degree theory, the total number of the collinear equilibrium points for the five binary systems were obtained and then, their positions were determined numerically. The linear stability of these equilibrium points was also examined and found to be unstable in the Lyapunov sense. An analytical approximation of three-dimensional periodic solutions around them was established via the Lindstedt–Poincaré local analysis. Finally, using the analytical solution to obtain starting orbits, the families of three-dimensional periodic orbits emanating from these equilibria have been continued numerically.


Introduction
In the circular planar restricted three-body problem, two point masses or primaries are fixed in a coordinate system rotating at the orbital angular velocity with the origin at the barycentre of the two primaries (as illustrated in figure 1). It is well known that in this rotating frame, there are five stationary points at which the infinitesimal body would remain fixed if placed there. Three of these stationary (or equilibrium) points lie on the line connecting the two primaries and are otherwise known as the collinear equilibrium points denoted as L 1 , L 2 , L 3 , while the other two which form triangular configuration with the primaries in the plane of primaries' motion are called triangular equilibrium points and are denoted by L 4 , L 5 .
Studies on periodic orbits around these equilibrium points have aided meaningful developments in the fields of celestial mechanics and space explorations. Researchers have utilized different approaches in order to examine periodic orbits around the equilibrium points and in particular, the collinear equilibrium points (Richardson 1980, Howell 1984, Kalantonis et al 2001, Lara and Pelàez 2002, Hou and Liu 2009, Liu et al 2014, Jiang 2015, Zotos 2015, Abouelmagd et al 2016. The stability of planar periodic orbits in the vicinity of the collinear equilibrium points either for in-plane or out-of-plane perturbations have been studied by Hénon (1965Hénon ( , 1973, Ragos and Zagouras (1991), Jain et al (2006). Recently, Singh et al (2016) examined periodic motions around the collinear points of the restricted three-body problem where the primary is a triaxial rigid body and the secondary is an oblate spheroid together with perturbation in the Coriolis and centrifugal forces.
Studies on periodic orbits with applications to the stellar systems are very important in Astrophysics. The masses of the stars in a binary system can be indirectly obtained from the calculation of their orbits. As a result, other stellar parameters like radius and density can be estimated. The masses of single stars can also be estimated by the determination of their empirical mass-luminosity relationship. Nagel and Pichardo (2008) gave a simple analytical formulation for periodic orbits in binary stars. Bosanac et al (2015) investigated periodic motions near a large mass ratio binary in the restricted three-body problem where stability analysis was used to evaluate the effect of the mass ratio on the structure of families of periodic orbits. For an investigation on periodic orbits around the triangular points with applications to the binary stellar systems: Kepler-34, Kepler-35, Kepler-413, and Kepler-16, Mia and Kushvah (2016) applied the Fourier series method to obtain them semi-analytically.
In the present work we consider the five binary systems Sirius, Procyon, Luhman 16, α-Centuari and Luyten 726-8 modeled as being radiating sources and sufficiently oblate in the framework of the spatial restricted threebody problem. In our investigation we firstly determine the number and positions of the collinear equilibrium points for all five binary systems and then examine their linear stability. We then apply the Lindstedt-Poincaré method in order to obtain analytical expressions up to second order terms for three-dimensional periodic orbits around the collinear equilibria. The analytical investigation has been also used to obtain starting points in order to compute numerically the families of three-dimensional periodic orbits emanating from these points. The stability of these orbits has also been examined.

Equations of motion and variation
The equations of motion of the infinitesimal body in the dimensionless synodic coordinate system Oxyz i.e., taking the units of mass, length and time such that the sum of the masses of the primaries, the distance between them and the gravitational constant are all unity, with the bigger and smaller stars of the binary systems having masses m 1 and m 2 , with the mass parameter being m m m , 2 1 2 m = + ( )and radiation pressure and oblateness parameters taken as q 1 , q 2 q i 1, 1, 2 i  = ( ) and A 1 , ) respectively, are (Singh and Ishwar 1999): where the potential function is given by: are the distances of the infinitesimal mass body from the two stars, respectively. In the six-dimensional phase space, equation (1) can be written in the form: where we have set: while the distances are now expressed by: 1 .
The Jacobi integral w.r.t. the equations of motion (2) is given by: The coordinates of the third body in the six-dimensional phase space depend uniquely along any solution on the initial conditions and the time, i.e.: x x x x x t i , , ; , 1, 2, , 6, i i 01, 02 06 and their partial derivatives with respect to the initial conditions satisfy the equations of variation (see, e.g., Markellos 1977, Jain et al 2006): The variational equations can be written as follows:

Numerical values of the physical parameters
In table 1, we present the physical parameters of the binary systems. The parameters M A and M B are the masses of the more massive and less massive stars in each binary system as compared to the mass of the Sun, with an exception of the binary Luhman 16 system where comparison is being to the mass of Jupiter. The luminosity of the binary systems denoted by L A and L B respectively is obtained from the relation (Mia and Kushvah 2016): where L S and M S are the luminosity and mass of the Sun. Radiation pressure has had a key effect on the formation of stars and shaping of clouds of dust and gases on a wide range of scales. The mass reduction factor is represented as q F F i 1 , 1 , 2 i p g = -= (F p and F g are the radiation pressure and the gravitational attraction forces being exerted by the binary systems on objects around them) or q i 1 , 1, 2 or on the basis of the Stefan-Boltzmann's law (Xuetang and Lizhong 1993) as: where M, L, and k are the mass, luminosity, and radiation pressure efficiency factor of a star. Also, a and r are the radius and density of the dust grain particles moving in the binary systems while A The determination of the exact number of roots of equation (6), at each one of the above mentioned open intervals, has been established by an approach based on the topological degree theory and can be briefly described as follows: if the function F(x): ] is two times continuously differentiable in this interval, then the total number N of roots of the equation F(x)=0 is obtained by the following scheme (Picard 1892, Kalantonis et al 2001, Singh andBegha 2011): where γ is a small positive real constant. For all the considered binary systems, i.e. for the parameter values corresponding to the specific systems, we have shown that each one of the aforementioned open intervals contains a unique real root; these roots correspond to the collinear equilibrium points L 1 , L 2 and L 3 , respectively. Note, that the involved integral in the above formula has been computed numerically by using Romberg integration. Since we have determined the exact number of the collinear equilibrium points we are able to solve numerically (6) so as to obtain their positions accurately. These are shown in table 2 for each binary system. Also, in figure 2 we present the effect of the oblateness coefficients on the positions of all collinear equilibrium points of the Sirius binary system while in figure 3 we show the effect of radiation factors on these points of the same system.

Stability of the collinear equilibrium points
We examine the linear stability of the collinear equilibrium points in the plane of motion of the primaries for each one of the five binary systems by solving numerically the following characteristic equation: arising from the linearization of the equations of motion around these equilibria. The coefficients 1 P and 1 Q in the above equation are given in the next section. The roots of equation (7) are presented in table 3 for each one of the five binary systems. As we can see from this table, for all systems and at each collinear equilibrium point four roots exist two of which are real and two pure imaginary. Therefore, the collinear equilibrium points are unstable in the Lyapunov sense. However, by eliminating the hyperbolic component of the general solution of the linearized system around the collinear equilibria, we are able to obtain analytical approximations of periodic motions around these points.

Three-dimensional periodic motion around the collinear equilibrium point
In order to study the motion of the infinitesimal body near any of the collinear equilibrium points L i , i=1, 2, 3, we write , a n d ,

Analytical approximation of periodic solutions
Expanding equation (8) into Taylor series up to second order terms we obtain the following system:   À = À + À I Figure 3. The effect of the radiation coefficient of the binary Sirius system on the collinear equilibrium points L 1 , L 2 (top row, left and right frame, respectively) and L 3 (second row).
The symbols 1 J and 2 J represent the signs of r x where ε is a small orbital parameter. Using them in system (9) we find that: which has non-zero solution only when 0.
which has been chosen in order to have initial velocity on the À axis, so as to obtain the motion out of the orbital plane. Also, by setting:

I R
From (13) and for t=0 we obtain the initial conditions of a periodic orbit in the form: x , , , , , , 0, 0, 0, 2 , . 14 t 0 0 11 12 2 2 22 a a e we b we while the value of the orbital parameter ε is arbitrarily chosen.
In figure 4, we present the analytical three-dimensional solution (blue curve denoted by A) emanating from the collinear equilibrium point L 1 of the α-Centauri binary system, as this is obtained from formula (14) of our analysis and by varying the orbital parameter ε in the range of values [0, 1.2]. For comparison reasons and for establishing the validity of our analysis, we also give in the same plot the numerical solution of the threedimensional periodic orbits (black curve denoted by N) emanating from the same equilibrium point as this is determined by the numerical integration of the full equations of motion (2). As we can observe from figures 4(a)-(c), our analysis is valid for small values of ε, in all the projections of the space of initial conditions x y z , , , , . À =     ( ) ( ) I R

Numerical approximation of periodic solutions
The figure eight shaped orbits occur as a result of the gravitational pull of each star of the binary system on the motion of the infinitesimal body. Since three-dimensional periodic orbits emanating from the collinear equilibrium points are of figure eight shape they are of double symmetry with respect to the Ox-axis and the Oxz plane. So, and for economy on the computations, it suffices to compute them at the quarter of their period. In particular, for the numerical determination of a three-dimensional periodic orbit of double symmetry w.r.t. the Ox-axis and the Oxz plane we integrate numerically the equations of motion with initial conditions of the form x x x , 0,0,0, , , ( ) i.e. the numerical integration starts on the Ox-axis and seek a perpendicular crossing of the orbit with the Oxz plane at which obviously the condition x x x x , 0,0,0, , 0 2 01 05 06 = ( ) is fulfilled. Therefore, we look for the following two periodicity conditions: Since two equations with three unknown components x x , 01 05 and x 06 of the initial state vector have to be satisfied we have to fix one unknown and apply well-known differential corrections procedures for the remaining two (see, e.g., Perdiou et al 2013). So, for choosing, e.g., x const., 06 = and by linearization of (15) we obtain the corrector system: The stability of a three-dimensional periodic orbit can be determined by integrating numerically the variational equation (4). Such an orbit will be stable if simultaneously the following conditions hold (Bray andGoudas 1967, Zagouras andMarkellos 1977): where V is the variational matrix. For stability of a three-dimensional periodic orbit in the restricted problem we also refer to Perdios (2007). x y , , plane.
In figure 5, we present the characteristic curves, in the space of initial conditions x y z , , ,   ( ) of the families of three-dimensional doubly symmetric periodic orbits emanating from all collinear equilibrium points of the α-Centauri binary system, as these were computed by the aforementioned procedure. In the right frames of this figure, we plot a typical figure eight member orbit for each family in the physical space while, in figure 6, where we plot the quarter of period versus the stability parameters P and Q as these were defined by (17), we give the stability diagram of the corresponding families. As we can see from this figure, stability of three-dimensional periodic orbits occurs only for families emanating from the equilibrium points L 2 and L 3 while we also observe that the period of these orbits goes to zero in the evolution of all the computed families.
This was the criterion used to stop computing them. In other words, such a member orbit with period near to zero has been considered as a termination point of any determined family. For example, the full period of the last computed three-dimensional periodic orbit of the family emanating from the collinear equilibrium point L 1 This member orbit is presented in figure 7 where we observe that it almost appears to exist on the Oxz plane since the values of the y coordinate are near zero.
In tables 4-6, we give initial conditions for all families of three-dimensional periodic orbits about the collinear equilibrium points of all five considered binary systems. In particular, we give the values of the orbit on the Ox-axis, i.e. x x, 1 = y x x,  the period up to its vertical intersection with the Oxz plane, i.e. at the quarter of the full period, as well as the value of the Jacobi constant C. We refer here to the work by Tsirogiannis et al (2006) where they have presented results for three-dimensional periodic orbits around the collinear equilibrium points of the restricted three-body problem in the case when the larger primary is a source of radiation and the smaller one is an oblate spheroid.

Conclusions
In this investigation, the primary and secondary bodies of the restricted three-body problem have been taken as oblate spheroids and radiation sources. The theory has been applied to five binary systems: Sirius, Procyon,   Luhman 16, α-Centuari and Luyten 726-8. The physical parameters of each binary system were obtained and used to calculate their respective mass parameters and mass reduction factors. The number and positions of the collinear equilibrium points of each binary system were obtained numerically by combining the topological degree theory and the numerical determination of roots of nonlinear algebraic equations. The effect of the respective oblateness and radiation coefficients on the collinear equilibrium points of the binary Sirius system were shown graphically. The linear stability of the collinear equilibrium points of the five binary systems was also examined and found to be unstable in the Lyapunov sense just like in the classical case. An analytical approximation of three-dimensional periodic motion around the collinear points was obtained by utilizing the Lindstedt-Poincaré method. The analytical solution was used in order to find suitable starting points for the numerical computation of the respective families of three-dimensional periodic orbits about the collinear points. The stability of the obtained periodic orbits for the binary α-Centuari system was also examined. Our results showed that only families emanating from L 2 and L 3 contain stable parts. Finally, we found that all the computed families, in their evolution comprised by member orbits which have decreasing period. Our termination criterion for their continuation was that the value of the full period of an orbit will be less than the value 0.0005. Using this criterion none of the determined families terminates at the physical plane which is not the usual termination of this kind of families in the restricted problem where they all go to a bifurcation with a coplanar periodic orbit, such as in the study by Tsirogiannis et al (2006) and Singh et al (2016). The choice to take the primaries to be sufficiently oblate bodies seems to be the reason where we obtain the different result in the evolution of three-dimensional periodic orbits emanating from the collinear equilibrium points of the binary α-Centuari system.