Motion about equilibrium points in the Jupiter-Europa system with oblateness

The restricted three body problem (RTBP) describes the motion of a restricted mass (also known as infinitesimal mass), which is moving under the effect of gravitational field of two massive bodies (called as primaries), without disturbing the actual motion of these primaries about their common centre of mass [25]. The RTBP is a demonstrative model of many problems in celestial mechanics, and it is used as a base model to describe the dynamics of small bodies (asteroid, comet, satellite, spacecraft, etc.) in solar as well as in extra-solar planetary systems [3, 11].

In the study of RTBP, the perturbing forces such as radiation pressure, Poynting-Robertson drag, solar wind drag, oblateness of massive bodies, etc. have an effective influence on celestial bodies. A number of researchers have contributed to the study of motion of infinitesimal mass under the frame of RTBP, Chermnykh-like problem etc. through different approaches and techniques in the presence of different perturbations in the form of oblateness, radiation pressure force, disc's effect, drag forces etc. [4, 6, 12, 13, 24, 29]. AbdulRaheem and Singh [1] have studied the combined effect of Coriolis force, centrifugal force, radiation pressure force and oblateness of the primaries on the equilibrium points and found that only Coriolis force acts as a stabilising factor, whereas centrifugal force, radiation pressure force and oblateness act as destabilising factors that cause a contraction in the stability range of the triangular equilibrium points. A few authors [15, 23, 28] have discussed the effects of radiation pressure, triaxial rigid bodies, oblateness of massive bodies and presence of a disc during the analysis of equilibrium points and their stability. They have found a deviation in the positions of equilibrium points and stability region because of these perturbations.

In the solar system dynamics, the analysis of periodic orbits produces many important ideas to understand the general properties of different dynamical systems. Companys et al. [8] have determined the existence of quasi-periodic orbits near the equilibrium points of the Earth-Moon system. The quasi-periodic orbits play an essential role in the designing of new space missions due to which fuel consumption in the spacecraft reduces. Motion of asteroids in the surroundings of triangular equilibrium points L₄ and L₅ in the Sun-Jupiter system can be more understandable by studying the tadpole and horseshoe orbits about the L₄ and L₅ points, respectively [9, 16]. Rabe [22] and Taylor [26] have determined approximate initial conditions to compute a series of periodic orbits around the L₄ and L₅ points of the Sun-Jupiter system and they have found stable horseshoe-shaped periodic orbits.

In order to understand and to analyse the evolution of periodic orbits, Winter [27] studied the evolution of periodic and quasi-periodic orbits and their stability property in the phase space of the Earth-Moon system using a numerical technique known as Poincaré surfaces of section (PSS) and found that these are sensitive dependence on the initial conditions. Many authors [10, 14, 21] have analysed the nature and stability of the periodic and quasi-periodic orbits with the help of PSS under the frame of perturbed restricted three-body problem and have noticed the effects of perturbing parameters. Abouelmagd et al. [5] have investigated the analytical derivations to study periodic solutions for the two-body problem perturbed by the first zonal harmonic parameter by using different methods. Further, Abouelmagd and Ansari [3] have analysed the cases, which agree that Jacobian and energy conversation are constants, in the Sun-perturbed Earth-Moon system and have illustrated the equilibrium points and their stability and PSS. Again, Pathak et al. [17, 19] have studied interior higher order resonant orbits in the perturbed photo-gravitational restricted three-body problem and have found that there exist periodic orbits for seventh- and ninth-order resonance, which are passing around the Earth. They also analysed the effect of radiation pressure and oblateness on these orbits.

In this paper, we are interested to determine the tadpole and horseshoe orbits and to observe the dynamical behaviour of orbits using PSS technique under the frame of RTBP with oblate primaries. The paper is organised as follows: Section 2 presents the equations of motion of the infinitesimal mass in the presence of oblateness. Sections 3 and 4 are devoted to the determination of equilibrium points and their linear stability test under the effect of oblateness, respectively. Tadpole and horseshoe orbits in the neighbourhood of equilibrium points are discussed in Section 5. Computation of PSS and the analysis of the effect of perturbing parameters on the PSS are presented in Section 6. Finally, the paper is concluded in Section 7. Numerical and semi-analytic computations are performed by using MATHEMATICA software (Version 12.0). Accuracy and precision goals of the values are achieved by considering or taking the values up to seven decimal places.

