Halo Orbits have been the topic of interest for the past few decades among the research society. It is a periodic, three-dimensional orbit around the Lagrangian points in a three-body problem. The name “Halo” was first used by Farquhar  in his doctoral thesis (1968). It was first proposed for the Apollo mission for placing the satellite in an orbit around L2 Lagrangian point in the Earth-Moon system. An enormous contribution was made by Richardson  towards the construction of Halo orbits using Lindstedt-Poincare method. His method explained the general idea about Halo orbits. The first and simplest periodic exact solution to the three-body problem is the motion on collinear ellipses found by Euler . Euler studied the motion of the Moon assuming the Earth and the Sun orbit each other on circular orbits and the Moon to be massless. This approach is known as the restricted three-body problem. The Lagrangian points are an equilibrium solution to the restricted three-body problem. Szebehely  has an excellent thesis on the Lagrangian points and the restricted three-body problem. The classical model of the three-bodies is not accurate for studying any system in the solar family as it does not consider the perturbing forces such as oblateness and radiation pressure as given by Dutt and Sharma . Some of the significant works in the Photogravitational Restricted Three-Body Problem were done by Radzievskii , Bhatnagar and Chawla , Eapen and Sharma . The first person to study about the solar radiation pressure was Radzievskii and he proposed that the maximum force due to radiation pressure that acts in the radial direction can be given by
where q is a pressure reduction factor which is defined in terms of particle radius a, density δ and radiation pressure efficiency factor x as
is a variable, depending upon the specification of the third body which is also a radiating body.
A star radiates or shines due to nuclear fusion of hydrogen into helium in its core which then releases energy that traverses through the interior of the star and then radiates into space. The eventual existence of star is directly dependent on its mass. The characteristics of a star like its diameter and temperature change over time while its rotation and movement are dependent on its environment. When two or more stars that are bound by gravity and move around each other in a stable orbit in the same plane is called a binary or multi-star system.
This paper considers that both the primaries as source of radiation and the following calculations and results are made. So, there will be two pressure radiation terms viz. q1 and q2 representing the radiation of the larger and smaller primaries, respectively.
Lindstedt-Poincare method is used to find the solution for the equations of motion of CRTBP through successive approximations. Numerical methods are used with this analytical solution as the initial value to generate the required Halo orbits. Changes in the Halo orbits for various values of radiation pressure are also analyzed.
2. Circular Restricted Three-Body Problem
In the Circular Restricted Three-Body Problem, two bodies revolve around their center of masses in circular orbits under the influence of their mutual gravitational attraction and a third body moves in a plane defined by the motion of the two revolving bodies. The problem of the motion of the third body is called the Circular Restricted Three-Body Problem, henceforth referred to as CRTBP. In the CRTBP, the mass of the third body is considered negligible compared to the primaries and its presence does not disturb the circular motion of the primaries (Figure 1).
Figure 1. Circular restricted three-body problem in normalized units.
The restricted three-body problem has five equilibrium points, called Lagrangian or Libration Points. These points are the points of zero velocities and the objects placed in these points remains stable. Out of the five Lagrangian points, three are collinear (L1, L2, L3) and the other two points (L4, L5) form the equilateral triangle. Although the Lagrangian point is just a point in empty space, its peculiar characteristic is that it can be orbited.
The circular restricted three-body problem (CRTBP) is a special case of the restricted three-body problem where bigger and smaller primaries move in circular motion around their common center of mass. The five equilibrium points of CRTBP are known as the Lagrangian points where the gravitational forces due to two primaries and the centrifugal force on a spacecraft are balanced.
Another assumption is made that the universe is strictly limited to these three bodies under consideration, so there is no other gravitational influence in the system. Though the limitations that have been introduced seem to restrict the problem to the point of impracticality, there are many situations in the solar system to which the CRTBP applies.
3. Equations of Motion
The analytical approximation of three dimensional, periodic orbits about collinear points is performed. The equations of motion are developed from a Lagrangian approach described in Richardson . The analytical approximation solution for the equations of motion is calculated through computer-based program.
The equations of motion for the RTBP is written as (Szebehely , Sharma , Simmons  )
Here the oblateness is neglected so and q1 and q2 are the radiation pressure terms of larger and smaller primary,
where m1 and m2 are masses of larger and smaller primary respectively.
Here and .
Location of Lagrangian Points
These are the equations for the location of the Lagrangian points including the perturbing force of radiation pressure of both the primaries. The equations for the location of Lagrangian points have been referred from Nishanth and Sharma .
The location of Lagrangian points varies with the variation of radiation pressures. This can be observed after solving the equations. The equations will give the positions of Lagrangian points.
4. Analytical Computation of Halo Orbits
For the computation of the Halo orbits, the origin should be transferred to the Lagrangian points L1 and L2. The transformation is given by
Hence the equation of motion can be again written as
upper sign in the above equation depict L1 and lower sign depict L2. The distance between smaller primary and larger primary is considered to be unity.
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 as given by Koon .
where , .
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 steps by defining a new variable, cm after expanding up to m = 2, we get
Neglecting the non-linear higher order terms in Equations (3.1)-(3.3), 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. Meanwhile, the motion in XY-plane is coupled.
A fourth-degree polynomial is obtained which gives two real and two imaginary roots as Eigen values.
Thus, the solution to the linearized Equations (3.4)-(3.6) will be of the given form as derived by Thurman and Worfolk  as
Here are arbitrary constants. Since we are interested in periodic solutions we take .
As mentioned earlier, there is a necessity to introduce frequency and amplitude terms to perform the Lindstedt-Poincare method. Then the solution to the linear problem can be written in terms of amplitudes, phases and the frequencies (λ and ), as
Here an assumption is made that .
4.1. Amplitude Constraint
When we consider the periodic 3D orbits, the amplitudes Ax and Az are related by non-linear algebraic relationship given by Richardson .
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 (3.6).
4.2. Phase Constraint
For periodic 3D orbits, the in-plane phase (ϕ) and out-of-plane phase (ψ) 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 northern Halo (z < 0).
4.3. Lindstedt-Poincare Method
We can find better approximations to the non-linear problem in a neighborhood of the equilibrium point solutions by using the perturbation techniques of Lindstedt-Poincare. 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 as discussed in Equation (5.1). Here on, refers to .
We choose the value of such that the secular terms get removed with successive approximations.
Under the following assumptions and , most of the secular terms are removed.
The coefficients and are given in Appendix. The equations of motion changes again.
We continue the perturbation analysis by assuming the solutions of the form, where is a small parameter which takes care of the non-linearity.
We substitute these Equations (3.10)-(3.12) in Equations (3.7)-(3.9) and equating the coefficients of the same order of , we get the first, second and third order equations respectively.
4.3.1. First-Order Equations
The analytical solution can be found using mathematical software like Maxima and Mathematica. The first-order equation is obtained by taking the coefficients of the term . It is given as
The periodic solution to the equations above is
where, k = k2.
4.3.2. Second-Order Equations
To get the second-order equations we segregate the terms of ν2,
We need to set the value , to remove the secular terms. The particular solution of the second-order equations is found using mathematical software like Maxima and Mathematica.
The coefficients are given in Appendix I.
4.3.3. Third-Order Equations
To get the third-order equations we segregate the terms of ν3,
From the above equations, it is not possible to remove the secular terms by setting . Hence, the phase relation is used here to remove the secular terms.
Hence, the solution of the third-order equation will be
The coefficients are given in Appendix I. Thus, the third-order analytical solution is developed.
4.3.4. Final Approximation
Applying the mapping, and , will remove from all
equations. We now combine the solutions up to third order to get the final solution.
The coefficients are given in the Appendix I.
5. Analytical Construction of Halo Orbits
The output from the analytical computation of Halo orbit is used to construct Halo orbits through computer-based codes and is plotted using software such as MATLAB. Input parameters are given such as mass ratio (µ), mean distance between the two primaries, radiation pressures (q1 and q2), amplitude in Z-direction (Az).
Halo orbits are plotted for the classical case of Antares-Sun-Proxima Centauri system around L1 Lagrangian point. µ = 0.0011314, q1 = q2 = 1, Az = 110,000 km (Figures 2-5).
Figure 2. X vs Y—classical.
Figure 3. X vs Z—classical.
Figure 4. Y vs Z—classical.
Figure 5. 3D halo orbit—classical.
6. Numerical Computation
Differential Correction (DC) schemes use the STM for targeting purposes. One application is the iterative process to isolate a trajectory arc that connects two points in solution space. Differential Correction techniques can be used to quickly obtain a solution with the desired parameters in a wide range of problems. A sufficiently accurate guess for the initial state is always required. Differential corrections are often used to obtain periodic solutions to the non-linear differential equations in the CRTBP. A common assumption, making use of the symmetry in this problem, is that the desired solution is symmetric about the x-z plane.
The initial guess for finding the solution is taken from the analytic solution. Since the halo orbits are symmetric about x-z plane, (y = 0), and they intersect this plane perpendicularly i.e. ( ) , the initial state vector takes the form,
The final state vector which again lies on the same plane, takes the form as given below and crosses the x-z plane perpendicularly.
Then the orbit will be periodic with period . The transition matrix at can be used to adjust the initial values of a nearby periodic orbit. Using Runge-Kutta fourth order 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 as . The orbit is considered periodic if and are nearly zero at . If this is not the case, and can be reduced by correcting two of the three initial conditions and integrate again.
6.1. Numerical Computation of Halo Orbits
The following plots are 2D and 3D plots of the halo orbits for perturbed models. The variation in size, shape of the halo orbits is observed by comparing it with the plots of classical case and the effects of radiation pressures on the halo orbit is observed. The numerical computation of halo orbits is referred from Chidambararaj P. and Sharma R. K. .
The input parameters are µ = 0.0011314, q1 = 0.98, q2 = 0.97 (Figures 6-10).
Figure 6. X vs Y.
Figure 7. Y vs Z.
Figure 8. X vs Z.
Figure 9. 3D halo orbit (q1 > q2).
Figure 10. 3D halo orbit (q2 > q1).
6.2. Perturbing Effects on the Orbit
The effects due to radiation pressure can be seen from Figure 11 below and it is seen that with an increase in radiation pressure, the orbit moves towards the larger source of radiation irrespective of the size of the primary. The oblateness coefficient are taken as A2 = 0, μ = 0.0011314 and the two cases of radiation pressures are studied when q1 = 0.98, q2 = 0.97 and when q1 = 0.96, q2 = 0.97. With increasing perturbing forces, it is seen that the orbit will move towards the primary with larger perturbing force and thus it doesn’t remain Halo orbit. Hence it is to be seen that the value of Az should be taken care such that the mission objectives are optimally obtained (Table 1).
It has been observed that the orbit length shortens and the width increases closer to the primary whose radiation pressure is larger. The size factor also comes into effect on the Halo orbit. This can be observed from the 3D plots above as both the plots are not mirror images of each other. This tells that the orbit doesn’t vary with the radiation pressure only. The variation of the orbit shape and size is affected by the mass of the primary as we increase the radiation pressure.
Case: When the radiation pressure for both the primaries is equal and not equal to 1.
Around the L1 Lagrangian point for the Sun-Proxima system.
For µ = 0.109509, q1 = q2 = 0.98 Halo period = 2.4662.
Figure 11. 3D halo orbit.
Table 1. Radiation pressure parameters, mass ratio and halo period.
7. Summary and Conclusions
The radiation pressures of larger as well as smaller primary are added. The Lindstedt-Poincare method gives approximate good initial values for the numerical solution. The numerical solution makes use of the adaptive Runge-Kutta fourth-order method as integrator. In the vicinity of L1 point, it has been seen that with increase in radiation pressure, the Halo orbit increases in orbital period. It has been seen that the Halo orbits are very sensitive to initial conditions.
Coefficients for the second and third-order equations and solution