On the perturbed restricted three-body problem

The three-body problem has a special relevance, particularly in astrophysics and astrodynamics. In general, the three-body problem is classified into two types:

Thus, the general problem has some applications in celestial mechanics such as the dynamics of triple star systems (and only a very few in space dynamics and solar system dynamics), whereas the restricted problem plays an important role in studying the motion of artificial satellites. It can be used also to evaluate the motion of the planets, minor planets and comets. The restricted problem gives an accurate description not only regarding the motion of the Moon but also with respect to the motion of other natural satellites. Furthermore, the restricted problem has many applications not only in celestial mechanics but also in physics, mathematics and quantum mechanics, to name a few. In quantum mechanics, a general form of the restricted problem is formed to solve the Schrödinger equation of helium-like ions. Furthermore, in modern solid state physics, the restricted problem can be used to discuss the motion of an infinitesimal mass affected not only by the gravitational field but also by light pressure from one (or both) of the primaries, which is called the photogravitational problem.

In the last decades, the restricted three-body problem has been enhanced by a great number of researches. Many of them deal with the effects of the perturbed forces such as the lack of sphericity, photogravitational force, Coriolis and centrifugal forces, variation of masses, the Pointing-Robertson and Yarkovsky effects and the atmospheric drag and solar wind. The Kirkwood gaps in the ring of the asteroid’s orbits lying between the orbits of Mars and Jupiter are examples of the perturbation produced by Jupiter on an asteroid. This motivates many researchers to study the restricted three-body problem under the effects of small perturbations in Coriolis and centrifugal forces when the three involved objects are oblate, . . . , etc.

The significance of the restricted problem when describing actual physical situations can be judged by the results obtained when these are compared with observations. It is worth mentioning that the utility might be prejudged by order of magnitude evaluations regarding the masses and the distances of the participating bodies. In this way, a classic example in space dynamics is the Sun-Earth-Moon system.

Some important contributions related to the libration points in the restricted three-body problem with one or both primaries are oblate spheroids when the equatorial plane is coincident with the plane of motion were studied by Subbarao and Sharma [28], Sharma and Subbarao [23], and Markellos et al. [20]. Abouelmagd [9] also studied the effects of oblateness J₂ and J₄ for the more massive primary in the planar restricted three-body problem on the locations of the triangular points and their linear stability. He found that these locations are affected by the coefficients of oblateness. Furthermore, he investigated that the triangular points are stable for 0 < μ < μ_c and unstable when μ_c ≤ μ ≤ 1/2, where μ_cis the critical mass parameter which depends on the coefficients of oblateness.

The existence of libration points and their linear stability as well as periodic orbits around these points when the more massive primary is radiating and the smaller is an oblate spheroid were studied by Abouelmagd and Sharaf [3]. Their study also includes the effects of oblateness of J¯2i:i=1,2$\begin{array}{} \displaystyle \overline{J}_{2i}:i=1,2 \end{array}$, with respect to the smaller primary in the restricted three-body problem.

The restricted problem when the three participating bodies are oblate spheroids was also studied by Elipe and Ferrer [15], and El-Shaboury and El-Tantawy [18]. When one or two of the primaries are triaxial bodies this problem was introduced by El-Shaboury et al. [17], Khanna and Bhatnagar [19], and Sharma et al. [24].

Various researchers made studies in the restricted three-body problem under the effects of small perturbations in centrifugal and Coriolis forces such as in Szebehely [29], Bhatnagar and Hallan [12], Devi and Singh [14], and Shu and Lu [25], to quote some of them.

The effect of small perturbations ε,ε′ in Coriolis and centrifugal forces with variable mass in the restricted three-body problem has been studied by Singh [26]. He found that in the nonlinear sense the triangular points are stable for all mass ratios in the range of linear stability except for three mass ratios which depend on ε,ε′ and the constant β due to the variation of mass governed by Jeans’ law.

Mittal et al. [21] have studied periodic orbits generated by Lagrangian solutions of the restricted three-body problem when the bigger body is a source of radiation and the smaller is an oblate spheroid. They used the definition of Karimov and Sokolsky for mobile coordinates to determine these orbits, they also used the predictor method to draw them.

Singh and Begha [27] have studied the existence of periodic orbits around the triangular points in the restricted three-body problem when the bigger primary is a triaxial and the smaller primary is considered as an oblate spheroid in the range of linear stability with the perturbed Coriolis and centrifugal forces. They deduced that long and short periodic orbits exist around these points and their periods, orientation and eccentricities are affected by the non sphericity and the perturbation in Coriolis and centrifugal forces.

The existence of libration points and their linear stability when the three participating bodies are oblate spheroids and the primaries are radiation source as well were studied by Abouelmagd and El-Shaboury [2]. They found that the collinear points are still unstable while the triangular points are stable for 0 < μ < μ_c, and unstable for μ_c ≤ μ ≤ 1/2, where μ_c ∊ (0,1/2). They also deduced that for these points the range of stability will decrease. In addition they studied the periodic orbits around these points in the range 0 < μ < μ_c.