Let m₁ and m₂ be the masses of bigger and smaller oblate primaries, respectively. Assume that the primaries rotate about their common centre of mass in the xy-plane due to mutual gravitational attraction and the infinitesimal mass moves under the gravitational field of the primaries in the orbital plane without influencing their motion as shown in Figure 1. Assume that the sum of masses of the primaries is the unit of mass and distance between the primaries is the unit of distance. The unit of time is such that G(m₁ + m₂) = 1. Since the primaries are considered as oblate spheroid therefore, the angular speed (mean motion) n is written [2, 28] as n=1+3(A21+A22)2 n = \sqrt {1 + {{3({A_{21}} + {A_{22}})} \over 2}} , where A_2i, i = 1, 2 represent the oblateness coefficients of the bigger and smaller primary [28]. The mass ratio is defined as μ=m2m1+m2 \mu = {{{m_2}} \over {{m_1} + {m_2}}} so that 0<μ≤12 0 < \mu \le {1 \over 2} . Let the coordinates of the bigger and smaller primary be (−μ, 0) and (1 − μ, 0), respectively and (x, y) be the position coordinate of the infinitesimal mass in the synodic frame Oxy. Therefore, the equations of motion of the infinitesimal mass P(x, y) in the xy-plane [18, 25] are written as (1) x¨−2ny˙=Ωx, \ddot x - 2n\dot y = {\Omega _x}, (2) y¨+2nx˙=Ωy, \ddot y + 2n\dot x = {\Omega _y}, where the effective potential Ω, experienced by the infinitesimal mass is given as (3) Ω=n22(x2+y2)+(1−μ)r1(1+A212 r12)+μr2(1+A222 r22) \Omega = {{{n^2}} \over 2}({x^2} + {y^2}) + {{(1 - \mu)} \over {{r_1}}}\left({1 + {{{A_{21}}} \over {2{\kern 1pt} r_1^2}}} \right) + {\mu \over {{r_2}}}\left({1 + {{{A_{22}}} \over {2{\kern 1pt} r_2^2}}} \right) and the respective distances r₁ and r₂ of the infinitesimal mass from the bigger and smaller primaries are defined as r1=(x+μ)2+y2 {r_1} = \sqrt {{{(x + \mu)}^2} + {y^2}} and r2=(x+μ−1)2+y2 {r_2} = \sqrt {{{(x + \mu - 1)}^2} + {y^2}} .

Fig. 1

Now, multiplying Eq. (1) by ẋ and Eq. (2) by ẏ, respectively and then adding them, we get x˙ x¨+y˙ y¨=Ωx x˙+Ωy y˙. \dot x{\kern 1pt} \ddot x + \dot y{\kern 1pt} \ddot y = {\Omega _x}{\kern 1pt} \dot x + {\Omega _y}{\kern 1pt} \dot y.

On integration, we obtain (4) x˙2+y˙22=∫∂ Ω (x, y)∂ t−C⇒x˙2+y˙2=2 Ω−C, {{{{\dot x}^2} + {{\dot y}^2}} \over 2} = \int {{\partial {\kern 1pt} \Omega {\kern 1pt} (x,{\kern 1pt} y)} \over {\partial {\kern 1pt} t}} - C \Rightarrow {\dot x^2} + {\dot y^2} = 2{\kern 1pt} \Omega - C, where, C is constant of integration. Eq. (4) represents a well-known integral called as Jacobi integral in which C is known as Jacobi constant.

In this section, we shall discuss the existence of those points at which the velocity of the infinitesimal mass become zero. These points can be obtained by solving Ω_x = 0, Ω_y = 0, simultaneously i.e. (5) Ωx=n2 x−(1−μ)(x+μ)r13(1+3 A212 r12)−μ(x+μ−1)r23(1+3 A222 r22)=0, {\Omega _x} = {n^2}{\kern 1pt} x - {{(1 - \mu)(x + \mu)} \over {r_1^3}}\left({1 + {{3{\kern 1pt} {A_{21}}} \over {2{\kern 1pt} r_1^2}}} \right) - {{\mu (x + \mu - 1)} \over {r_2^3}}\left({1 + {{3{\kern 1pt} {A_{22}}} \over {2{\kern 1pt} r_2^2}}} \right) = 0, (6) Ωy=n2 y−(1−μ) yr13(1+3 A212 r12)−μ yr22(1+3 A222 r22)=0, {\Omega _y} = {n^2}{\kern 1pt} y - {{(1 - \mu){\kern 1pt} y} \over {r_1^3}}\left({1 + {{3{\kern 1pt} {A_{21}}} \over {2{\kern 1pt} r_1^2}}} \right) - {{\mu {\kern 1pt} y} \over {r_2^2}}\left({1 + {{3{\kern 1pt} {A_{22}}} \over {2{\kern 1pt} r_2^2}}} \right) = 0, for space coordinates x and y.

