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 circular- restricted problem (CRTBP) and found that an integral of motion exists. Plummer  , 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 in mid 1960s initiated an 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  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  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.  , Richardson  , Huber et al.  , Gomez et al.   , Rausch  , Nakamiya et al.  , Koon et al.  and Nath & Ramanan  . Considering the motion in the vicinity of the collinear points in CRTBP, Breakwell and Brown  numerically extended the work of Farquhar and Kamel  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   and Howell et al.  .
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   , Subba Rao and Sharma  -  and Sharma  -  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  , Papadakis  , Kalantonis et al.  , Raheem et al.  , Stuchi et al.  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  (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  ). We utilize Newton’s method of differential correction (Mireles  ; Eapen and Sharma  ; Nishanth and Sharma  ) to compute the halo orbits numerically about the Lagrangian points L1 and L2 in seven of the Saturn-Satellites systems.
2. Equations of Motion
The equations of motion for the restricted three-body problem are considered with, and as the positions of three bodies with masses (more massive oblate spheroid-Saturn), (Satellite of Saturn) and (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 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..
The parameter defines the position of the larger and smaller primaries as and, respectively (Szebehely  ; Sharma  ; Bhatnagar and Chawla  ).
The three-dimensional equations of motion are:
where is the pseudo potential of the system and is given as
In the above expression and are the position vectors from the more massive and smaller primaries to the particle, respectively.
Figure 1. Diagram of the circular restricted three-body problem in normalized units
3. 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. or. 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  by considering the more massive primary as an oblate spheroid with the help of Lindstedt-Poincaré method.
3.1. Analytical Approximation
For obtaining an analytical solution, following Tiwary and Kushvah  , the origin is transferred to the Lagrangian points L1 and L2 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.   and  .
The non-linear terms are expanded by using the following formula as given by  :
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 cm 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 c2 > 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  , is
are arbitrary constants. Since we are concentrating on constructing a halo orbit, which is periodic, we consider.
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 (Ax and Az) and phases (in-plane phase, ϕ and out-of-plane phase, ψ) and the frequencies (λ and), with an assumption that, as
The amplitudes and are constrained by a non-linear algebraic relationship given by Richardson as
where l1 and l2 depend upon the roots of the characteristic equation of the linear equation. The correction term, arises due to the addition of frequency term in Equation (10).
Hence, any halo orbit can be characterized by specifying a particular out-of-plane amplitude of the solution to linearized equations of motion. Both analytical and numerical developments employ this scheme. From the above expression, we can find the minimum permissible value of Ax to form the halo orbit (Az > 0).
The phases ϕ and ψ are related as
When Ax is greater than certain value, the third-order solution bifurcates. This bifurcation is manifested through the phase-angle constraint. The solution branches are obtained according to the value of m. For m = 1, Az is positive and we have the northern halo (z > 0) and for m = 3, Az is negative and we have the southern halo (z < 0).
Lindstedt-Poincaré method involves successive adjustments of the frequencies to avoid secular terms and allows one to obtain approximate periodic solution. The equations of motion with non-linear terms up to third-order approximation as in Richardson and Thurman and Worfolk are
A new independent variable τ = ωt is introduced, where refers to.
It is to be noted that the values of are chosen in such a way that the secular terms get removed with successive approximations. Most of the secular terms are removed by the following assumptions
The coefficients and are given in Appendix. The equations of motion are:
We continue the perturbation analysis by assuming the solutions of the form:
where is a small parameter which takes care of the non-linearity.
Substituting Equations (14)-(16) in Equations (11)-(13), and equating the coefficients of the same order of, we get the first-, second- and the third-order equations, respectively.
3.1.1. 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 = k2.
3.1.2. Second-Order Equations
By collecting the terms of, we get the second-order equations
To remove the secular terms, we need to set the value ω1 = 0. The particular solution of the second-order equations are obtained with the help of Maxima software as
The coefficients are given in the Appendix.
3.1.3. Third-Order Equations
By collecting the terms of, we get the third-order equations
From the above equations, it is not possible to remove the secular terms by setting ω2 = 0. Hence, the phase relation is used here to remove the secular terms.
The solution of the third-order equation is obtained as
The coefficients are given in the Appendix. Thus, the third-order analytical solution is developed.
3.1.4. Final Approximation
The mapping, and, will remove ν from all the equations. We now combine the solutions up to third-order to get the final solution:
The coefficients are given in the Appendix.
3.2. 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    . 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 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. The tolerance is considered to be of the order of 10−12. The orbit is considered periodic if and are zero at. If this is not the case, and 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.
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
Figure 2. Numerical and analytical solutions for Saturn-Mimas system.
Figure 3. Numerical and analytical solutions for Saturn-Enceladus system.
Figure 4. Numerical and analytical solutions for Saturn-Tethys system.
Figure 5. Numerical and analytical solutions for Saturn-Dione system.
Figure 6. Numerical and analytical solutions for Saturn-Rhea system.
Figure 7. Numerical and analytical solutions for Saturn-Titan system.
Cassini made the next four discoveries: Iapetus (1671), Rhea (1672), Dione (1684), and 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  .
Figure 8. Numerical and analytical solutions for Saturn-Iapetus system.
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
Table 1. Systems and their parameters.
Table 2. Initial conditions for Saturn-Mimas system.
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.
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-dimen- sional time period of the halo orbit increases at L1 and decreases at L2.
Figure 10. Saturn-Mimas L2 halo orbit with and without oblateness.
Table 3. Initial conditions for Saturn-Enceladus system.
Figure 11. Saturn-Enceladus L1 halo orbit with and without oblateness.
Figure 12. Saturn-Enceladus L2 halo orbit with and without oblateness.
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 1 and Table 4. Figure 13 and Figure 14 illustrate the L1 and L2 halo orbits. We observe that, the initial conditions change for different values of A1 and the halo orbits shift towards Saturn by 101 km at L1 and 3 km at L2.
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.
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
Table 4. Initial conditions for Saturn-Tethys system.
Figure 13. Saturn-Tethys L1 halo orbit with and without oblateness.
Figure 14. Saturn-Tethys L2 halo orbit with and without oblateness.
Figure 15. Saturn-Dione L1 halo orbit with and without oblateness.
Table 5. Initial conditions for Saturn-Dione system.
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.
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, dunes, mountains, lakes, and rivers. Titan’s atmosphere is approximately 95 percent nitrogen with traces of methane. While Earth’s atmosphere extends about 60 km into space, Titan’s extends nearly by 600 km (10 times that of Earth’s atmosphere) into space. Effect of oblateness of Saturn on Titan in the frame work of restricted three-body problem was previously studied by Beevi and Sharma   . Here we compute L1 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.
Table 6. Initial conditions for Saturn-Rhea system.
Figure 17. Saturn-Rhea L1 halo orbit with and without oblateness.
Figure 18. Saturn-Rhea L2 halo orbit with and without oblateness.
Table 7. Initial conditions for Saturn-Titan system.
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.
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.
Figure 20. Saturn-Titan L2 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.
5. 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  . 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.
Table 8. Initial conditions for Saturn-Iapetus system.
Figure 21. Saturn-Iapetus L1 halo orbit with and without oblateness.
Figure 22. Saturn-Iapetus L2 halo orbit with and without oblateness.
Table 9. Variation in location and time period of halo orbit with oblateness effect.
Coefficients for the second and third-order equations and solution:
Submit or recommend next manuscript to SCIRP and we will provide best service for you:
Accepting pre-submission inquiries through Email, Facebook, LinkedIn, Twitter, etc.
A wide selection of journals (inclusive of 9 subjects, more than 200 journals)
Providing 24-hour high-quality service
User-friendly online submission system
Fair and swift peer-review system
Efficient typesetting and proofreading procedure
Display of the result of downloads and visits, as well as the number of cited articles
Maximum dissemination of your research work
Submit your manuscript at: http://papersubmission.scirp.org/
Or contact firstname.lastname@example.org