From a physical point of view, it is unreasonable to consider all objects as being point masses with no physical dimensions. This is conflict with the real cases for the celestial bodies. On the other hand the effect of rotation causes deformation in the shape of the objects at the poles as might be expected. For this reason most objects may be treated to a good approximation as oblate spheroids, in this context there are a number of contributions established by: Abouelmagd [5, 6], Abouelmagd et al. [7, 8], and Abouelmagd et al. [9–11]. Therefore in this paper we will generalize the restricted three-body, in which the three participating bodies are oblate spheroids as well as the existence of small perturbations in the Coriolis and centrifugal forces. The existence of libration points and their linear stability and the periodic orbits around these points will be studied. This model could be applicable in astrodynamics or astrophysics.

Now we assume that the three participating bodies are oblate spheroids. Therefore we refer to oblateness parameters of the bigger, smaller primaries and infinitesimal body by A₁,A₂, and A, respectively. We also denote the small perturbations in Coriolis and centrifugal forces by ε_i : i = 1,2, where 0 < A₁ ≪ 1, 0 < A₂ ≪ 1, 0 < A ≪ 1, and ε_i ≪ 1.

It is well known that there are five exact solutions in three-body problem which the three masses maintain a constant configuration which revolves with constant angular velocity. An important specialization of the three-body problem is the restricted three-body problem in which m is infinitesimal body, and m₁ and m₂ are the bigger and smaller primaries, respectively. They move in circular orbits around their barycenter. The smallness of m means that it does not influence the motion of m₁ and m₂. For many purposes, it is convenient to describe the motion of m in a coordinate system which is attached to both m₁ and m₂. In this rotating coordinate system, the five Lagrangian solutions show up as five fixed points at which m would be stationary if placed there with zero velocity.

If the object has axial symmetry, it can be shown from potential theory that the gravitational potential experienced by the satellite can be written as in (Murray and Dermott, [22]), namely,

V=−Gm0r[1−∑n=2∞Jn(R0r)npn(sinδ)],$$\begin{array}{} \displaystyle V=-\dfrac{Gm_0}{r}\Bigg[1-\sum_{n=2}^{\infty} J_n \Bigg(\dfrac{R_0}{r}\Bigg)^n p_n(\sin \delta)\Bigg], \end{array}$$(1)

where

If the motion of the satellite and the object are in the same plane (δ = 0), then Eq. (1) becomes

V=−Gm0r[1−∑n=2∞Jn(R0r)npn(0)],$$\begin{array}{} \displaystyle V=-\dfrac{Gm_0}{r}\Bigg[1-\sum_{n=2}^{\infty} J_n \Bigg(\dfrac{R_0}{r}\Bigg)^n p_n(0)\Bigg], \end{array}$$(2)

where p2n(0)=(−1)n(2n)!22n(n!)2,p2n+1(0)=0$\begin{array}{} \displaystyle p_{2n}(0)=\dfrac{(-1)^{n}(2n)!}{2^{2n}(n!)^2}, p_{2n+1}(0)=0 \end{array}$. In addition, let us consider the first zonal harmonic. Hence, Eq. (6) can be reduced to

V=−Gm0r[1r+J2R022r3].$$\begin{array}{} \displaystyle V=-\dfrac{Gm_0}{r}\Bigg[\dfrac{1}{r}+\dfrac{J_2 R_0^2}{2r^3}\Bigg]. \end{array}$$(3)

Now we assume that m₁, m₂ and m are the masses of the bigger, smaller and infinitesimal bodies, r₁ and r₂ are the distances from m₁ and m₂ to m, respectively.

Let m₁ and m₂ have a circular orbit about their common center of mass and m is moving under their gravitational field in the same plane but it is so small that it does not affect their motion. Furthermore let:

Using Eq. (3), then the potential V₁₂ between m₁ and m₂ can be written as

V12=−Gm1m2[1R+A1+A22R3].$$\begin{array}{} \displaystyle V_{12}=-Gm_1m_2\Bigg[\dfrac{1}{R}+\dfrac{A_1+A_2}{2R^3}\Bigg]. \end{array}$$(4)

Thus, the equation of motion for m₂ is

R¨_=−(m1+m2m1m2)∇V12,$$\begin{array}{} \displaystyle \underline{\ddot{R}}=-\Bigg(\dfrac{m_1+m_2}{m_1m_2}\Bigg)\nabla V_{12}, \end{array}$$(5)

where

∇=i∂∂u+j∂∂ν.$$\begin{array}{} \displaystyle \nabla=i\dfrac{\partial}{\partial u}+j\dfrac{\partial}{\partial \nu}. \end{array}$$

Hence Eq. (5) has the following particular solutions:

