On the triangular points within frame of the restricted three–body problem when both primaries are triaxial rigid bodies

Restricted three–body problem (RTBP) plays very important role in celestial mechanics and space science. [12, 15 and 23] are very good books in celestial mechanics which explains importance of RTBP in space dynamics. Classical RTBP is explained in detailed in [15]. [2, 3, 6, 8] have studied RTBP with different perturbations like solar radiation pressure, oblateness, air drug etc. With both primaries as point masses. [7] have studied numerical integration with Lie series in the case of RTBP.

It is well known that the classical planar restricted three body problem (CRTBP) possesses five liberation or stationary points. These points are known as Lagrangian points. Out of these five points three points are collinear which are unstable where as two points are triangular which are stable in nature. [1, 4, 5, 11] studied existence and stability of these Lagrangian points for the perturbed RTBP. Recently, [16, 17, 18-19] studied different family of periodic orbits and its stability using Poincaré surface section for perturbed RTBP. In recent times many perturbing forces i.e., oblateness and radiation forces of the primaries, Coriolis and centrifugal forces, variation of masses of the primaries and of the infinitesimal mass etc., have been included in the study of the restricted three body problem.

For the case, where the bigger primary is an oblate spheroid whose equatorial plane coincides with the plane of motion, [22] have studied the stability of the liberation points. A similar problem has been studied by [13]. [10] have studied existence and stability of the equilibrium points of the triaxial rigid body which is moving around another triaxial rigid body. [9] have studied the non-linear stability of triangular point L₄ in the RTBP when the bigger primary is a triaxial rigid body with its equatorial plane coincident with the plane of motion.

[14] have studied the problem when the smaller primary is a triaxial rigid body. Also [20, 21] have studied the problem when both the primaries are triaxial rigid bodies in the case of stationary rotational motion (θ_i = ψ_i = φ_i = 0). In this paper we consider the restricted three body problem when both the primaries are triaxial rigid bodies in two cases of stationary rotational motion: (θi=ϕi=π2,ψi=0) $\left(\theta _{i} = \phi _{i} = \dfrac{\pi}{2}, \psi _{i} = 0 \right)$ and (θi=0,ψi+ϕi=π2) $\left(\theta _{i} =0, \psi _{i} +\phi _{i} =\dfrac{\pi}{2}\right)$ . Also for Euler’s angles (θi=π2,ϕi=ψi=0) $\left(\theta _{i} = \dfrac{\pi}{2}, \phi _{i} = \psi _{i} = 0 \right)$ , and (θ_i = 0, ψ_i + φ_i = 0) we can do similar kind of analysis from Case–I and Case–II respectively.

We shall adopt the notation and terminology of [23]. As a consequence, the distance between the primaries does not change and is taken equal to one; the sum of masses of the primaries is also taken one. The unit of time is chosen so as to make the gravitational constant unity. Besides this the principle axes of the primaries are oriented to the synodic axes by Euler’s angels (θ_i, ψ_i, φ_i(i = 1, 2)). Since the axes are supposed to rotate with the same angular velocity as that of the rigid bodies and the bodies are moving around their center of mass without rotation, the Euler’s angles remain constant throughout the motion. Using dimensionless variables, the equations of motion of the infinitesimal mass m₃ in a synodic coordinate system (x, y) are

x¨−2ny˙=∂Ω∂x,y¨+2nx˙=∂Ω∂y.$$ \begin{equation} \begin{array}{rl} \ddot{x}- 2n\dot{y} = \dfrac{\partial \Omega}{\partial x},\\ \ddot{y}+ 2n\dot{x} = \dfrac{\partial \Omega}{\partial y}. \end{array} \end{equation} $$(1)

where,

Ω=n22[(1−μ)r12+μr22]+(1−μ)r1+μr2+(1−μ)2m1r13[I1+I2+I3−3I]+μ2m2r23[I′1+I′2+I′3−3I′],$$ \begin{equation} \begin{array}{rl} \Omega = &\dfrac{n^{2}}{2} \left[\left(1-\mu \right)r_{1}^{2} +\mu r_{2}^{2} \right]+\dfrac{\left(1-\mu \right)}{r_{1}}\\ &+\dfrac{\mu}{r_{2}} +\dfrac{\left(1-\mu \right)}{2m_{1} r_{1}^{3}} \left[I_{1} + I_{2} + I_{3} -3I \right]\\ &+\dfrac{\mu}{2m_{2} r_{2}^{3}} \left[I'_{1} +I'_{2} +I'_{3} -3I'\right], \end{array} \end{equation} $$(2)r12=(x−μ)2+y2,r22=(x+1−μ)2+y2.$$ \begin{equation} \begin{array}{rl} r_{1}^{2} = (x - \mu)^{2} + y^{2},\\ r_{2}^{2} = (x + 1 - \mu)^{2} + y^{2}. \end{array} \end{equation} $$(3)

