Halo Orbits around Sun-Earth L1 in Photogravitational Restricted Three-Body Problem with Oblateness of Smaller Primary

This paper deals with generation of halo orbits in the three-dimensional photogravitational restricted three-body problem, where the more massive primary is considered as the source of radiation and the smaller primary is an oblate spheroid with its equatorial plane coincident with the plane of motion. Both the terms due to oblateness of the smaller primary are considered. Numerical as well as analytical solutions are obtained around the Lagrangian point L1, which lies between the primaries, of the Sun-Earth system. A comparison with the real time flight data of SOHO mission is made. Inclusion of oblateness of the smaller primary can improve the accuracy. Due to the effect of radiation pressure and oblateness, the size and the orbital period of the halo orbit around L1 are found to increase.


Introduction
The subject of halo orbits in restricted three-body problem (RTBP) has received considerable attention in the last four decades.A halo orbit is a periodic, three-dimensional orbitnear the collinear Lagrangian points in the three-body problem.Analytical solutions to generate the halo orbits around any of the collinear points can be obtained.Farquhar first used the name "halo" in his doctoral thesis [1] and proposed an idea for placing a satellite in an orbit around L 2 of the Earth-Moon RTBP.Had this idea been implemented, one could continuously view the Earth and the dark side of the moon at the same time.A communication link satellite is not necessary, if satellites are placed around Lagrangian points of the respective systems.Breakwell & Brown [2] generated halo orbits for the Earth-Moon system.Richardson [3] contributed significantly to obtain the halo orbits, by finding an analytical solution using Lindstedt-Poincaré method.A good amount of research has been carried out till date in this interesting area of halo orbits.Some of the important contributions are by Howell [4], Howell and Pernicka [5], Folta and Richon [6] and Howell et al. [7].
The classical model of the restricted three-body problem are relatively less accurate for studying the motion of any Sun-planet system as they do not account for the effect of the perturbing forces such as oblateness of planets, cosmic rays, magnetic field, radiation pressure etc. Cosmic rays are immensely high-energy radiation, mainly originating outside the solar system.They produce showers of secondary particles that penetrate and impact the Earth's atmosphere and sometimes even reach the surface.They interact with gaseous and other matter at high altitudes and produce secondary radiation.The combination of both contributes to the space radiation environment.In addition to the particles originating from the Sun are particles from other stars and heavy ion sources such as nova and supernova in our galaxy and beyond.In interplanetary space these ionizing particles constitute the major radiation threat.These particles are influenced by planetary or Earth's magnetic field to form radiation belts, which in Earth's case are known as Van Allen Radiation belts, containing trapped electrons in the outer belt and protons in the inner belt.The composition and intensity of the radiation varies significantly with the trajectory of a space vehicle.Anomalies in communication satellite operation have been caused by the unexpected triggering of digital circuits by the cosmic rays.In our present study, we restrict ourselves with the solar radiation and oblateness of the planet.
Radzievskii [8] was the first one to study the effect of solar radiation pressure.He found out that the maximum force due to the radiation pressure acts in the radial direction, given by ( ) where q is defined in terms of particle radius a, density δ and radiation pressure efficiency factor x as (c.g.s.units) ε is a variable, depending upon the nature of the third body (satellite).The value of q can be considered as a constant, if the fluctuations in the beam of solar radiation and the effect of planet's shadow are neglected.Using the model of [1], Dutt and Sharma [9] studied periodic orbits in the Sun-Mars system using the numerical technique of Poincare surface of sections and found out more than 74 periodic orbits.
Sharma and Subba Rao [10] introduced the oblateness of the more massive primary in the three-dimensional restricted three-body problem.It has two terms, one with a z term in the numerator.Sharma [11] studied the periodic orbits around the Lagrangian points in the planar RTBP by considering Sun as source of radiation and smaller primary as an oblate spheroid with its equatorial plane coincident with the plane of motion.Tiwary and Kushvah [12] followed the model of Sharma [11] to study the halo orbits around the Lagrangian points L 1 and L 2 analytically.However, they did not consider the z term in oblateness in their study.In the present work, we have considered both the terms due to oblateness of the smaller primary in the photogravitational restricted three-body problem to study the halo orbits analytically as well as numerically around L 1 .
The first satellite placed in the halo orbit at Sun-Earth L 1 point was International Sun-Earth Explorer-3 (ISEE-3), launched in 1978.Solar Heliospheric Observatory (SOHO), launched in 1975 by NASA succeeded ISEE-3.We have taken data of the path of the SOHO mission from the mission website over a period of January-June 2008, to validate our analytical and numerical solutions.

