Combined effect of Poynting-Robertson (P-R) drag, oblateness and radiation on the triangular points in the elliptic restricted three-body problem

This study investigates the motion of a test particle around triangular equilibrium points in the elliptic restricted three-body problem (ER3BP) under the influence of the two oblate and radiating primaries having Poynting-Robertson (P-R) drag. It is observed that the position of triangular points of the problem is affected by oblateness, radiation pressure, eccentricity, semi-major axis and Poynting-Robertson (P-R) drag. The stability of these points is demonstrated analytically by the Routh-Hurwitz criterion. It is seen that they are unstable under the combined effect of involved parameters. The effect of these parameters on the position of triangular points is examined numerically using the binary systems, 61 Cygni and Archird. The results obtained by these binary systems can be used to broaden the scope of interest in astronomy, astrophysics, space science and celestial mechanics in general.


Derivation of equations of motion
Consider a rotating frame of reference 0xyz where 0 is the origin at the center of mass of the two primaries m 1 and m 2 (Szebehely 41 ).The x − axis lies along the line joining the two finite masses m 1 and m 2 .The y − axis lies in the plane of the orbits of the finite bodies and is perpendicular to the x − axis .The z − axis is perpendicular to the orbital plane of the finite bodies.Let r 1 and r 2 be the respective distances of m 1 andm 2 of the infinitesimal body m from the bigger and smaller primaries positioned at (x 1 , 0, 0) and (x 2 , 0, 0) respectively as shown in Fig. 1.Then, the equations of motion of the tiny mass in the sidereal coordinate according to Singh and Umar 8 , can be given as; where k 2 is the gravitational constant, we now rotate the rectangular frame of reference about the z − axis so that the primaries always lie on the x − axis .The new coordinates are related to the old ones by transformation formulae.
where υ is the true anomaly.Then, the equations of motion are given in Eq. (1) in the synodic reference system are given by (1) x = X cos υ − Y sin υ, y = X sin υ + Y cos υ, z = Z www.nature.com/scientificreports/Considering the force function U in Eq. ( 2) we have that where R is the distance between the primaries.We change the equations of motion from the true anomaly to the eccentric anomaly, we use the following relation by Singh and Umar 8 ; where e, E, andυ are respectively the eccentricity, eccentric anomaly and the true anomaly of the orbit of the primaries.Then from Eq. ( 5) we have Differentiating Eq. ( 5) for time t we obtain By Kepler's equation, we obtain where n and T are, the mean motion and the period of the elliptic orbit respectively.Then, the distances between the primaries are given as; where ' a ' is the semi-major axis of the elliptic orbit of m 2 around m 2 .Differentiating Eq. ( 8) With t and consider- ing Eq. ( 9) we obtain  www.nature.com/scientificreports/using Eqs.( 6), ( 9) and (10) in Eq. ( 7) yields With the introduction of dimensionless pulsating coordinates using the following transformation by Singh and Umar 8 , Differentiating Eqs.(11) with t and denoting the differentiation with respect to E by the prime, we obtain Also, Considering the values of X, Ẋ, Ẍ, Y , Ẏ , Ÿ , Z, υ and ϋ in Eq. ( 3) we obtain Transforming U in terms of the new coordinate system (ξ , η, ζ ) denoted by U , we have Comparing Eqs. ( 12) and ( 13) we obtain With the force function (1 − e cos E) 2 ; ϋ = −2en 2 1 − e 2 1 2 sin E ρ 4  and ρ = R a Vol.:(0123456789)  15) that the force function depends on the eccentric anomaly E , that is it varies for different positions of the primaries.We consider the average force function as; Now, integrating Eqs. ( 14) With respect to E and averaging, the system of equations reduces to the following; With the force function Equations ( 17) and ( 18) represent the equations of motion of the general elliptic restricted three-body problem (ER3BP).

Force due to radiation
Since the radiation pressure force, F p changes with distance by the law as the gravitational attraction force F g and acts opposite to it, it is possible that this force may lead to a reduction of the mass of the massive particle.Since this reduction of mass depends on the properties of the particles, then the resulting force on the particle is F g is the mass reduction factor.We denote the radiation factor as q i (i = 1, 2) for both primaries such that 0 < 1 − q i ≪ 1 , hence, instead of mass m 1 , m 1 q 1 will appear on the force function, and instead of m 2 , m 2 q 2 will appear on the force function.