Here μ is the ratio of mass of the smaller primary to the total mass of primaries and 0≤μ≤12 $0\le \mu \le \dfrac{1}{2}$ , i.e., μ = m2m1+m2≤12 $\dfrac{m_{2}}{m_{1} +m_{2}} \le \dfrac{1}{2}$ with m₁ ≥ m₂ being the masses of the primaries. I₁, I₂, I₃ are the principal moments of inertia of the triaxial rigid body of mass m₁ at its center of mass, with a, b, c as its axes. I is the moment of inertia about a line joining the center of the rigid body of mass m₁ and the infinitesimal body of mass m₃ and is given by

I=I1l′12+I2m′12+I3n′12$$ \begin{equation} I=I_{1} l _{1}^{'2} + I_{2} m_{1}^{'2} + I_{3} n_{1}^{'2} \end{equation} $$(4)

where l′₁, m′₁ and n′₁ are the directional cosines of the line respect to its principal axes. I′₁, I′₂, I′₃ are the principal moments of inertia of the triaxial rigid body of mass m₂ at its center of mass, with a′, b′, c′ as its axes. I′ is the moment of inertia about a line joining the center of the rigid body of mass m₂ and the infinitesimal body of mass m₃ and is given by

I′=I′1l′22+I′2m′22+I′3n′22$$ \begin{equation} I^{'} =I_{1}^{'} l _{2}^{'2} + I_{2}^{'} m_{2}^{'2} + I_{3}^{'} n_{2}^{'2} \end{equation} $$(5)

where l′₂, m′₂ and n′₂ are the directional cosines of the line respect to its principal axes. We denote the unit vectors along the principle axes at p₁(or p₂) by i, j, k and the unit vectors parallel to the synodic axes by I, J, K with the help of Euler’s angles (θ_i, ψ_i, φ_i), (i = 1, 2).

I=a1ii^+b1ij^+c1ik^,J=a2ii^+b2ij^+c2ik^,K=a3ii^+b3ij^+c3ik^,$$ \begin{equation} \begin{array}{rl} & I=a_{1i}\hat{i}+ b_{1i}\hat{j} + c_{1i}\hat{k},\\ & J=a_{2i} \hat{i}+ b_{2i}\hat{j} + c_{2i}\hat{k},\\ & K=a_{3i} \hat{i}+ b_{3i}\hat{j} + c_{3i}\hat{k}, \end{array} \end{equation} $$(6)

(i = 1, 2), where

a1i=−sinϕisinψi+cosθicosϕicosψi,a2i=cosϕisinψi+cosθisinϕicosψi,a3i=−sinθicosψi,b1i=−sinϕicosψi−cosθicosϕisinψi,b2i=cosϕicosψi−cosθisinϕisinψi,b3i=sinθisinψi,c1i=sinθicosϕi,c2i=sinθisinϕi,c3i=cosθi,$$ \begin{equation} \begin{array}{rl} & a_{1i} = -\sin\phi _{i} \sin\psi _{i} + \cos\theta _{i} \cos\phi _{i} \cos\psi _{i},\\ & a_{2i} = \cos\phi _{i} \sin \psi _{i} + \cos\theta _{i} \sin\phi _{i} \cos\psi _{i},\\ & a_{3i} = -\sin\theta _{i} \cos\psi _{i}, \\ & b_{1i} = -\sin \phi _{i} \cos \psi _{i} -\cos \theta _{i} \cos\phi_{i} \sin\psi _{i}, \\ & b_{2i} = \cos\phi _{i} \cos\psi _{i} -\cos\theta _{i} \sin\phi _{i} \sin\psi _{i}, \\ & b_{3i} = \sin\theta _{i}\sin\psi _{i}, \\ & c_{1i} = \sin\theta _{i}\cos\phi _{i},\\ & c_{2i} = \sin\theta _{i}\sin\phi _{i}, \\ & c_{3i} =\cos\theta _{i}, \end{array} \end{equation} $$(7)

(i = 1, 2).

The axes O(xyz) have been defined by [23]. Now, Ω in equation (2) can be written as

I=a1ii^+b1ij^+c1ik^,J=a2ii^+b2ij^+c2ik^,K=a3ii^+b3ij^+c3ik^,$$ \begin{equation} \begin{array}{rl} & I=a_{1i}\hat{i}+ b_{1i}\hat{j} + c_{1i}\hat{k},\\ & J=a_{2i} \hat{i}+ b_{2i}\hat{j} + c_{2i}\hat{k},\\ & K=a_{3i} \hat{i}+ b_{3i}\hat{j} + c_{3i}\hat{k}, \end{array} \end{equation} $$(8)

Ω=n221−μr12+μr22+1−μr1+μr2+1−μ2r132A1+A2+A3−31r12A2+A3a11x−μ+a21y2+A1+A3b11x−μ+b21y2+A2+A1c11x−μ+c21y2+μ2r232A1′+A2′+A3′−31r22A2′+A3′a12x+1−μ+a22y2+A1′+A3′b12x+1−μ+b22y2+A2′+A1′c12x+1−μ+c22y2,$$ \begin{equation} \begin{array}{rl} \Omega =& \dfrac{n^{2}}{2} \left[\left(1-\mu \right)r_{1}^{2} +\mu r_{2}^{2} \right] +\dfrac{\left(1-\mu \right)}{r_{1} } +\dfrac{\mu }{r_{2} } \\ &+\dfrac{\left(1-\mu \right)}{2r_{1}^{3} } \left[2\left(A_{1} +A_{2} +A_{3} \right)-3\dfrac{1}{r_{1}^{2} } \left\{\begin{array}{l} {\left(A_{2} +A_{3} \right)\left(a_{11} \left(x-\mu \right)+a_{21} y\right)^{2} }\\ {+\left(A_{1} +A_{3} \right)\left(b_{11} \left(x-\mu \right)+b_{21} y\right)^{2} }\\ {+\left(A_{2} +A_{1} \right)\left(c_{11} \left(x-\mu \right)+c_{21} y\right)^{2} }\\ \end{array}\right\}\right]\\ &+\dfrac{\mu }{2r_{2}^{3} } \left[2\left(A'_{1} +A'_{2} +A'_{3} \right)-3\dfrac{1}{r_{2}^{2} } \left\{\begin{array}{l} {\left(A'_{2} +A'_{3} \right)\left(a_{12} \left(x+1-\mu \right)+a_{22} y\right)^{2} }\\ {+\left(A'_{1} +A'_{3} \right)\left(b_{12} \left(x+1-\mu \right)+b_{22} y\right)^{2} }\\ {+\left(A'_{2} +A'_{1} \right)\left(c_{12} \left(x+1-\mu \right)+c_{22} y\right)^{2} } \end{array}\right\}\right], \end{array} \end{equation} $$(9)

where

A1=a25R2,A2=b25R2,A3=c25R2,$$ \begin{equation} A_{1} =\dfrac{a^{2}}{5R^{2}}\, , \quad A_{2} =\dfrac{b^{2}}{5R^{2}}\,, \quad A_{3} =\dfrac{c^{2}}{5R^{2}}, \end{equation} $$(10)A′1=a′25R2,A′2=b′25R2,A′3=c′25R2,$$ \begin{equation} \begin{array}{rl} A_{1}^{'} =\dfrac{a^{'2}}{5R^{2}}\, , \quad A_{2}^{'} =\dfrac{b^{'2}}{5R^{2}}\,, \quad A_{3}^{'} =\dfrac{c^{'2}}{5R^{2}}, \end{array} \end{equation} $$(11)

and R is the distance between the primaries. The mean motion, n is given by,

n2=1+32[2(A1+A2+A3)−3a112(A2+A3)−3b112(A1+A3)−3c112(A2+A1)]+32[2(A′1+A′2+A′3)−3a122(A′2+A′3)−3b122(A′1+A′3)−3c122(A′2+A′1)].$$ \begin{equation} \begin{array}{rl} n^{2} =&1+\dfrac{3}{2} \left[2\left(A_{1} +A_{2} +A_{3} \right)-3a_{11}^{2} \left(A_{2} +A_{3} \right)-3b_{11}^{2} \left(A_{1} +A_{3} \right)-3c_{11}^{2} \left(A_{2} +A_{1} \right)\right] \\ &+\dfrac{3}{2} \left[2\left(A'_{1} +A'_{2} +A'_{3} \right)-3a_{12}^{2} \left(A'_{2} +A'_{3} \right)-3b_{12}^{2} \left(A'_{1} +A'_{3} \right)-3c_{12}^{2} \left(A'_{2} +A'_{1} \right)\right]. \end{array} \end{equation} $$(12)

Equation (1) permit an integral analogous to Jacobi integral

x˙2+y˙2−2Ω+C=0,f(x,y,x˙,y˙)=x˙2+y˙2−2Ω+C=0.$$ \begin{equation} \begin{array}{rl} &\dot{x}^{2}+\dot{y}^{2}-2\Omega + C =0, \\ & f \left(x, y,\dot{x}, \dot{y} \right) = \dot{x}^{2} +\dot{y}^{2}- 2\Omega + C = 0. \end{array} \end{equation} $$(13)

The liberation points are the singularities of the manifold Therefore, these points are the solutions of the equations Ω_x = 0, Ω_y = 0 . We have Ω_x and Ω_y are established by [11, 12].

The particular solutions may be obtained with any desired degree of numerical accuracy. The next question is the examination of orbits in the vicinity of these particular solutions. In general, let X_o, Y_o be the coordinates correspond to any one of the particular solutions. They satisfy equations (1).

X¨0−2nY˙0=(∂Ω∂X)0,Y¨0+2nX˙=(∂Ω∂Y)0.$$ \begin{equation} \begin{array}{rl} &\ddot X_{0} - 2n\dot Y_{0} =\left(\dfrac{\partial \Omega}{\partial X} \right)_{0}, \\ & \ddot Y_{0} + 2n \dot X_{} =\left(\dfrac{\partial \Omega}{\partial Y} \right)_{0}. \end{array} \end{equation} $$(28)

In which the subscript o indicates that the particular derivatives of Ω must be evaluated for X = X_o, Y = Y_o. If now X = X₀ + ξ, Y = Y₀ + η are substituted into the equations, there results,

ξ¨−2nη˙=Ωxx0ξ+Ωxy0η,η¨+2nξ˙=Ωyx0ξ+Ωyy0η.$$ \begin{equation} \begin{array}{rl} & \ddot{\xi} - 2n \dot{\eta} = \Omega^0 _{xx} \xi +\Omega^0 _{xy} \eta, \\ & \ddot{\eta} + 2n \dot{\xi} = \Omega^0 _{yx}\xi +\Omega^0 _{yy} \eta. \end{array} \end{equation} $$(29)

The super script of second partial derivatives of Ω refers to their values at X = X₀ and Y = Y₀. Specific numerical values may be obtained for any of the particular solutions. Let a solution of equations (29) be

ξ=Aexp(λt),η=Bexp(λt).$$ \begin{equation} \begin{array}{rl} &\xi =A \exp \left(\lambda t\right),\\ &\eta = B \exp \left(\lambda t\right). \end{array} \end{equation} $$(30)

Substituting (30) into (29) gives for the equations that must be satisfied by the coefficients A and B

λ4+(4n2−Ωxx0−Ωyy0)λ2+Ωxx0Ωyy0−(Ωxy0)2=0$$ \begin{equation} \lambda ^{4} +\left(4n^{2} -\Omega^0_{xx} -\Omega^0_{yy} \right)\lambda ^{2} +\Omega^0_{xx} \Omega^0_{yy} -(\Omega^0_{xy})^{2} =0 \end{equation} $$(31)

The character of the solution of the differential equations depends on the character of the solution for λ² from this quadratic equation. The solution is stable only if the quadratic has two unequal negative roots for λ² . Now, we consider two cases of stationary rotational motion of the primaries.

To obtain numerical solution of the given system for equations of motion of infinitesimal mass, we will use Runge–Kutta–Gill method with constant interval. This method is correct of order 4. Before applying method, we are shifting location of both primaries towards positive x–axis. So, new location of bigger primary is at(x, y) = (−μ, 0) and that of smaller primary is at (1−μ, 0).

Since the standard Runge-Kutta–Gill fourth order method for first order differential equation ẏ = f(x, y) is as follows.

yn+1=yn+16[k1+(2−2)k2+(2+2)k3+k4].$$ \begin{equation} y_{n+1} = y_{n} + \dfrac{1}{6}\left[k_{1} + (2 -\sqrt{2})k_{2}+ (2 +\sqrt{2})k_{3} + k_{4}\right]. \end{equation} $$(37)

where

k1=hf(xn,yn),k2=hf(xn+h2,yn+k12),k3=hf[xn+h2,yn+k12(−1+2)+(1−122)k2],k4=hf[xn+h,yn−k222+(1+122)k3].$$ \begin{equation} \begin{array}{l} k_{1} = hf(x_{n}, y_{n}),\\ k_{2} = hf(x_{n}+\dfrac{h}{2}, y_{n}+\dfrac{k_{1}}{2}),\\ k_{3} = hf\left[x_{n}+\dfrac{h}{2}, y_{n}+ \dfrac{k_{1}}{2}(-1+\sqrt{2})+(1-\dfrac{1}{2}\sqrt{2})k_{2}\right],\\ k_{4} = hf\left[x_{n}+h, y_{n}- \dfrac{k_{2}}{2}\sqrt{2}+(1+\dfrac{1}{2}\sqrt{2})k_{3}\right]. \end{array} \end{equation} $$(38)

Here, h is the step size. The system of equations (1) is second order ordinary differential equations. To convert equations (1) in to system of first order differential equations, consider, y₁ = x, y₂ = ẋ, y₃ = y and y₄ = ẏ. Thus, equations (1) can be written as,

y1˙=y2,y2˙=2ny4+n2y1−(1−μ)(y1+μ)r13−μ(y1+μ−1)1r23−1.5(1−μ)(y1+μ)r15((2A2+4A3−A1)−5r12(A3((y1+μ)2)+A2y32))−1.5μ(y1+μ−1)r25((2A2′+4A3′−A1′)−5r22(A3′(y1−1+μ)2+A2′y32)),y3˙=y4,y4˙=−2ny2+n2y3−(1−μ)y3r13−μy31r23−1.5(1−μ)y3r15((4A2+2A3−A1)−(5r12)(A3((y1+μ)2)+A2y32))−1.5μy3r24((4A2′+2A3′−A1′)−(5r22(A3′(y1−1+μ)2+A2′y32))).$$ \begin{equation} \begin{array}{rl} &\dot{y_{1}} = y_{2},\\ &\dot{y_{2}} = 2ny_{4} + n^{2}y_{1} - \dfrac{(1-\mu)(y_{1} + \mu)}{r_{1}^{3}} - \mu (y_{1} + \mu -1)\left[\dfrac{1}{r_{2}^{3}}\right]\\ &\qquad -\dfrac{1.5(1-\mu)(y_{1}+\mu)}{r_{1}^{5}}((2A_2+4A_3- A_1)- \dfrac{5} {r_{1}^{2}} (A_3((y_{1}+ \mu)^2)+ A_2 y_{3}^2))\\ &\qquad-\dfrac{1.5 \mu (y_{1}+\mu-1)}{r_{2}^{5}} ((2{A_2}^{'}+4{A_3}^{'}-{A_1}^{'}) - \dfrac{5} {r_{2}^{2}} ({A_3}^{'} (y_{1}-1+\mu)^2 + {A_2}^{'}y_{3}^2)),\\ &\dot{y_{3}} = y_{4},\\ &\dot{y_{4}} = - 2ny_{2} + n^{2}y_{3} -\dfrac{(1-\mu)y_{3}}{r_{1}^{3}} - \mu y_{3}\left[\dfrac{1}{r_{2}^{3}}\right]\\ &\qquad-\dfrac{1.5(1-\mu)y_{3}}{r_{1}^{5}}((4A_2+2A_3-A_1)-(\dfrac{5}{r_{1}^{2}}) (A_3((y_{1}+\mu)^2)+A_2 y_{3}^2))\\ &\qquad-\dfrac{1.5 \mu y_{3}} {r_{2}^{4}} ((4{A_2}^{'}+2{A_3}^{'}-{A_1}^{'}) -(\dfrac{5}{r_{2}^{2}} ({A_3}^{'}(y_{1}-1+\mu)^2+{A_2}^{'} y_{3}^2))). \end{array} \end{equation} $$(39)

We use Runge–Kutta–Gill fourth order method to integrate the system of first order differential equations. The algorithm for this method is as follows:

The system of equations are given by

y1.=f1(y1,y2,y3,y4),y2.=f2(y1,y2,y3,y4),y3.=f3(y1,y2,y3,y4),y4.=f4(y1,y2,y3,y4).$$ \begin{eqnarray*}\dot{y_{1}} = f_{1}(y_{1}, y_{2}, y_{3}, y_{4}), \\ \dot{y_{2}} = f_{2}(y_{1}, y_{2}, y_{3}, y_{4}),\\ \dot{y_{3}} = f_{3}(y_{1}, y_{2}, y_{3}, y_{4}),\\ \dot{y_{4}} = f_{4}(y_{1}, y_{2}, y_{3}, y_{4}). \end{eqnarray*} $$