Circular Restricted Three-Body Problem
The circular RTBP consists of two primary masses revolving in circular orbits around their centre of mass under the influence of their mutual gravity.The third body of infinitesimal mass moves under the gravitational effect of these two primaries (Figure 1).RTBP has five equilibrium points, called Lagrangian or libration points.These points are the points of zero velocities and an object placed in these points remains there.Out of the five Lagrangian points, three are collinear (L 1 , L 2 , L 3 ) and the other two points (L 4 , L 5 ) form equilateral triangles with the primaries.Although the Lagrangian point is just a point in empty space, its peculiar characteristic is that it can be orbited.Halo orbits are three-dimensional orbits around the collinear points.
The equations of motion for the RTBP (Szebehely [13]) including radiation pressure and oblateness of the smaller primary is written in accordance with Sharma & Subba Rao [10], Sharma [11] as 2 , x x ny x . where , where m 1 and m 2 are masses of larger and smaller primary respectively.The perturbed mean motion, n of the primaries due to oblateness is given by . 5 AE, AP being the dimensional equatorial and polar radii of the smaller primary and R is the distance between the primaries.The two terms occurring in Ω due to oblateness of the smaller primary were introduced by Sharma and Subba Rao [10].

Computation of Halo Orbits
For the computation of the halo orbits, the origin is transferred to the Lagrangian points L 1 and L 2 .The trans-formation is given by 1, The equation of motion can be written as The upper sign in the above equations depicts the Lagrangian point L 1 and the lower sign corresponds to L 2 .The distance between these Lagrangian points and the smaller primary is considered to be the normalized unit as in Koon et al. [14] and Tiwary and Kushvah [12].
The usage of Legendre polynomials can result in some computational advantages, when non-linear terms are considered.The non-linear terms are expanded by using the following formula given in Koon et al. [14]: 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 by some algebraic manipulations by defining a new variable c m after expanding up to m = 2 become ( ) ( ) , with Neglecting the higher-order terms in Equations (( 4)-( 6)), 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 in Thurman and Worfolk [15], is are arbitrary constants.Since we are concentrating on constructing a halo orbit, which is periodic, we consider 1 2 0 a a A A = = .Frequency and amplitude terms are introduced to find solution through Lindstedt-Poincaré method.The solution of the linearized equations is written in terms of the amplitudes (A x and A z ) and the phases (In-plane phase, ϕ and out-of-plane phase, ψ) and the frequencies (λ and 2 c ), with an assumption that

Amplitude Constraint
For halo orbits, the amplitudes x A and z A are constrained by a non-linear algebraic relationship given by Richardson [3]:  9).Hence, any halo orbit can be characterized by specifying a particular out-of-plane amplitude z A of the solu- tion to linearized equations of motion.Both analytical and numerical methods employ this scheme.From the above expression, we can find the minimum permissible value of A x to form the halo orbit (A z > 0).

Phase Constraint
For halo orbits, the phases ϕ and ψ are related as When A x 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, A z is positive and we have the northern halo (z > 0).For m = 3, A z is negative and we have the southern halo (z < 0).