Therefore, the perturbed mean motion is given by the following relation

n2=G(m1+m2)[1R3+3(A1+A2)2R5].$$\begin{array}{} \displaystyle n^2=G(m_1+m_2)\Bigg[\dfrac{1}{R^3}+\dfrac{3(A_1+A_2)}{2R^5}\Bigg]. \end{array}$$(6)

Now we assume that XY is the sidereal frame while xy are the synodic frames such that they have the same origin at the center of mass for the primaries. We also chose the synodic frame to rotate with angular velocity n in positive direction (see Fig. (1)). Let the coordinates of m₁,m₂, and m in the sidereal frame be (X₁, Y₁), (X₂, Y₂), and (X, Y ), respectively, but in the synodic frame are (x₁, y₁),(x₂, y₂), and (x, y), respectively, too.

Fig. 1

Based on the above discussion, the sidereal coordinates are related to the synodic coordinates via the following relation:

X=x cos nt−y sin nt.$$\begin{array}{} \displaystyle X=x\cos nt-y\sin nt. \end{array}$$(7)

Y=xsinnt+ycosnt.$$\begin{array}{} \displaystyle Y=x\sin nt+y\cos nt. \end{array}$$(8)

Using Newton’s law, the equations of motion for infinitesimal body in a sidereal frame are as follows:

mX¨=∂V∂X.$$\begin{array}{} \displaystyle m\ddot{X}=\dfrac{\partial V}{\partial X}. \end{array}$$(9)

mY¨=∂V∂Y.$$\begin{array}{} \displaystyle m\ddot{Y}=\dfrac{\partial V}{\partial Y}. \end{array}$$(10)

From Eq. (3), the potential V₁ between m and m₁ is given by

V1=−Gmm1[1r1+A1+A2r13],$$\begin{array}{} \displaystyle V_{1}=-Gmm_1\Bigg[\dfrac{1}{r_1}+\dfrac{A_1+A}{2r_1^3}\Bigg], \end{array}$$(11)

and the potential V between m, m₁ and m₂ is also given by

V2=−Gmm2[1r2+A2+A2r23],$$\begin{array}{} \displaystyle V_{2}=-Gmm_2\Bigg[\dfrac{1}{r_2}+\dfrac{A_2+A}{2r_2^3}\Bigg], \end{array}$$(12)

while the total potential V between m, m₁, and m₂ is controlled by

V=−Gmm1[1r1+A1+A2r13]−Gmm2[1r2+A2+A2r23],$$\begin{array}{} \displaystyle V=-Gmm_1\Bigg[\dfrac{1}{r_1}+\dfrac{A_1+A}{2r_1^3}\Bigg]-Gmm_2\Bigg[\dfrac{1}{r_2}+\dfrac{A_2+A}{2r_2^3}\Bigg], \end{array}$$(13)

once the notation and terminology of Szebehely [30] has been adopted. As a consequence, the distance between the primaries does not change and is equal to one, the sum of masses of the primaries is also taken as one, the unit time is chosen as to make the unperturbed mean motion and the gravitational constant unity, then it follows that m₁ = 1 − μ, whose coordinate is (μ, 0), m₂ = μ ≤ 1/2 and located at (μ − 1,0).

Substituting Eqs. (7-8) and (11-13) into (9-138) that gives the equations of motion in a synodic coordinate system (x − y) for infinitesimal body with dimensionless variables as

where

U=12n2[(1−μ)r12+μr22]+(1−μ)[1r1+A1+A2r13]+μ[1r2+A2+A2r23].$$\begin{array}{} \displaystyle U=\dfrac{1}{2}n^2\Bigg[(1-\mu)r_1^2+\mu r_2^2\Bigg]+(1-\mu)\Bigg[\dfrac{1}{r_1}+\dfrac{A_1+A}{2r_1^3}\Bigg]+\mu\Bigg[\dfrac{1}{r_2}+\dfrac{A_2+A}{2r_2^3}\Bigg]. \end{array}$$(14)

n2=1+32(A1+A2).$$\begin{array}{} \displaystyle n^2=1+\dfrac{3}{2}(A_1+A_2). \end{array}$$(15)

r12=(x−μ)2+y2.$$\begin{array}{} \displaystyle r_1^2=(x-\mu)^2+y^2. \end{array}$$(16)

r22=(x−μ+1)2+y2.$$\begin{array}{} \displaystyle r_2^2=(x-\mu+1)^2+y^2. \end{array}$$(17)

Eqs. (i)-(ii) and (14) represent the equation of motion for the infinitesimal body in a rotating coordinate system when the three participating bodies are oblate spheroids.