Let the co-ordinate of the equilibrium point be denoted by (x₀,y₀). Suppose α and β denote small displacements of the infinitesimal mass from the equilibrium point (x₀,y₀), then position of infinitesimal mass at any future time is given as x = x₀ + α, y = y₀ + β. We assume that the displacements α and β are sufficiently small so that the Taylor's expansions of Ω_x and Ω_y about the equilibrium point (x₀, y₀) up to linear order terms in α, β, can be written as Ωx=α Ωxx0+β Ωxy0,Ωy=α Ωyx0+β Ωyy0. \matrix{{{\Omega _x} = \alpha {\kern 1pt} \Omega _{xx}^0 + \beta {\kern 1pt} \Omega _{xy}^0,} \hfill \cr {{\Omega _y} = \alpha {\kern 1pt} \Omega _{yx}^0 + \beta {\kern 1pt} \Omega _{yy}^0.} \hfill \cr}

Therefore, the equations of motion (1) and (2) in the neighbourhood of (x₀,y₀) can be written by using Taylor series expansion in the linear order as (15) α¨−2nβ˙=α Ωxx0+β Ωxy0, \ddot \alpha - 2n\dot \beta = \alpha {\kern 1pt} \Omega _{xx}^0 + \beta {\kern 1pt} \Omega _{xy}^0, (16) β¨+2nα˙=α Ωyx0+β Ωyy0, \ddot \beta + 2n\dot \alpha = \alpha {\kern 1pt} \Omega _{yx}^0 + \beta {\kern 1pt} \Omega _{yy}^0, where Ωxx0 \Omega _{xx}^0 , Ωxy0 \Omega _{xy}^0 and Ωyy0 \Omega _{yy}^0 are second order partial derivatives at equilibrium point (x₀, y₀) and the expressions of Ω_xx, Ω_xy, Ω_yy are given by Ωxx=n2+(1−μ)r13[6 (x+μ)2−3 A212 r12+15 A21 (x+μ)22 r14−1]+μr23[6 (x+μ+1)2−3 A222 r22+15 A22(x+μ+1)22 r24−1],Ωyy=n2+(1−μ)r13[6 y2−3 A212 r12+15 A21 y22 r14−1]+μr23[6 y2−3 A222 r22+15 A22 y22 r24−1],Ωxy=(1−μ)r15[3 y(x+μ)+15 A21 (x+μ) y2 r12]+μr25[15 A22 (x+μ−1) y2 r22+3 (x+μ−1) y]. \matrix{{{\Omega _{xx}} = {n^2} + {{(1 - \mu)} \over {r_1^3}}\left[ {{{6{\kern 1pt} {{(x + \mu)}^2} - 3{\kern 1pt} {A_{21}}} \over {2{\kern 1pt} r_1^2}} + {{15{\kern 1pt} {A_{21}}{\kern 1pt} {{(x + \mu)}^2}} \over {2{\kern 1pt} r_1^4}} - 1} \right]} \cr {+ {\mu \over {r_2^3}}\left[ {{{6{\kern 1pt} {{(x + \mu + 1)}^2} - 3{\kern 1pt} {A_{22}}} \over {2{\kern 1pt} r_2^2}} + {{15{\kern 1pt} {A_{22}}{{(x + \mu + 1)}^2}} \over {2{\kern 1pt} r_2^4}} - 1} \right],} \cr {{\Omega _{yy}} = {n^2} + {{(1 - \mu)} \over {r_1^3}}\left[ {{{6{\kern 1pt} {y^2} - 3{\kern 1pt} {A_{21}}} \over {2{\kern 1pt} r_1^2}} + {{15{\kern 1pt} {A_{21}}{\kern 1pt} {y^2}} \over {2{\kern 1pt} r_1^4}} - 1} \right] + {\mu \over {r_2^3}}\left[ {{{6{\kern 1pt} {y^2} - 3{\kern 1pt} {A_{22}}} \over {2{\kern 1pt} r_2^2}} + {{15{\kern 1pt} {A_{22}}{\kern 1pt} {y^2}} \over {2{\kern 1pt} r_2^4}} - 1} \right],} \cr {{\Omega _{xy}} = {{(1 - \mu)} \over {r_1^5}}\left[ {3{\kern 1pt} y(x + \mu) + {{15{\kern 1pt} {A_{21}}{\kern 1pt} (x + \mu){\kern 1pt} y} \over {2{\kern 1pt} r_1^2}}} \right] + {\mu \over {r_2^5}}\left[ {{{15{\kern 1pt} {A_{22}}{\kern 1pt} (x + \mu - 1){\kern 1pt} y} \over {2{\kern 1pt} r_2^2}} + 3{\kern 1pt} (x + \mu - 1){\kern 1pt} y} \right].} \cr}