Lindstedt-Poincaré Method
Lindstedt-Poincaré method involves successive adjustments of the frequencies to avoid secular terms and provides approximate periodic solution.The equations of motion with non-linear terms up to third-order approximation as in Richardson [3] and Thurman and Worfolk [15] are The values of n ω are chosen in such a way that the secular terms are removed with successive approxima- tions.Most of the secular terms are removed by the following assumptions: The coefficients 1 s and 2 s 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.

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

Second-Order Equations
Coefficients of the term 2 ν provide the second-order equations ( ) ; .
To remove the secular terms, we set ω 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.

Third-Order Equations
By collecting the coefficients of 3 ν , we get the third-order equations ( )

= + =
The solution of the third-order equation is obtained as ( ) ( ) The coefficients are given in the Appendix.Thus, the third-order analytical solution is obtained.

Final Approximation
The mapping, , will remove ν from all the equations.We now combine the so- lutions up to third-order to get the final solution as ( ) ( ) The coefficients are given in the Appendix.

Input Parameters
The input parameters to generate a halo orbit include mass ratio (μ), mean distance between the two primaries, radiation pressure (q or ε), oblateness coefficient (A 2 ), amplitude in z-direction (A z ).Once all the input parameters are given, the co-ordinates for the halo orbit are computed.

Halo Orbit for Classical Case
The halo orbits are generated with the following data for Sun-Earth L 1 : Mass ratio = 0.000003, mean distance = 149,600,000 km, q = 1, A 2 = 0.The amplitude in z-direction is taken to be A z = 110,000 km (Amplitude of the ISEE-3 mission around the Sun-Earth L 1 point).Figures 2-5 show the different views of the halo orbits for the classical case.

Halo Orbit for Different Values of Radiation Pressure
The halo orbits are generated for different values of radiation pressurein the vicinity of L 1 .Figure 6 shows the variation in orbit for different values of radiation pressure.The orbital period increases and the orbits increase in size with the increase in the radiation pressure of the more massive primary.
Table 1 shows the change in the distance of the Lagrangian point L 1 due to radiation pressure.It is noticed that the Lagrangian point L 1 moves towards the more massive primary with increase in radiation pressure.It is also noticed that L 1 moves further towards the more massive primary for q < 0.8.

Numerical Computation of Halo Orbits
The third-order analytical solution provides a good initial estimate to compute the halo orbits numerically.Analytical approximation must be combined with a suitable numerical method to obtain accurate halo orbits around the Lagrangian points.The method of differential correction is a powerful application of Newton's method that employs the State Transition Matrix (STM) to solve various kinds of boundary value problems.In order to propagate an orbit in RTBP, the equations of motion are numerically integrated by using adaptive fourth-order Runge-Kutta method, until the desired orbit is obtained.
The initial guess for finding the numerical solution is taken from the analytic solution.Since the halo orbits are symmetric about xz-plane (y = 0), and they intersect this plane perpendicularly i.e. ( ) , the initial state vector takes the form ( ) ( ) , 0, , 0, , 0 .X t x z y =  The final state vector which lies on the same plane, takes the form as given below and crosses the xz-plane perpendicularly:   , 0, , 0, , 0 .
 Then the orbit will be periodic with period 1 2

. T t =
The state transition matrix at 1 2 t can be used to adjust the initial values of a nearby periodic orbit.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 around 10 −12 .The orbit is considered periodic if x  and z  are nearly 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 then integrate the equations of motion again.

Numerically Generated Halo Orbits
The numerically generated orbits which are more accurate, are usually considered for mission purposes.shows numerically generated orbit due to variation in radiation pressure.The numerically generated obits show a close resemblance to the analytically generated orbits.Figure 8 shows the comparison between them.
The analytical solution which includes the radiation pressure of more massive primary and oblateness of the smaller primary can improve the accuracy over the existing solution as given by Tiwary and Kushvah [12].
The time period of SOHO orbit was around 179 days and its dimensions were approximately 4.3 × 2.7 × 3.7 m.The computed value of radiation pressure is roughly 0.99997.The flight data of SOHO mission from January to June 2008 was taken from the mission website and is plotted.From Figure 9 and Figure 10, it is clear that the addition of the extra term due to oblateness shows improvement in the orbit computation.However, the deviation in the orbit is due to the actual scenario, where the real satellite in space undergoes the effect of various other perturbations.