Force due to oblateness
According to McCuskey 42 , the force due to oblateness of the bigger primary with mass m 1 and smaller primary m 2 is given by where AE and AP are respectively the dimensional equatorial and polar radii of the primaries.Let F be the potential due to oblateness, then; Comparing Eqs.(19) and (20) and then integrating with respect to r 1 and r 2 we have Therefore, the force function due to the combined effects of oblateness and radiation on both primaries becomes;

Mean motion
According to Singh and Umar 8 , the distance between the primaries in the elliptic case is given as; r = a(1−e 2 ) 1+e cos µ and the mean distance between them is given by Vol:.(1234567890) www.nature.com/scientificreports/Szebehely 45 asserted that the orbits of m 1 and m 2 with respect to the center of mass, with semi-major axes a 1 = am 1 and a 2 = am 2 respectively have the same eccentricity in the two-body problem.That is; Adding the Eqs.(22) we have Now, we choose the unit of time such as to make k 2 = 1 and the unit of mass such that the sum of the mass of the finite bodies is taken as unity.For this, we take

Scientific
are the masses of the bigger and smaller primaries and their coordinates are and (1 − µ, 0, 0) respectively.A i << 1 (i = 1,2) depict the oblateness of the bigger ( A 1 ) and smaller ( A 2 ) primaries,q 1 , q 2 are the radiation factors for both primaries; a is the semi-major axis of the orbit; r i (i = 1,2) are their distances from the infinitesimal mass; e is the eccentricity of the orbits Where µ is the mass ratio.Then m 1 = 1 − µ .Also, ξ 1 = −µ and ξ 2 = 1 − µ .Considering these terms in Eqs. ( 17), ( 18), ( 21) and ( 23) we obtain the equations of motion in the elliptic restricted three-body problem with oblate and radiating primaries as; with the force function The mean motion, n, is given as

Force due to Poynting Robertson drag
Following Singh and Amuda 39 , the components of Poynting-Robertson(P-R) drag of the bigger and smaller primaries can be written as; Simplifying Eq. ( 28), the equations of motion of the infinitesimal mass (test particle) in the elliptic restricted three-body problem when both primaries are oblate and luminous with Poynting-Robertson drag can be written using the dimensionless variables and a barycentric synodic coordinate system ( ξ,η , ζ ) as: where primes denote differentiation with respect to the eccentric anomaly.where C d characterizes the dimensionless velocity of light and q is a factor known as the radiation effect (Radzievsky 32 ).Thus, these equations are affected by the oblateness, radiation pressure, semi-major axis, eccentricity, and Poynting-Robertson drags of the primaries.

Location of triangular equilibrium points
The equilibrium points are those points at which the velocity and acceleration are zero.Therefore, these points are the solutions to equations of the system (29). where This implies that they are solutions of equations From the last equation one can obtain ζ = 0 This shows that the coplanar stationary or equilibrium points exist.
The location of triangular libration points is obtained by the solutions of the first two equations taking into consideration η = 0 , ζ = 0.
Thus, from Eq. ( 32) we have For a triangular point, the motion occurs in the orbital plane ( ξη ), on solving the two Eq.( 33), we obtain; Using Eq. ( 34) in the photogravitational elliptic case, that is, when oblateness and P-R drag are absent, 34) reduces to Which provide us If we accept that both primaries are oblate and radiating with P-R drag, then Let ε 1 and ε 2 be the perturbations caused by the presence of radiation, oblateness, and P-R drag, then the values of r 1 and r 2 will change slightly by ε 1 , ε 2 respectively so that where |ε 1 |, |ε 2 | << 1.
On using only linear terms in A 1 ,A 2 , ande 2 and neglecting their products, Eq. ( 26) becomes; If A 1 , A 2 = 0, we get Substituting the value of Eq. ( 35) in Eq. ( 36), we obtain Using equations.( 36) and (37) in Eq. ( 34) and neglecting second and higher-order terms of small quantities; Substituting the values of ε 1 and ε 2 from Eq. (38) in Eq. (37), we obtain; then Substituting values of r 2 1 and r 2 2 from Eq. ( 39) into Eq.( 31) and then ignoring the second and higher order terms for ξ and η , we obtain