Now, we investigate the motion of infinitesimal body about the equilibrium point (x₀, y₀) in xy-plane, for this we suppose α = Ae^{λ t}, β = Be^{λ t} where A, B, λ are parameters to be determined. Thus Eqs (15) and (16) reduces to A (λ2−Ωxx0)+B (−2nλ−Ωxy0)=0A (2nλ−Ωyx0)+B (λ2−Ωyy0)=0. \matrix{{A{\kern 1pt} ({\lambda ^2} - \Omega _{xx}^0) + B{\kern 1pt} (- 2n\lambda - \Omega _{xy}^0) = 0} \hfill \cr {A{\kern 1pt} (2n\lambda - \Omega _{yx}^0) + B{\kern 1pt} ({\lambda ^2} - \Omega _{yy}^0) = 0.} \hfill \cr}

The above equations have a non-trivial solution for A and B if |λ2−Ωxx0−2 n λ−Ωxy02 n λ−Ωyx0λ2−Ωyy0|=0, \left| {\matrix{{{\lambda ^2} - \Omega _{xx}^0} & {- 2{\kern 1pt} n{\kern 1pt} \lambda - \Omega _{xy}^0} \cr {2{\kern 1pt} n{\kern 1pt} \lambda - \Omega _{yx}^0} & {{\lambda ^2} - \Omega _{yy}^0} \cr}} \right| = 0, which gives a bi-quadratic equation in λ known as characteristic equation (17) λ4+P1 λ2+Q1=0, {\lambda ^4} + {P_1}{\kern 1pt} {\lambda ^2} + {Q_1} = 0, where P1=4 n2−Ωxx0−Ωyy0 {P_1} = 4{\kern 1pt} {n^2} - \Omega _{xx}^0 - \Omega _{yy}^0 and Q1=Ωxx0Ωyy0−(Ωxy0)2 {Q_1} = \Omega _{xx}^0\Omega _{yy}^0 - {(\Omega _{xy}^0)^2} . The characteristic roots of characteristic Eq (17) are given as (18) λ=±−P1±P12−4Q12, \lambda = \pm \sqrt {{{- {P_1} \pm \sqrt {P_1^2 - 4{Q_1}}} \over 2}}, which are function of μ; A₂₁ and A₂₂. For the linear stability of equilibrium points, the nature of roots play a key role, which depends on the sign of discriminant P12−4 Q1 P_1^2 - 4{\kern 1pt} {Q_1} of Eq. (17). Hence, in general, the solution of Eqs (15) and (16) are expressed as (19) α(t)=A1eλ1 t+A2eλ2 t+A3eλ3 t+A4eλ4 t, \alpha (t) = {A_1}{e^{{\lambda _1}{\kern 1pt} t}} + {A_2}{e^{{\lambda _2}{\kern 1pt} t}} + {A_3}{e^{{\lambda _3}{\kern 1pt} t}} + {A_4}{e^{{\lambda _4}{\kern 1pt} t}}, (20) β(t)=B1eλ1 t+B2eλ2 t+B3eλ3 t+B4eλ4 t, \beta (t) = {B_1}{e^{{\lambda _1}{\kern 1pt} t}} + {B_2}{e^{{\lambda _2}{\kern 1pt} t}} + {B_3}{e^{{\lambda _3}{\kern 1pt} t}} + {B_4}{e^{{\lambda _4}{\kern 1pt} t}}, where coefficients B_j are related to A_j by the relation as Bj=(λj2−Ωxx02 n λj)Aj ; j=1,2,3,4. {B_j} = \left({{{\lambda _j^2 - \Omega _{xx}^0} \over {2{\kern 1pt} n{\kern 1pt} {\lambda _j}}}} \right){A_j}{\kern 1pt} ;\quad j = 1,2,3,4.

