Oblateness Effect of Saturn on Halo Orbits of L1 and L2 in Saturn-Satellites Restricted Three-Body Problem

The Circular Restricted Three-Body Problem (CRTBP) with more massive primary as an oblate spheroid with its equatorial plane coincident with the plane of motion of the primaries is considered to generate the halo orbits around L1 and L2 for the seven satellites (Mimas, Enceladus, Tethys, Dione, Rhea, Titan and Iapetus) of Saturn in the frame work of CRTBP. It is found that the oblateness effect of Saturn on the halo orbits of the satellites closer to Saturn has significant effect compared to the satellites away from it. The halo orbits L1 and L2 are found to move towards Saturn with oblateness.


Introduction
Three-Body Problem formulated by Newton provided route to the analysis of closed form analytical solution. This solution remains elusive even today, as one has never been found for the three-body problem. Euler developed the restricted problem using a rotating frame in the 1770s and located collinear points. Along with Euler, Lagrange considered this form of the three-body problem and calculated the locations of equilateral points, often known as libration or Lagrange points. Jacobi studied the circularrestricted problem (CRTBP) and found that an integral of motion exists. Plummer [1], using an approximate, second-order analytical solution to the differential equations in the circular restricted three-body problem, produced a family of two-dimensional periodic orbits near the collinear libration points. Farquhar  analytical investigation into a class of periodic three-dimensional trajectories around the collinear points known as halo orbits. These trajectories are associated with the collinear points, and are a special case of the more general libration point orbits frequently designated as Lissajous orbits. Kamel and Farquhar [2] developed analytical approximations for quasi-periodic solutions associated with L2, the Earth-Moon libration point on the far side of the Moon. Richardson and Cary [3] derived a third-order approximation for motion near the interior Sun-Earth liberation point in the restricted problem.
Mission design utilizing the halo orbits become more challenging and several works had been established since then in this interesting area. Some of the important contributions are by Farquhar et al. [4], Richardson [5], Huber et al. [6], Gomez et al. [7] [8], Rausch [9], Nakamiya et al. [10], Koon et al. [11] and Nath & Ramanan [12]. Considering the motion in the vicinity of the collinear points in CRTBP, Breakwell and Brown [13] numerically extended the work of Farquhar and Kamel [2] to produce a family of periodic halo orbits. The discovery of a set of stable orbits in Breakwell and Brown's halo family motivated a search for stable orbits in the families associated with all the three collinear points by Howell [14] [15] and Howell et al. [16].
The subject of periodic solutions of the CRTBP has received enormous attention in the past few decades. Since the late twentieth century until today, enormous amount of research has enriched the study of CRTBP, but the influence of the various perturbing forces has not been studied in many of such interesting problems. The classical model does not account for some of the perturbing forces such as oblateness, solar radiation pressure, Poynting-Robertson drag effects and variation in the mass of the primaries.
Some of significant works in RTBP with oblateness effects are done by Sharma and Subba Rao [17] [18], Subba Rao and Sharma [19]- [21] and Sharma [22]- [25] by considering the more massive primary as an oblate spheroid with its equatorial plane co-incident with the plane of motion of the primaries. Danby [26], Papadakis [27], Kalantonis et al. [28], Raheem et al. [29], Stuchi et al. [30] discussed the restricted three-body problem with one or two bodies as oblate spheroids. The perturbing force due to the oblateness of Saturn is comparable with the perturbing force due to the gravitational attraction of the Sun in the Saturn-Satellites systems. The inclusion of oblateness effect has shown significant improvement in the theories of motion of certain satellites in the solar system [31] (Oberti and Vienne).
In the present study, we consider the restricted three-body problem by considering the more massive primary as an oblate spheroid with its equatorial plane coincident with the plane of motion (Sharma and Subba Rao [17]). We utilize Newton's method of differential correction (Mireles [32]; Eapen and Sharma [33]; Nishanth and Sharma [34]) to compute the halo orbits numerically about the Lagrangian points L1 and L2 in seven of the Saturn-Satellites systems.