In the previous equations the perturbed effects of Coriolis and centrifugal forces are switched off. In physics the Coriolis effect is a deflection of moving objects when they are viewed in a rotating reference frame. In a reference frame with clockwise rotation the deflection is to the left of the motion of the object, in one with counter-clockwise rotation the deflection is to the right. The Newton’s laws of motion govern the motion of an object in a non-accelerating inertial frame of reference. When Newton’s laws are transformed to a uniformly rotating frame of reference, the Coriolis and centrifugal forces appear.

Both forces are proportional to the mass of the object. The Coriolis force is proportional to the rotation rate and the centrifugal force is proportional to its square. The Coriolis force acts in a direction perpendicular to the rotation axis and to the velocity of the body in the rotating frame and is proportional to the object’s speed in the rotating frame. The centrifugal force acts outwards in the radial direction and is proportional to the distance of the body from the axis of the rotating frame.

It is important to note that these forces do not arise from any physical agency: they arise solely as a result of rotation for the coordinate system. Thus, if the angular velocity is reduced to zero, both the Coriolis and the centrifugal forces will vanish. For this reason, the Coriolis and centrifugal forces are referred to as “fictitious” or “inertial” forces. These additional forces allow the application of Newton’s laws to a rotating system. They are correcting factors that do not exist in a non-accelerating or inertial reference frame.

Anyway, if the perturbed effects of Coriolis and centrifugal forces are considered, let φ, ψ denote these forces, respectively. Assume that ε₁ represents the perturbation in Coriolis force, while ε₂ is the accompanying perturbation in the centrifugal force such that

φ=1+ε1,ε1≪1.$$\begin{array}{} \displaystyle \varphi=1+\varepsilon_1,\qquad \varepsilon_1\ll 1. \end{array}$$(18)

ψ=1+ε2,ε2≪1.$$\begin{array}{} \displaystyle \psi=1+\varepsilon_2,\qquad \varepsilon_2\ll 1. \end{array}$$(19)

Hence, Eqs. (i)-(ii) and (14) can be written as follows:

x¨−2φny˙=Ωx.$$\begin{array}{} \displaystyle \ddot{x}-2\varphi n\dot{y}=\Omega_x. \end{array}$$(20)

y¨+2φnx˙=Ωy,$$\begin{array}{} \displaystyle \ddot{y}+2\varphi n\dot{x}=\Omega_y, \end{array}$$(21)

where

Ω=12n2ψ[(1−μ)r12+μr22]+(1−μ)[1r1+A1+A2r13]+μ[1r2+A2+A2r23].$$\begin{array}{} \displaystyle \Omega=\dfrac{1}{2}n^2\psi\Bigg[(1-\mu)r_1^2+\mu r_2^2\Bigg]+(1-\mu)\Bigg[\dfrac{1}{r_1}+\dfrac{A_1+A}{2r_1^3}\Bigg]+\mu\Bigg[\dfrac{1}{r_2}+\dfrac{A_2+A}{2r_2^3}\Bigg]. \end{array}$$(22)

It is worth mentioning that the subscripts on the right handside denote the partial derivation with respect to x and y, respectively.

The libration points are singular points of the differential equations of motions in the restricted problem of three bodies. They are also equilibrium points since the gravitational force on a mass placed in such a point are balanced by the centrifugal force. The three libration points are called the collinear points which are found on the line connecting with the primaries. The other two are called the triangular points which are symmetrical with respect to this line and form triangles, see Fig. (1).

From Eqs. (i)-(91), we can obtain the Jacobi integral as

x¨2+y¨2−2Ω+c=0,$$\begin{array}{} \displaystyle \ddot{x}^2+\ddot{y}^2-2\Omega+c=0, \end{array}$$(23)

where c is a constant of integration.

The locations of these points will be determined by solving the following expressions:

Ωx=Ωx˙=Ωy=Ωy˙=0,$$\begin{array}{} \displaystyle \Omega_x=\Omega_{\dot{x}}=\Omega_y=\Omega_{\dot{y}}=0, \end{array}$$(24)

where the first partial derivatives of Ω will be governed by

Ωx=(x−μ)f1(r1)+(x−μ+1)f2(r2).$$\begin{array}{} \displaystyle \Omega_x=(x-\mu)f_1(r_1)+(x-\mu+1)f_2(r_2). \end{array}$$(25)

Ωy=y[f1(r1)+f2(r2)].$$\begin{array}{} \displaystyle \Omega_y=y[f_1(r_1)+f_2(r_2)]. \end{array}$$(26)

f1(r1)=(1−μ)[ψn2−(1r13+3(A1+A)2r15)].$$\begin{array}{} \displaystyle f_1(r_1)=(1-\mu)\Bigg[\psi n^2-\Bigg(\dfrac{1}{r_1^3}+\dfrac{3(A_1+A)}{2r_1^5}\Bigg)\Bigg]. \end{array}$$(27)

f2(r2)=μ[ψn2−(1r23+3(A2+A)2r25)].$$\begin{array}{} \displaystyle f_2(r_2)=\mu\Bigg[\psi n^2-\Bigg(\dfrac{1}{r_2^3}+\dfrac{3(A_2+A)}{2r_2^5}\Bigg)\Bigg]. \end{array}$$(28)