4.2

Linear Stability of Triangular Equilibrium Points

In case of triangular equilibrium points, y ≠ 0. Therefore, the critical value of mass ratio can be obtained by solving the discriminant D=P12−4 Q1 D = P_1^2 - 4{\kern 1pt} {Q_1} by equating it to zero, i.e. (23) [2n2−(1−μ)(2r13+6 A21r15)+μ(2r23+6 A22r25)]2−[1+2 n2+(1−μ)(2r13+6 A21r15)+μ(2r23+6 A22r25)]=0. \matrix{{{{\left[ {2{n^2} - (1 - \mu)\left({{2 \over {r_1^3}} + {{6{\kern 1pt} {A_{21}}} \over {r_1^5}}} \right) + \mu \left({{2 \over {r_2^3}} + {{6{\kern 1pt} {A_{22}}} \over {r_2^5}}} \right)} \right]}^2} -} \hfill \cr {\left[ {1 + 2{\kern 1pt} {n^2} + (1 - \mu)\left({{2 \over {r_1^3}} + {{6{\kern 1pt} {A_{21}}} \over {r_1^5}}) + \mu ({2 \over {r_2^3}} + {{6{\kern 1pt} {A_{22}}} \over {r_2^5}}} \right)} \right] = 0.} \hfill \cr}

In absence of the oblateness, i.e. when A₂₁ = 0, A₂₂ = 0, the Eq. (23) gives the critical mass ratio μ_c = 0.038529... [25], whereas in presence of oblateness (in particular at A₂₁ = 0.000285, A₂₂ = 0.007189), the critical value of mass ratio becomes μ_c = 0.024758. Thus, there arise three possible case as follows:

Case 1. When μc<μ<12 {\mu _c} < \mu < {1 \over 2} ; i.e. P² − 4Q < 0

Let D=+P2−4Q D = + \sqrt {{P^2} - 4Q} , then Eq. (18) infer that four roots of the form a + ib. Due to positive real part, solutions (19) and (20) will contain exponential term and hence, it is clear that these type of roots make the triangular points instable.

Case 2. When μ = μ_c; i.e. P² − 4Q = 0

This case implies P12−4 Q1=0 P_1^2 - 4{\kern 1pt} {Q_1} = 0 , thus Eq. (18) yields λ1,2=−i P12, λ3,4=i P12 {\lambda _{1,2}} = - i{\kern 1pt} \sqrt {{{{P_1}} \over 2}},\quad {\lambda _{3,4}} = i{\kern 1pt} \sqrt {{{{P_1}} \over 2}} which show that the roots are multiples of pure imaginary and hence, there appears secular terms in the solution. Therefore, the triangular points are unstable in this case too.

Case 3. When 0 < μ < μ_c; i.e. P² − 4Q > 0

In this case, we found that the roots are of the form λ1,2=±i ω1, λ3,4=±i ω2, {\lambda _{1,2}} = \pm i{\kern 1pt} {\omega _1},\quad {\lambda _{3,4}} = \pm i{\kern 1pt} {\omega _2}, where ω1,2=−P1±P12−4 Q2, {\omega _{1,2}} = \sqrt {{{- {P_1} \pm \sqrt {P_1^2 - 4{\kern 1pt} Q}} \over 2}}, are real numbers. Thus, the resulting motion will be a combination of two periodic motions. Motion corresponds to the period T=2 π|λ1,2|=2 πω1 T = {{2{\kern 1pt} \pi} \over {|{\lambda _{1,2}}|}} = {{2{\kern 1pt} \pi} \over {{\omega _1}}} is known as short periodic motion and that corresponds to period T=2 π|λ3,4|=2 πω2 T = {{2{\kern 1pt} \pi} \over {|{\lambda _{3,4}}|}} = {{2{\kern 1pt} \pi} \over {{\omega _2}}} is called long periodic motion.