Equations of Motion
The equations of motion for the restricted three-body problem are considered with 1 r , 2 r and r as the positions of three bodies with masses 1 m (more massive oblate sphe-roid-Saturn), 2 m (Satellite of Saturn) and m (spacecraft).
The origin of the co-ordinate system is the barycentre of the two primaries with the more massive primary lying to the left of the origin and the smaller primary to the right as shown in Figure 1. For scaling purpose, the distance between the two primaries is taken as unity, the sum of masses of the primaries ( ) 1 2 m m + is assumed unity, and the gravitational constant is assumed to be unity. Then the mean motion of the primaries is AE and AP are dimensional equatorial and polar radii of the more massive primary and R is the distance between the primaries.
The non-dimensional mass ratio ( ) µ is defined as the ratio of the mass of the smaller primary to the sum of masses of the primaries i.e. , respectively (Szebehely [35]; Sharma [17]; Bhatnagar and Chawla [36]).
The three-dimensional equations of motion are: where U is the pseudo potential of the system and is given as In the above expression 1 r and 2 r are the position vectors from the more massive and smaller primaries to the particle, respectively.

Liberation Points and Halo Orbits
From the equations of motion (1)-(3), it is apparent that an equilibrium solution exists relative to the rotating frame when the partial derivative of the pseudo potential function are all zero, i.e. 0 U ∆ = or , , 0 x y z =    . These points correspond to the positions in the rotating frame at which the gravitational force and the centrifugal force associated with the rotation of the synodic reference frame cancel, with the result that a particle positioned at one of these points appears stationary in the synodic frame. Also at the collinear points, Richardson's third-order approximation provides a deep qualitative insight. The approximate solution is sufficient for generating accurate motion near L1 and L2. Analytical approximation need to be combined with numerical techniques to generate a halo orbit accurate enough for mission design. In the present study to generate the halo orbits, we use analytical approximation as the first guess for the differential correction process, we have modified the third-order approximation of Thurman and Worfolk [37] by considering the more massive primary as an oblate spheroid with the help of Lindstedt-Poincaré method.

Analytical Approximation
For obtaining an analytical solution, following Tiwary and Kushvah [37], the origin is transferred to the Lagrangian points L 1 and L 2 and the transformation is given by The equations of motion can be written as The upper sign in the above equations depicts the Lagrangian point L1 and the lower sign corresponds to L2.
The usage of Legendre polynomials can result in some computational advantages, when non-linear terms are considered. The distance between these Lagrangian points and the smaller primary is considered to be the normalized unit as in Koon et al. [11] [38] and [39].
The non-linear terms are expanded by using the following formula as given by [11]: The above formula is used for expanding the non-linear terms in the equations of motion. The equations of motion after substituting the values of the non-linear terms and carrying out some algebraic manipulations by defining a new variable c m after expanding up to m = 2 become ( ) Neglecting the non-linear higher-order terms in Equations (5)-(7), we get ( ) It is clear that the z-axis solution, obtained by putting X = Y = 0, does not depend upon X and Y and c 2 > 0. Hence we can conclude that the motion in Z-direction is simple harmonic. The motion in XY-plane is coupled. A fourth degree polynomial is obtained which gives two real and two imaginary roots as eigenvalues: The solution of the linearized Equations (7)-(9), as derived in [36] There is a necessity to introduce frequency and amplitude terms to perform the Lindstedt-Poincaré method. The solution of the linearized equations is again written in terms of amplitudes (A x and A z ) and phases (in-plane phase, ϕ and out-of-plane phase, ψ) and the frequencies (λ and 2 c ), with an assumption that The amplitudes x A and z A are constrained by a non-linear algebraic relationship given by Richardson as We continue the perturbation analysis by assuming the solutions of the form: where ν is a small parameter which takes care of the non-linearity.

First-Order Equations
The first-order equations are obtained by taking the coefficients of the term ν . These are given as ( ) The periodic solution to the above equations is where k = k 2 .

Second-Order Equations
The coefficients are given in the Appendix.