From Eqs. (25-26), two cases yield. They will be investigated in the following subsection.

3.1

The collinear points

The collinear points L_i : i = 1, 2, 3 are the solutions of Eqs. (24) and (25-26) when y = 0. Therefore,

ν=(−1)ir24(r15−pr12−a)r14(r25−pr22−b):i=1,2,3,$$\begin{array}{} \displaystyle \nu=(-1)^{i}\dfrac{r_2^4(r_1^5-pr_1^2-a)}{r_1^4(r_2^5-pr_2^2-b)}:\qquad i=1,2,3, \end{array}$$(29)

where i indicates the value of ν at L_i : i = 1, 2, 3. Further,

ν=μ1−μ,p=1q,q=ψn2,$$\begin{array}{} \displaystyle \nu=\dfrac{\mu}{1-\mu},\quad p=\dfrac{1}{q},\quad q=\psi n^2, \end{array}$$(30)

a=3(A1+A)2q,b=3(A2+A)2q.$$\begin{array}{} \displaystyle a=\dfrac{3(A_1+A)}{2q},\quad b=\dfrac{3(A_2+A)}{2q}. \end{array}$$(31)

Lagrange developed a useful method for inverting series expansions which could be applicable in the present work. He investigated that if a variable χ can be expressed as a function of 𝒴 as in the below relation (Murray and Dermott, [22]):

Y=χ+ϑϕ(Υ),ϑ<1.$$\begin{array}{} \displaystyle \mathscr{Y}=\chi+\vartheta\phi(\Upsilon),\quad \vartheta<1. \end{array}$$(32)

Therefore, ϒ can be also expressed as a function of χ as follows:

Y=χ+∑k=1∞ϑkk!dk−1dχk−1[ϕ(χ)]k.$$\begin{array}{} \displaystyle \mathscr{Y}=\chi+\sum_{k=1}^{\infty}\dfrac{\vartheta^k}{k!}\dfrac{d^{k-1}}{d \chi^{k-1}}[\phi(\chi)]^k. \end{array}$$(33)

3.1.1

Location of L₁

The solution of L₁ lies beyond the smaller mass as in Fig. (2). Therefore, r₁ − r₂ = 1, and hence, r₂ = μ − x − 1, r₁ = μ − x, which implies that

Fig. 2

∂r1∂x=∂r2∂x=−1,$$\begin{array}{} \displaystyle \dfrac{\partial r_1}{\partial x}=\dfrac{\partial r_2}{\partial x}=-1, \end{array}$$

provided that we take r₂ = r, r₁ = 1 + r, and substitute them into Eq. (29). After some simple calculations the value of ν will be written as

ν=a4r4+a5r5+a6r6+a7r7+O[r8],$$\begin{array}{} \displaystyle \nu=a_4r^4+a_5r^5+a_6r^6+a_7r^7+\mathscr{O}[r^8], \end{array}$$(34)

where the values a₄, a₅, a₆ and a₇ are given as in Appendix I.

Now, we use Lagrange’s inversion method to invert the above series and hence to express r as a function of ν which is written as

r=(1a4)14ν14−14(1a4)32a5ν+132(1a4)114(7a52−8a4a6)ν34+−a53+2a4a5a6−a42a74a44ν+O[ν54].$$\begin{array}{} \displaystyle r=\Bigg(\dfrac{1}{a_4}\Bigg)^{\frac{1}{4}} \nu^{\frac{1}{4}}-\dfrac{1}{4}\Bigg(\dfrac{1}{a_4}\Bigg)^{\frac{3}{2}}a_5 \sqrt{\nu}+\dfrac{1}{32}\Bigg(\dfrac{1}{a_4}\Bigg)^{\frac{11}{4}}(7a_5^2-8a_4a_6)\nu^{\frac{3}{4}}+\dfrac{-a_5^3+2a_4a_5a_6-a_4^2a_7}{4a_4^4}\nu+\mathscr{O}[\nu^{\frac{5}{4}}]. \end{array}$$(35)

Consequently, the location of L₁ will be given by

x1=−1+ν1+ν−(1a4)14ν14+14(1a4)32a5ν−132(1a4)114(7a52−8a4a6)ν34−−a53+2a4a5a6−a42a74a44ν.$$\begin{array}{} \displaystyle x_1=-1+\dfrac{\nu}{1+\nu}-\Bigg(\dfrac{1}{a_4}\Bigg)^{\frac{1}{4}}\nu^{\frac{1}{4}}+\dfrac{1}{4}\Bigg(\dfrac{1}{a_4}\Bigg)^{\frac{3}{2}}a_5\sqrt{\nu}-\dfrac{1}{32}\Bigg(\dfrac{1}{a_4}\Bigg)^{\frac{11}{4}}(7a_5^2-8a_4a_6)\nu^{\frac{3}{4}}-\dfrac{-a_5^3+2a_4a_5a_6-a_4^2a_7}{4a_4^4}\nu. \end{array}$$(36)