Further, we found that all characteristic roots at A₂₁ = 0.000285, A₂₂ = 0.007198 are purely imaginary for mass ratio 0 < μ < 0.0247585 (Figure 2) which indicates that L_4,5 are stable there. Also, we have observed a variation in the nature of characteristic roots with respect to oblateness A₂₁, A₂₂ as given in Figure 3 and noticed that the oblateness of Europa affects more than that of oblateness of Jupiter on the stability region.

Fig. 2

Fig. 3

Suppose, the coefficients A_j and B_j in solutions (19) and (20) are of the form A_j = a_j + ib_j and B_j = c_j + id_j where j = 1,2,3,4. For the motion of infinitesimal body, the displacement components α, β and the velocity components α˙ \dot \alpha , β˙ \dot \beta must be real. Since, the coefficients of exponential terms in solutions (19) and (20) exist in complex conjugate pairs [16] as follows Re(A1,A2)=a1, Re(A3,A4)=a2,Im(A1,−A2)=b1, Im(A3,−A4)=b2,Re(B1,B2)=c1, Re(B3,B4)=c2,Im(B1,−B2)=d1, Im(A3,−A4)=d2. \matrix{{Re({A_1},{A_2}) = {a_1},\quad Re({A_3},{A_4}) = {a_2},} \cr {Im({A_1}, - {A_2}) = {b_1},\quad Im({A_3}, - {A_4}) = {b_2},} \cr {Re({B_1},{B_2}) = {c_1},\quad Re({B_3},{B_4}) = {c_2},} \cr {Im({B_1}, - {B_2}) = {d_1},\quad Im({A_3}, - {A_4}) = {d_2}.} \cr}

Therefore, after simplification solutions (19) and (20) reduces to (24) α(t)=2 a1 cosω1 t+2 a2 cosω2 t−2 b1 sinω1 t−2 b2 sinω2 t, \alpha (t) = 2{\kern 1pt} {a_1}{\kern 1pt} \cos {\omega _1}{\kern 1pt} t + 2{\kern 1pt} {a_2}{\kern 1pt} \cos {\omega _2}{\kern 1pt} t - 2{\kern 1pt} {b_1}{\kern 1pt} \sin {\omega _1}{\kern 1pt} t - 2{\kern 1pt} {b_2}{\kern 1pt} \sin {\omega _2}{\kern 1pt} t, (25) β(t)=2 c1 cosω1 t+2 c2 cosω2 t−2 d1 sinω1 t−2 d2 sinω2 t. \beta (t) = 2{\kern 1pt} {c_1}{\kern 1pt} \cos {\omega _1}{\kern 1pt} t + 2{\kern 1pt} {c_2}{\kern 1pt} \cos {\omega _2}{\kern 1pt} t - 2{\kern 1pt} {d_1}{\kern 1pt} \sin {\omega _1}{\kern 1pt} t - 2{\kern 1pt} {d_2}{\kern 1pt} \sin {\omega _2}{\kern 1pt} t.

In particular at L₄ (0.4965504, 0.8638773), μ = 0.01, A₂₁ = 0.000285 and A₂₂ = 0.007198, the perturbed solution α and β with initial conditions α(0) = 10⁻⁴, β (0) = 10⁻⁴ and α˙(0)=0 \dot \alpha (0) = 0 , β˙(0)=0 \dot \beta (0) = 0 , are expressed as (26) α(t)≈2×10−4 (1.91 cosω1 t−1.41 cosω2 t−35.96 sinω1 t−19.45 sinω2 t), \alpha (t) \approx 2 \times {10^{- 4}}{\kern 1pt} \left({1.91{\kern 1pt} \cos {\omega _1}{\kern 1pt} t - 1.41{\kern 1pt} \cos {\omega _2}{\kern 1pt} t - 35.96{\kern 1pt} \sin {\omega _1}{\kern 1pt} t - 19.45{\kern 1pt} \sin {\omega _2}{\kern 1pt} t} \right), (27) β(t)≈2×10−4 (−11.89 cosω1 t−11.24 cosω2 t−18.84 sinω1 t−8.78 sinω2 t), \beta (t) \approx 2 \times {10^{- 4}}{\kern 1pt} \left({- 11.89{\kern 1pt} \cos {\omega _1}{\kern 1pt} t - 11.24{\kern 1pt} \cos {\omega _2}{\kern 1pt} t - 18.84{\kern 1pt} \sin {\omega _1}{\kern 1pt} t - 8.78{\kern 1pt} \sin {\omega _2}{\kern 1pt} t} \right), where, ω_1,2 ≈ ±0.355 and ω_3,4 ≈ ±0.948. The bounded and periodic nature of solutions (26) and (27) infer that motion in the neighbourhood of L_4,5 are stable. The semi analytical solution α(t) and β (t) graphically can be seen in Figure 4.