Third-Order Equations
By collecting the terms of 3 ν , we get the third-order equations ( ) The coefficients are given in the Appendix. Thus, the third-order analytical solution is developed.

Final Approximation
The coefficients are given in the Appendix.

Numerical Computation of Halo Orbits
The method of differential correction is a powerful application of Newton's method that employs the state transition matrix (STM) to solve various boundary value problems. Differential correction method is used to determine the initial conditions of the halo orbits from the initial guess [32] [34] [40]. Taking advantage of the fact that halo orbits are symmetric about xz-plane, the initial state vector takes the form The equations of motion and state transition matrix are integrated numerically until the trajectory crosses the xz-plane again. The desired final condition is of the form It is obtained by modifying the known parameters of the initial state vector. Then the orbit will be periodic with period The state transition matrix at 1 2 t can be used to adjust the initial values of a nearby periodic orbit. Using fourth-order Runge-Kutta method, the equations of motion are integrated until y changes sign. Then the step size is reduced and the integration goes forward again. This is repeated until y becomes almost zero, and the time at this point is defined to be 1 2 t . The tolerance is considered to be of the order of 10 −12 . The orbit is considered periodic if x  and z  are zero at 1 2 t . If this is not the case, x  and z  can be reduced by correcting two of the three initial conditions and by integrating again. Figures 2-8 provide numerical as well as analytical solution of L1 and L2 halo orbits about Saturn-Satellite systems. The variation in halo orbits due to oblateness of Saturn on its satellites is studied in section 4.

Halo Orbits in the Saturn-Satellites System
Saturn, the sixth planet from the Sun, is home to a vast array of intriguing and unique satellites. It is also the largest oblate body in the solar system. Hence, studying the halo orbits about the moons of Saturn is interesting to observe the oblateness effect on the periodic orbits about the collinear points of the Saturn-Satellite systems. Christian Huygens discovered Titan, the first known moon of Saturn in 1655. Jean-Dominique      Tethys (1684). Mimas and Enceladus were both discovered by William Herschel in 1789. We have used these known moons of Saturn for our study. The mass parameter μ and the oblateness coefficient A1 of the systems under study are presented in Table 1 [17].

Saturn-Mimas
Mimas is the closest moon to Saturn which is considerably larger than other moons closer to Saturn. Mimas has an enormous crater on one side, the result of an impact that nearly split the moon apart. We have computed the initial guesses for the differential correction process provided in Table 2. The initial guess for the halo orbit computation about L1 and L2 of Saturn-Mimas system is subjected to differential correction process and then the equations of motion are numerically integrated with Adams-Bashforth-Moulton multistep method with relative tolerance of 2.5e−10 and absolute tolerance of 1.0e−14. Figure 9 and Figure 10 provide the variation of Saturn-Mimas L1 and L2 halo orbits with and without oblateness effect. We observe that the initial conditions change for different values of A1 and hence the orbits about the collinear points L1 and L2 shift. It is observed that the effect of Saturn's oblateness on the halo orbits is significant  Figure 9. Saturn-Mimas L1 halo orbit with and without oblateness. and the L1 halo orbit shifts towards Saturn by 260 km and L2 halo orbit shifts by 11.6 km. We also find that the non-dimensional time period of the halo orbits increases at L1 and decreases at L2 of the Saturn-Mimas system.

Saturn-Enceladus
Enceladus is the second closest moon next to Mimas at a distance of 237,948 km from Saturn. It displays evidence of active ice volcanism. Cassini observed warm fractures where evaporating ice evidently escapes and forms a huge cloud of water vapour over the South Pole. Similar to Saturn-Mimas system, we compute the halo orbits about the L1 and L2 of Saturn-Enceladus system with the initial guesses provided in Table 3 and the parameters taken from Table 1. We refine these conditions with differential correction method and compute halo orbits numerically. The oblateness effect on this system is lesser than on Saturn-Mimas system. Figure 11 and Figure 12 provide the L1 and L2 halo orbits with and without oblateness for Saturn-Enceladus system.
We observe that the halo orbits about L1 and L2 shift towards Saturn with oblateness by 153 km and 25 km, respectively. Similar to Saturn-Mimas system, the non-dimensional time period of the halo orbit increases at L1 and decreases at L2.   Figure 11. Saturn-Enceladus L1 halo orbit with and without oblateness.