Linear stability of triangular equilibrium points
The stability of the test particle (infinitesimal mass) around an equilibrium point describes the coordinates of a given equilibrium point denoted by (ξ 0 , η 0 ) and small displacements ( σ , ω ) from the point.So that Then the variational equations of motion are given as Then the characteristic equation that relates to ( 42) is given as where Solving the second partial derivatives at the triangular libration point, we have Substituting the values 17), the system of equations becomes Now we wish to discuss analytically the stability of motion using the Routh-Hurwitz criterion 43-46 .For this, we first evaluate the principal diagonal minors of the Hurwitz matrix composed of the coefficients of polynomials (43).Substituting the values of a 1 ,a 2 ,a 3 , and a 4 from Eq. ( 46) in Hurwitz diagonal minors, we have Which on neglecting higher order terms, reduces to Since Hurwitz minor matrices are not all greater than zero, therefore, by the Routh-Hurwitz criterion, the triangular points are unstable.Even the value of a 3 shows that the triangular points are unstable.

Numerical application
The triangular points can be obtained numerically.These are done for two binary systems 61 CYGNI and ARCHIRD, which have oblate and radiating primaries with Poynting Robertson drag.The numerical data relating to the system stated are shown in Table 1 below.
Using the value q = 1 − AκL rρM by Stefan-Boltzmann's law we compute the radiation pressure factor q by taking κ = 1 , M as the mass of the star, L as the luminosity of a star, r as the radius, ρ as the density of a moving body and κ as the efficiency factor of a star; A = 3 16πCG is a constant.In the C.G.S. system, A = 2.9838 × 10 −5 , r = 0.02cm and ρ = 1.4gcm −3 for some dust grain particles in the systems.The mass ratio is µ = m 2 m 1 +m 2 .The dimensionless velocity of light (C d ) for the binary stars according to 39,47 is computed as Au Here A u is the binary separation of the primaries, c is the velocity of light, and γ is the gravitational constant; m 1 , m 2 are masses of the bigger and smaller primary respectively.In C.G.S. system, c = 2.99792458 × 10 10 cms −1 and γ = 6.6743 × 10 −8 cm 2 g −1 s −2 .
Using these data, we compute the locations of triangular points, showing the effects of the parameters on the binary systems 61 CYGNI and ARCHIRD.The effect of the oblateness of A 1 in the absence of A 2, and then with a constant A 2 is shown in Tables 2 and 3 respectively.Then, in the absence of A 1 and keeping A 1 constant, the effect of A 2 is examined (Table 4 and 5 respectively) using the Mathematica 12.1 software.We show the effects of radiation pressure ( q ), eccentricity ( e ), and semi-major axis ( a ), on the triangular points (Tables 6, 7, 8, and  9) for arbitrary values.