Fig. 4

In the RTBP, there exist different types of motion of an infinitesimal mass in the vicinity of the L_4,5, valid only for small magnitude of displacement from the equilibrium point. Rabe [22] has computed the approximate initial conditions to analyse the periodic orbits close to L_4,5. He has found two set of initial conditions for horseshoe orbits and for long tadpole-shaped orbits around the L_4,5. Also, Taylor [26] briefly has described the smooth horseshoe orbits around L_4,5 for different initial conditions in the Sun-Jupiter system without perturbation. We are interested to compute the tadpole and horseshoe periodic orbits in the Jupiter-Europa system to analyse the effect of oblateness on these orbits. To determine the tadpole and horseshoe orbits in Jupiter-Europa system, we have followed the method described in [16, 26].

We assume that the infinitesimal mass P(x, y) lie on the straight line joining the Jupiter (bigger primary) and L₄ or L₅, which is the closest approach for P towards L₄ or L₅. A good approximation for periodic orbit about L_4,5 is that the magnitude of velocity V of the infinitesimal mass changes slowly in such a way that ∂V∂t {{\partial V} \over {\partial t}} and ∂2V∂2t {{{\partial ^2}V} \over {{\partial ^2}t}} approaches to zero. Let the distance from the infinitesimal body P to the point L₄ be k at any instant, then the coordinates of infinitesimal mass in the rotating frame are x=(1+k)2−μ x = {{(1 + k)} \over 2} - \mu and y=3(1+k)2 y = {{\sqrt 3 (1 + k)} \over 2} . For different value of k, we generate the set of initial conditions. Thus, integrating the equations of motion (1) and (2) in rotating frame with these initial conditions under the presence of oblateness of primaries in Jupiter-Europa system, we obtain tadpole and horseshoe orbits. Numerical results for different values of k are shown in Figures 5 and 6, in which we have acquired tadpole and horseshoe orbits around the L₄ and L₅ in the Jupiter-Europa system in presence and absence of oblateness of Jupiter and Europa. For the initial conditions corresponding to the value of k = 0.0001, we have found a narrow tadpole orbits around L₄ (also, similar tadpole orbits about L₅ are obtained which are not presented here) at μ = 0.0000251 and A₂₁ = A₂₂ = 0. We have analysed the effect of oblateness on tadpole and horseshoe orbits as plotted in Figures 5(b,c,d) and 6(b,c,d). We have found that tadpole orbits gets elongated gradually due to oblateness of Jupiter (A₂₁ = 0.0285), but a slight increment in the oblateness of Europa (A₂₂ = 0.0007189) results in more elongation of the tadpole orbits (Figures 5b, 5c, 5d) with small deviation in amplitude of loops. On the other hand, in Figure 6(a–d), the horseshoe orbits about L₄ and L₅ at k = 0.005 are displayed. From Figure 6(a–d), it is noticed that in the presence of oblateness, the number of loops in the horseshoe orbits get varied. Thus, due to oblateness of the primaries, the shape and size of the tadpole and horseshoe orbits in the Jupiter-Europa system deviated from the classical one [16, 26].

Fig. 5

Fig. 6

The position coordinate (x, y) and corresponding velocity (ẋ, ẏ) at any instant corresponds to a point in a four-dimensional phase space, which can be reduced to three-dimensional phase space due to the existence of Jacobi constant C Eq. (4). In order to compute surfaces of section, we fix the plane y = 0 in the three dimensional space to get the values of x and ẋ which can be plotted with respect to time whenever the trajectory crosses the plane y = 0. This technique is called PSS or Poincaré map [20]. The information, whether the trajectory is regular or chaotic in nature, can be obtained with the help of PSS. Also, we can identify the periodic and quasi-periodic orbits in the phase space by using PSS technique. In the surface of section, the formation of chain of islands corresponds to a single quasi-periodic trajectory. The evolution of dusty regions between the curves and the chain of islands represents chaotic zones.