Saturn-Tethys
Tethys is the third closest and interesting satellite of Saturn at a distance of 294,619 km.
Tethys has a huge rift zone called Ithaca Chasma that runs nearly three-quarters of the way around the moon. Telesto and Calypso are the two moons occupying the two Lagrangian points of Saturn-Tethys system. With the help of the initial guess, we compute the halo orbits about L1 and L2 of Saturn-Tethys system using the parameters of Table

Saturn-Dione
Dione is the fourth closest moon of Saturn with considerable higher mass than Tethys, and is 377,396 km from Saturn. Similar to Tethys, Helene and Poly deuces are the two other moons occupy the corresponding Lagrangian points of Dione in Saturn-Dione system. Using Table 1 and Table 5, we compute the L1 and L2 halo orbits of Saturn-Dione system. These are shown in Figure 15 and Figure 16, respectively. The L1 halo orbit moves towards Saturn by 77.15 km and L2 halo orbit moves towards it by 0.48 km with oblateness.

Saturn-Rhea
Rhea is the moon of Saturn orbiting about 527,108 km from Saturn. It experiences the oblateness effect of Saturn very less compared to the previous moons. Using Table 1 Figure 13. Saturn-Tethys L1 halo orbit with and without oblateness.   and Table 6, the L1 and L2 halo orbits of Saturn-Rhea system are computed and are plotted in Figure 17 and Figure 18 with and without oblateness. L1 halo orbit moves towards Saturn by 50.3 km and L2 halo orbit moves towards it by 0.3 km.

Saturn-Titan
Titan is the solar system's second-largest moon with 5150 km diameter after Ganymede (5362 km) of Jupiter. Titan hides its surface beneath a thick, nitrogen-rich atmosphere.
Cassini's instruments have revealed that Titan possesses many parallels to Earth-clouds, and L2 halo orbits about Saturn-Titan system under the effects of Saturn's oblateness. Table 1 and Table 7 provide the initial conditions for numerical computation of these orbits. It is noticed that the effect of Saturn's oblateness on the halo orbit is quite less. Figure 16. Saturn-Dione L2 halo orbit with and without oblateness.  Figure 17. Saturn-Rhea L1 halo orbit with and without oblateness.  Oblateness attracts L1 halo orbit towards Saturn by 2.88 km and L2 halo orbit by 0.1 km. Figure 19 and Figure 20 show Saturn-Titan L1 and L2 halo orbits. As is expected, the effect of oblateness on the halo orbits decreases as the distance between the satellites of Saturn and the planet (Saturn) increases.

Saturn-Iapetus
Iapetus has one side as bright as snow and other side as dark as black velvet, with a huge Figure 19. Saturn-Titan L1 halo orbit with and without oblateness. ridge running around most of its dark-side equator. Iapetus is 3,560,820 km from Saturn and is smaller than Titan, Rhea and larger than Mimas, Enceladus. Initial conditions for computation of L1 and L2 halo orbits of Saturn-Iapetus system are given in Table 8. Figure 21 and Figure 22 show L1 and L2 halo orbits for this system. The effect of oblateness on the halo orbits of Saturn-Iapetus system is found to be negligible.

Results and Conclusion
Halo orbits in the vicinity of L1 and L2 collinear points in Saturn-Satellites systems in the frame work of circular restricted three-body problem with more massive primary Saturn as an oblate spheroid with its equatorial plane coincident with the plane of motion of the primaries are considered. The halo orbits with oblateness effect of Saturn for seven largest satellites of Saturn are computed through the differential correction method of Mireles [32]. It is found that oblateness effect on halo orbits of the satellites closer to Saturn has significant effect compared to the satellites away from it. The halo orbits L1 and L2 are found to move towards Saturn with oblateness effect and the results are tabulated in Table 9.