3.1.2

Location of L₂

The solution of L₂ lies between the two finite bodies as in Fig. (2), where r₁ + r₂ = 1. Hence, r₂ = x − μ + 1, r₁ = μ − x. This implies that

∂r1∂x=∂r2∂x=–1.$$\begin{array}{} \displaystyle \frac{{\partial {r_1}}}{{\partial x}} = \frac{{\partial {r_2}}}{{\partial x}} = -1. \end{array}$$

Thus, if we take r₂ = ρ, r₁ = 1 − ρ, and substitute them into Eq. (29), then the value of ν will be given by

[ν=b4ρ4+b5ρ5+b6ρ6+b7ρ7+O[ρ8],$$\begin{array}{} \displaystyle [\nu=b_4\rho^4+b_5\rho^5+b_6\rho^6+b_7\rho^7+\mathscr{O}[\rho^8], \end{array}$$

where the values b₄, b₅, b₆, and b₇ are governed as in Appendix I.

As mentioned previously, we use Lagrange’s inversion method in order to invert the above series and hence, express ρ as a function of ν, which can be written as

ρ=(1b4)14ν14−14(1b4)32b5ν+132(1b4)114(7b52−8b4b6)ν34+(−b53+2b4b5b6−b42b7)4b44ν+O[ν54].$$\begin{array}{} \displaystyle \rho=\left(\dfrac{1}{b_4}\right)^{\frac{1}{4}}\nu^{\frac{1}{4}}-\dfrac{1}{4}\left(\dfrac{1}{b_4}\right)^{\frac{3}{2}}b_5\sqrt{\nu}+\dfrac{1}{32}\left(\dfrac{1}{b_4}\right)^{\frac{11}{4}}(7b_5^2-8b_4b_6)\nu^{\frac{3}{4}}+\dfrac{(-b_5^3+2b_4b_5b_6-b_4^2b_7)}{4b_4^4}\nu+\mathscr{O}[\nu^{\frac{5}{4}}]. \end{array}$$(37)

Therefore, the location of L₂ could be written as

x2=−1+ν1+ν+(1b4)14ν14−14(1b4)32b5ν+132(1b4)114(7b52−8b4b6)ν34+(−b53+2b4b5b6−b42b7)4b44ν.$$\begin{array}{} \displaystyle x_2=-1+\dfrac{\nu}{1+\nu}+\left(\dfrac{1}{b_4}\right)^{\frac{1}{4}}\nu^{\frac{1}{4}}-\dfrac{1}{4}\left(\dfrac{1}{b_4}\right)^{\frac{3}{2}}b_5\sqrt{\nu}+\dfrac{1}{32}\left(\dfrac{1}{b_4}\right)^{\frac{11}{4}}(7b_5^2-8b_4b_6)\nu^{\frac{3}{4}}+\dfrac{(-b_5^3+2b_4b_5b_6-b_4^2b_7)}{4b_4^4}\nu. \end{array}$$(38)

3.1.3

Location of L₃

The solution of L₃ lies beyond the large mass as in Fig. (2), where r₂ − r₁ = 1. Hence, r₂ = 1 − μ + x, r₁ = x − μ. This implies that

∂r1∂x=∂r2∂x=1.$$\begin{array}{} \displaystyle \dfrac{\partial r_1}{\partial x}=\dfrac{\partial r_2}{\partial x}=1. \end{array}$$

Thus, if we take r₁ = 1 − +δ, r₂ = 2 + δ , and substitute them into Eq. (29), then v is as follows:

ν=c0+c1δ+c2δ2+O[δ3],$$\begin{array}{} \displaystyle \nu=c_0+c_1\delta+c_2\delta^2+\mathscr{O}[\delta^3], \end{array}$$(39)

where c₀, c₁ and c₂ are given as in Appendix I.

We also use Lagrange’s inversion method to invert the above series and to express δ as a function of v. Hence, the following expression yields:

δ=ν−c0c1−c2(ν−c0)2c13+O[(ν−c0)3].$$\begin{array}{} \displaystyle \delta=\dfrac{\nu-c_0}{c_1}-\dfrac{c_2(\nu-c_0)^2}{c_1^3}+\mathscr{O}[(\nu-c_0)^3]. \end{array}$$(40)

Thus, the location of L₃ will be controlled by

x3=1+ν1+ν+ν−c0c1−c2(ν−c0)2c13.$$\begin{array}{} \displaystyle x_3=1+\dfrac{\nu}{1+\nu}+\dfrac{\nu-c_0}{c_1}-\dfrac{c_2(\nu-c_0)^2}{c_1^3}. \end{array}$$(41)

Note that Eqs. (36), (38) and (41) represent semi-closed forms for the positions of collinear points. Furthermore, the small perturbation in the Coriolis force has no effect on these positions.

3.1.4

The triangular points

The solutions of Eqs. (24) and (25-26) when y ≠ 0 will give the triangular points L₄ and L₅. These solutions will lead to

f1(r1)=0=f2(r2).$$\begin{array}{} \displaystyle f_1(r_1)=0=f_2(r_2). \end{array}$$(42)

Therefore, we get

ψn2=(1r13+3(A1+A)2r15).$$\begin{array}{} \displaystyle \psi n^2=\left(\dfrac{1}{r_1^3}+\dfrac{3(A_1+A)}{2r_1^5}\right). \end{array}$$(43)

ψn2=(1r23+3(A2+A)2r25).$$\begin{array}{} \displaystyle \psi n^2=\left(\dfrac{1}{r_2^3}+\dfrac{3(A_2+A)}{2r_2^5}\right). \end{array}$$(44)

To study the effect of oblateness in the perturbed problem on the locations of the triangular points, we will assume that the three participating bodies are not oblate spheroids, namely, A₁ = A₂ = A = 0.

Hence Eqs. (43-44) will admit ri=1ψ3$\begin{array}{} \displaystyle r_i=\dfrac{1}{\sqrt[3]{\psi}} \end{array}$. Consequently if the three bodies are oblate, then the solution of Eqs. (43-44) will be controlled by

ri=1ψ3+πi,πi=1:i=1,2,$$\begin{array}{} \displaystyle r_i=\dfrac{1}{\sqrt[3]{\psi}}+\pi_i,\quad \pi_i=1: i=1,2, \end{array}$$(45)

where π_i is the additive displacement in r_i as a result of oblateness effects. Substituting Eq. (45) into (43-44) and restricting ourselves to the linear terms in π₁, π₂, A₁, A₂, and A, we obtain an appropriate approach to π_i as follows:

π1=−12(A1+A2)ψ−13+12(A1+A)ψ13.$$\begin{array}{} \displaystyle \pi_1=-\dfrac{1}{2}(A_1+A_2)\psi^{-\frac{1}{3}}+\dfrac{1}{2}(A_1+A)\psi^{\frac{1}{3}}. \end{array}$$(46)

π2=−12(A1+A2)ψ−13+12(A2+A)ψ13.$$\begin{array}{} \displaystyle \pi_2=-\dfrac{1}{2}(A_1+A_2)\psi^{-\frac{1}{3}}+\dfrac{1}{2}(A_2+A)\psi^{\frac{1}{3}}. \end{array}$$(47)

From Eqs. (16-17), the exact solutions of the triangular points could be written as

x=−12[1−2μ+r12−r22].$$\begin{array}{} \displaystyle x=-\dfrac{1}{2}\left[1-2\mu+r_1^2-r_2^2\right]. \end{array}$$(48)

y=±[(r12+r222)−(r12−r222)2−14]12.$$\begin{array}{} \displaystyle y=\pm\left[\left(\dfrac{r_1^2+r_2^2}{2}\right)-\left(\dfrac{r_1^2-r_2^2}{2}\right)^2-\dfrac{1}{4}\right]^{\frac{1}{2}}. \end{array}$$(49)

Substituting Eqs. (46-47) and (45) into both (48-49) and limiting ourselves to the previous assumptions, it holds that

x=μ−12−12(A1−A2).$$\begin{array}{} \displaystyle x=\mu-\dfrac{1}{2}-\dfrac{1}{2}(A_1-A_2). \end{array}$$(50)

y=±32[1−13(A1+A2−2A)−49(1+A1+A2+A)ε2].$$\begin{array}{} \displaystyle y=\pm \dfrac{\sqrt{3}}{2}\left[1-\dfrac{1}{3}(A_1+A_2-2A)-\dfrac{4}{9}(1+A_1+A_2+A)\varepsilon_2\right]. \end{array}$$(51)

Eqs. (50-51) state that the presence of parameters that represent the oblateness of the infinitesimal body, the small perturbations in Coriolis and centrifugal forces have no effect on the x−coordinate of the triangular points. While the y−coordinate of the triangular points may be affected by all the parameters of the perturbed forces except for the small perturbations in the Coriolis force.

The models which describe the restricted three-body problem are developed so far they provide a concise description regarding the third body motion. These models depend on the state variables that define the initial conditions and the variety of parameters that either affects the measurement process or dynamical motions, as well. There is a difficulty to solve these models directly for any choice of parameters from a given set of measurements, due to the great complexity underlying these models. This is the main reason for which our attention is focused on linearizing the dynamical system to obtain simplified expressions that could be handled more easily. To understand the possible motion of the third body in the vicinity of equilibrium points, the equations of motion should be linearized along the initial state vector. In other words, we will expand the right handside of motion equations around the equilibrium points. Consequently, the linear part in these equations is called the variational equations. We also assume that the variation vector is related to the initial state vector by the following expression:

r_=r_0+Δr_,$$\begin{array}{} \displaystyle \underline{r}=\underline{r}_0+\Delta \underline{r}, \end{array}$$(52)

where r ≡ (x, y),r₀ ≡ (x₀, y₀), and Δr ≡ (ξ, η), provided that ξ and η are small enough. Thus, substituting Eq. (52) into Eq. (20-21), and expanding them via a Taylor series, one obtain that

ξ¨−2nφη˙=Ωx0+11![ξ∂∂x+η∂∂y]Ωx0+12![ξ∂∂x+η∂∂y]2Ωx0+O[3].$$\begin{array}{} \displaystyle \ddot{\xi}-2n\varphi\dot{\eta}=\Omega_x^0+\dfrac{1}{1!}\left[\xi\dfrac{\partial}{\partial x}+\eta\dfrac{\partial}{\partial y}\right]\Omega_x^0+\dfrac{1}{2!}\left[\xi\dfrac{\partial}{\partial x}+\eta\dfrac{\partial}{\partial y}\right]^2\Omega_x^0+\mathscr{O}[3]. \end{array}$$(53)

η¨+2nφξ˙=Ωy0+11![ξ∂∂x+η∂∂y]Ωy0+12![ξ∂∂x+η∂∂y]2Ωy0+O[3],$$\begin{array}{} \displaystyle \ddot{\eta}+2n\varphi\dot{\xi}=\Omega_y^0+\dfrac{1}{1!}\left[\xi\dfrac{\partial}{\partial x}+\eta\dfrac{\partial}{\partial y}\right]\Omega_y^0+\dfrac{1}{2!}\left[\xi\dfrac{\partial}{\partial x}+\eta\dfrac{\partial}{\partial y}\right]^2\Omega_y^0+\mathscr{O}[3], \end{array}$$(54)

where 𝒪[3] stands from the third and higher order terms in ξ and η. When these terms and the second order terms are omitted, then Eqs. (53-54) give us the linear variational equations and ξ and η are called the variations. Indeed,

ξ¨−2nφη˙=ξΩxx0+ηΩxy0.$$\begin{array}{} \displaystyle \ddot{\xi}-2n\varphi\dot{\eta}=\xi\Omega_{xx}^0+\eta\Omega_{xy}^0. \end{array}$$(55)

η¨+2nφξ˙=ξΩxy0+ηΩyy0,$$\begin{array}{} \displaystyle \ddot{\eta}+2n\varphi\dot{\xi}=\xi\Omega_{xy}^0+\eta\Omega_{yy}^0, \end{array}$$(56)

where the subscripts x and y denote the second partial derivatives of Ω and superscript 0 means that these derivatives are calculated with respect to one of the five equilibrium points (x₀, y₀).

Eqs. (55-56) represent the linear differential system of second order in ℝ². Hence, we transform such a system in a first order system in ℝ⁴. To deal with, let us assume that ξ₁ = ξ, ξ₂ = η, ξ₃ = ξ₁, ξ₄ = ξ₂, and let also substitute them into Eqs. (55-56). Then the system can be written as follows:

ξ1.=ξ3.$$\begin{array}{} \displaystyle \dot{\xi_1}=\xi_3. \end{array}$$(57)

ξ2.=ξ4.$$\begin{array}{} \displaystyle \dot{\xi_2}=\xi_4. \end{array}$$(58)

ξ3.=Ωxx0ξ1+Ωxy0ξ2+2nφξ4.$$\begin{array}{} \displaystyle \dot{\xi_3}=\Omega_{xx}^0\xi_1+\Omega_{xy}^0\xi_2+2n\varphi\xi_4. \end{array}$$(59)

ξ4.=Ωxy0ξ1+Ωyy0ξ2−2nφξ3.$$\begin{array}{} \displaystyle \dot{\xi_4}=\Omega_{xy}^0\xi_1+\Omega_{yy}^0\xi_2-2n\varphi\xi_3. \end{array}$$(60)

The stability of libration points is a well-known problem explored alongn the classical literature. Therefore, in this section and also in forthcoming sections, we will follow the analysis based on Szebehely computations [29].

The periodic orbits have great significance generally in the celestial mechanics and specially in space dynamics. That significance makes necessity for studying the periodic orbits. There are some important reasons that motivated this study. First of all, note that the non-integrable dynamical systems prevent us to obtain a complete solution for any orbits unless they are asymptotic or periodic orbits. For any particular solution of the restricted three-body problem, one can always find out a periodic solution. We can also get periodic orbits via linearized solutions for the motion of infinitesimal body around the libration points. Periodic orbits play a crucial role in order to classify the different classes of orbits and also to reduce the dimensionality of the problem in the phase space. Additionally, they can be used as reference orbits.