We have computed the PSS at different values of mass ratio μ, Jacobi constant C_j and oblateness coefficients A₂₁, A₂₂ by the method described in [7,16,27]. The starting value of x are taken in such a way that ẏ > 0 for each value of C_j with y = ẋ = 0. For the range of Jacobi constant C_j between 2.55 and 3.30, we have integrated the equations of motion with starting conditions where the distance between the two consecutive starting conditions is chosen as 0.01 at each C_j. Different surfaces of section are generated as shown in Figures 7–10. From Figure 7, we can see that the chaotic zones increases with the decrease in the mass ratio μ. Also, we have plotted the surface of section for Jupiter-Europa system for μ = 0.0000251 (Figure 7). For 2.55 < C_j < 3.04, it is observed that chain of islands corresponds to a quasi-periodic trajectory increases gradually and formation of the chaotic zones starts (Figure 8). Further, increasing the value of C_j, there remains chain of islands corresponding to quasi-periodic orbits, but reduces the chaotic regions. As, the oblateness of primaries A₂₁ and A₂₂ varies, the structural region and islands changes in character and one particular structure may become more dominant as shown in Figures 9 and 10, but reduction in the chaotic zones is more due to oblateness of smaller primary A₂₂ than that of bigger primary A₂₁. In presence of oblateness, it is found that the curves or chain of islands representing quasi-periodic orbits predominates the chaotic regions.

Fig. 7

Fig. 8

Fig. 9

Fig. 10

We have considered the Jupiter-Europa system with oblateness of Jupiter and Europa and have found the equations of motion of infinitesimal mass. We have determined the equilibrium points and examined its linear stability under the influence of oblateness of both the primaries. It is found that due to oblateness of Jupiter and Europa, the stability range of the triangular equilibrium point is deviated from 0 < μ_c < 0.038529... [25] to 0 < μ_c < 0.024758... however, collinear equilibrium points are unstable for 0 < μ ≤ 1/2 as in classical case [25]. Moreover, positions of equilibrium points L_1,2,3,4,5 have shown deviations from their respective positions in classical cases (Table 1). Again, we have reckoned tadpole and horseshoe orbits in the Jupiter-Europa system under the influence of oblateness for different values of k, which are displayed in Figures 5–6. For the initial conditions corresponding to the value of k = 0.0001, we have found narrow tadpole orbits around L₄ at = 0.0000251 and A₂₁ = A₂₂ = 0. The effect of oblateness on tadpole and horseshoe orbits are analyzed (Figures 5b, 5c, 5d and 6b, 6c, 6d) and it is observed that tadpole orbits gets elongate gradually due to oblateness of Jupiter (A₂₁ = 0.0285), but a slight increment in the oblateness of Europa (A₂₂ = 0.0007189) results in more elongation of the tadpole orbits (Figures 5b, 5c, 5d) with small deviation in amplitude of loops. On the other hand, from the horseshoe orbits about L₄ and L₅ at k = 0.005 (Figures 6), it is noticed that in the presence of oblateness, a small difference in the loops occurs (Figures 6b, 6c, 6d). Because of the oblateness of the primaries, the results on tadpole and horseshoe orbits of the model in hand, shows a deviation from the results of classical case [16, 26]. Further, we have described the nature of orbits in presence of oblateness by using Poincaré surfaces of sections technique and computed the surfaces of section at different values of, C_j, A₂₁ and A₂₂ (Figures 7–10). It is noticed that the oblateness of the primaries minimize the chaotic regions and induce the formation of a chain of islands that corresponds to quasi-periodic orbits. Also, we have found that the oblateness parameters A₂₁ and A₂₂ affect the structure and nature of orbits. The PSS in our model can be degraded to that of PSS without the oblateness of primaries [7]. Finally, it is concluded that shape and size of the primaries play a significant role in the study of motion of the infinitesimal mass. These results may help to the study of more generalized problem of three bodies in space.