Discussion
We have studied the effects of oblateness, radiation pressure, eccentricity, and semi-major axis of the orbits of primaries together with P-R drag, on the motion of a test particle in the neighbourhood of triangular points in the ER3BP described by Eqs. ( 29)- (31).Equations ( 40) and ( 41) give the positions of triangular points similar to that of 6 in the absence of zonal harmonics J 4 and Poynting Robertson (P-R) drag (for This also validates the position of triangular points for other generalizations by suppressing certain parameters.
In the absence of radiation pressure force, oblateness and P-R drag effects, it reduces to the classical case of Szebehely 48 , i.e., when the primaries are spherical and non-emitters of radiation and without P-R drag.The location or position of the triangular points is computed numerically and shown in Tables 2, 3, 4, 5, 6, 7, 8 and 9, with each figure drawn and highlighted on the corresponding table.
It is seen from Tables 2 and 3 that, ξ increases and η decreases with A 1 in the case of 61 Cygni, while in the case of ARCHIRD, ξ increases with A 1 and η does not exist.It converts into collinear movement through ξ axis only as shown in Figs. 2, 3 and 4 respectively.In Tables 4 and 5, ξ decreases, and η decreases with A 2 in the case of 61 Cygni, while in the case of ARCHIRD, it moves towards the origin along the ξ axis only and η does not exist.It converts into collinear points as shown in Figs. 3 and 5 respectively.Decreasing eccentricity while radiation pressure, oblateness, and P-R drag remain constant causes a shift away from the origin in the negative direction and away from the line joining the primaries.This is shown in Table 7 and Figs.6, 7, ξ increasing in the negative direction and η moves away from ξ axis.It can be seen that for high eccentricity (e = 0.88), the triangular points  do not exist, thereby making it collinear.By comparing the present work with Singh and Amuda 39 , for circular restricted three-body problem (where e = 1) the triangular point ceases to exist.This is not the case with a decrease in the semi-major axis (Table 8 and Fig. 8).We also observe that, for a small value of the semi-major axis (a = 0.2), the triangular points cease to exist, thereby making it collinear.Also, increasing the radiation pressure with constant oblateness decreases ξ towards the origin and η increases away from the ξ axis as shown in Table 6 and Fig. 6.A change in the position of the triangular equilibrium points due to P-R drag forces has been found and it was observed that the system causes no change in the position of triangular points (as shown in Table 9 and Fig. 9).
In analyzing the stability of the triangular equilibrium points analytically, we have applied the Routh-Hurwitz criterion.Since Hurwitz minor matrices are not all greater than zero, therefore, by Routh-Hurwitz criterion, the triangular equilibrium points are unstable.Particularly, when A 1 = A 2 = W 1 = W 2 = e = α = β 1 = β 2 = 0 , then a 1 = 0, a 2 = 1, a 3 = 0 and a 4 = 27  4 µ(1 − µ) Szebehely 45 .This gives the classical case.The roots of the characteristics Eq. ( 43) are given in Eq. ( 47) where values of a 1 ,a 2 ,a 3 and a 4 are provided in Eq. (46).These equations all depend on the parameters of the perturbing factors.
It is observed that the nature of the value a 3 as compared to the condition necessary for stability by Routh Hurwitz, shows that the motion of an infinitesimal mass (test particle) around triangular points is unstable.The numerical values of the four roots of the characteristic Eq. ( 43) as given by Eq. ( 47) are shown in Tables 10, 11, 12. Tables 10 and 11 show that the system is stable in the presence of oblateness.But in the presence of both oblateness and P-R drag, it is unstable.Table 12, also shows that in the presence of radiation, the system is stable, but in the presence of PR-drag, whether radiation is present or absent, the system remains unstable.

Conclusion
We have studied the positions and stability of triangular equilibrium points in the elliptic restricted three body problem (ER3BP) under the speculation that the primary bodies are both oblate and luminous with Poynting Robertson (P-R) drag.It is seen that the equations of motion for the stated problem are affected by the oblateness, semi-major axis, radiation pressure, eccentricity of the orbits, and Poynting Robertson (P-R) drag.
The stability of the system has been investigated using the analytic and numerical approach.The locations and stability of triangular points are significantly affected by the aforesaid Parameters.The system is stable in the presence of oblateness or radiation or both of any of the primaries or both primaries.But it becomes unstable when the P-R drag of any or both of the primary bodies come into action, whether one or both of the primaries are oblate or luminous or both.Thus, Poynting Robertson (P-R) drag ruins the stability.Showing the effect of different values of P-R drag on the triangular points with A 1 = 0.02, A 2 = 0.04, q 1 = 0.79 , q 2 = 0.89,e = 0.8,a = 0.8 Table 10.Showing the roots of characteristic equation for 61 CYGNI with different values of oblateness and P-R drag with q 1 = 0.9997676, q 2 = 0.999856, a = 0.628, and e = 0.4900.

( 8 )Figure 1 .
Figure1.Bigger primary m 1 and smaller primary m 2 moving in an elliptic orbits about their common centre of mass unperturbed by the infinitesimal mass m.

Figure 2 .
Figure 2. Showing the effect of different values of A 1 on the triangular points of 61 CYGNI (left diagram) and ARCHIRD (right diagram) in the absence of A 2 .

Figure 3 .
Figure 3. Showing the effect of different values of A 2 on the triangular points of 61 CYGNI (left diagram) and ARCHIRD (right diagram) in the absence of A 1 .

Figure 4 .
Figure 4. Showing the effect of different values of A 1 on the triangular points of 61 CYGNI (left diagram) and ARCHIRD (right diagram) with constant A 2 .

Figure 5 .
Figure 5. Showing the effect of different values of A 2 on the triangular points of 61 CYGNI and ARCHIRD with constant A 1 .

Table 1 .
Showing numerical compilation for the binary systems used.Sources: Singh and Ashagwu 6 and Singh and Amuda 39 .

Table 2 .
Showing different values of A 1 in the absence of A 2 on the triangular points of 61 CYGNI and ARCHIRD.

Table 3 .
Showing different values of A 1 with constant A 2 on the triangular points of 61 CYGNI and ARCHIRD.

Table 4 .
Showing different values of A 2 in the absence of A 1 on the triangular points of 61 CYGNI and ARCHIRD.

Table 5 .
Showing different values of A 2 with constant A 1 on the triangular points of 61 CYGNI and ARCHIRD.