Improvement over the Existing Results
It can be seen from the Figure 11 that the results show improvement in orbit computation due to the addition of the new term due to oblateness.
The time period of the orbit increases by 30 minutes for q = 1 and A 2 = 0.000001 as given in Table 2. Thus, it is clear that the addition of the new term in the potential function has helped in predicting the orbit closer to the real time condition.

Perturbing Effects on the Orbit
The effects due to radiation pressure and oblateness can be seen from Figure 12.It is seen that with increase in radiation pressure as well as oblateness, the orbit moves towards the source of radiation i.e. the more massive primary.The oblateness coefficient are taken as A 2 = 0, 0.000001, 0.000002 and the radiation pressure, q = 1, 0.99, 0.98.With increase in perturbing forces, it is seen that the orbit moves towards the more massive primary.
Since the amplitude A z is bounded by amplitude constraint, the variation in it also affects the size and time period of the orbit, which can be seen in Figure 13.The time period of the orbit also changes with the perturbations and the amplitude A z .For the halo around L 1 , with increase in radiation pressure and oblateness, the time period increases.It can be noticed in Figure 14 and Table 3.
The amplitudes about the z-axis, A z was taken in a range of 80,000 to 140,000 km.There is a minimum A x for the feasibility of a halo orbit, when we fix the amplitude A z .

Conclusion
The analytical solution obtained using Lindstedt-Poincaré method provides good initial estimate to generate halo orbits around the collinear point L 1 with numerical integration in the photogravitational restricted three-body    problem, when the smaller primary is considered as an oblate spheroid with its equatorial plane coincident with the plane of motion.A comparison with the orbital data of SOHO mission is carried out.The radiation pressure and oblateness are found to play significant role in improving the accuracy of halo orbits computation.In the vicinity of L 1 , it is noticed that with increase in radiation pressure, the size of halo orbit increases with the increase in the orbital period.The results for L 1 coincide with the results obtained by Eapen and Sharma [16].The Lagrangian point L 1 moves towards the more massive primary with perturbations and the orbit tends to bend towards the smaller primary with the increase in oblateness.The future aspects of the project involve designing
where l 1 and l 2 depend upon the roots of the characteristic equation of the linear equation.The correction term the addition of frequency term in Equation (

Figure 2 .
Figure 2. X vs Y for Sun-Earth system.

Figure 3 .
Figure 3. Y vs Z for Sun-Earth system.

Figure 7 Figure 4 .
Figure 4. X vs Z for sun-earth system.

Figure 8 .
Figure 8.Comparison of numerically and analytically generated halo orbits.

Figure 9 .
Figure 9.Comparison of numerically generated halo orbit with SOHO Mission orbit.

Figure 10 .
Figure 10.Comparison of numerically generated halo orbit with SOHO mission orbit (enlarged view).

Figure 11 .
Figure 11.Effect of adding the new term on halo orbit at L 1 with q = 1, A 2 = 0.000001.
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/ From the above equations, it is not possible to remove the secular terms by setting ω 2 = 0. Hence, the phase relation is used to remove the secular terms.

Table 1 .
Variation in the location of Lagrangian point L 1 with radiation pressure q.

Table 2 .
Effect of new term on time period of the halo orbit.
q ε A2 Time period without new term (days) Time period considering new term (days) Difference (days)

Table 3 .
Time period of the halo orbit around L 1 (μ = 0.000003).Earth to Mars with the help of halo orbits.These trajectories can reduce the cost of the interplanetary and deep space missions.