In this paper, a third-order analytic approximation is described for computing the three-dimensional periodic halo orbits near the collinear L₁ and L₂ Lagrangian points for the photo gravitational circular restricted three-body problem in the Sun-Earth system. The constructed third-order approximation is chosen as a starting initial guess for the numerical computation using the differential correction method. The effect of the solar radiation pressure on the location of two collinear Lagrangian points and on the shape of the halo orbits is discussed. It is found that the time period of the halo orbit increases whereas the Jacobi constant decreases around both the collinear points taking into account the solar radiation pressure of the Sun for the fixed out-of-plane amplitude.

Keywords: photogravitational CRTBP, halo orbits, lindstedt-poincare method, Newton’s method, radiation pressure

From the last few decades, the space science community has shown considerable interest in missions which take place in the vicinity of the Lagrangian points in the restricted three-body problem (RTBP) of the Sun-Earth and the Earth-Moon systems.¹ Designing trajectories for these missions is a challenging task due to inadequacy of the conic approximations. The RTBP deals the situation where one of the three bodies has a negligible mass, and moves under the gravitational influence of two other bodies.^2–9 In the RTBP, the circular restricted three-body problem (CRTBP) is a special case where two massive bodies move in the circular motion around their common centre of mass.^10–14 The collinear Lagrangian point orbits have paid a lot of attentions for the mission design and transfer of trajectories.^15–22 When the frequencies of two oscillations are commensurable, the motion becomes periodic and such an orbit in the three-dimensional space is called halo.²³ Lyapunov orbits are the two-dimensional planar periodic orbits. These planar periodic orbits are not suitable for space applications since they do not allow the out-of-plane motion, e.g., a spacecraft placed in the Sun-Earth L₁ point must have an out-of-plane amplitude in order to avoid the solar exclusion zone (dangerous for the downlink); a space telescope around the Sun-Earth L₂ point must avoid the eclipses and hence requires a three-dimensional periodic orbit. Since the RTBP does not have any analytic solution, the periodic orbits are difficult to obtain because the problem is highly nonlinear and small changes in the initial conditions break the periodicity.²⁴ Farquhar²³ was the first person who introduced analytic computation of the halo orbit in his PhD thesis. In 1980,²⁵ introduced a third-order analytic approximation of the halo orbits near the collinear libration points in the classical CRTBP for the Sun-Earth system. Thurman & Worfolk¹ and Koon et al.,²⁶ found the halo orbits for the CRTBP with the Sun-Earth system in the absence of any perturbative force using Richardson method²⁵ up to third-order. Breakwell & Brown²⁷ and Howell²⁸ numerically obtained the halo orbits in the classical CRTBP Earth-Moon system using the single step differential correction scheme. Numerous applications of the halo orbits in the scientific mission design can be seen such as investigations concerning solar exploration and helio-spheric effects on planetary environments using the spacecraft placed in these orbits at different phases. ISEE-3 was the first mission in a halo orbit of the Sun-Earth system around L₁ to study the interaction between the Earth’s magnetic field and solar wind.²⁹ Solar and Heliospheric Observatory (SOHO) mission was second libration point mission launched jointly by ESA and NASA in a halo orbit around the Sun-Earth L₁ point, still operational till date, which was a virtual carbon copy of ISEE-3’s orbit.³⁰

The classical model of the RTBP does not account perturbing forces such as oblateness, radiation pressure and variations of masses of the primaries. The photogravitational RTBP arises from the classical RTBP if at least one of the bodies is an intense emitter of radiation. Radzievskii³¹ derived the photogravitational RTBP and discussed it for three specific bodies: the Sun, a planet and a dust particle. Recently, Eapen & Sharma³² discussed the planar photogravitational CRTBP including solar radiation pressure in the Sun-Mars system around L₁ point using the initial guess of the classical CRTBP.

Numerical computation of the periodic orbits requires an initial approximation to the orbit as an approximate analytic solution. However, the analytic solutions that are available do not generally include solar radiation pressure and other perturbing forces. Including these perturbed forces in the analytic approximation increases accuracy of the approximation and therefore, simplifies the numerical computations.³³ In this paper, we discuss analytic as well as numerical computations of the halo orbits around the libration points L₁ and L₂ in the CRTBP including solar radiation pressure of the Sun. The paper is organized as follows: Section 2 deals with the governing equations of motion considering the Sun as a radiating source. Section 3 describes the construction of a third-order analytic approximate solution for the periodic orbit using the Lindstedt-Poincare method. Section 4 illustrates numerical computation of the halo orbit using Newton’s method of differential correction. Results and discussion are given in Section 5 while Section 6 concludes our study.

We suppose that the CRTBP consists of the Sun, the Earth and the Moon, and an infinitesimal body such as a spacecraft having masses m₁, m₂ and m, respectively. Here the Earth and the Moon are clubbed as a single entity and we say this as the Earth. The spacecraft moves under the gravitational influence of the Sun and the Earth (Figure 1). The Sun is assumed as the radiating body contributing solar radiation pressure.

Figure 1 Geometry of the problem in the Sun-Earth system.

Let $\left(x,y,z\right),\left(-\mu ,0,0\right)$ and $\left(1-\mu ,0,0\right)$ denote the coordinates of the spacecraft, the Sun, and the Earth, respectively, where μ is the mass ratio parameter of the Sun and the Earth. The equations of motion of the spacecraft in the rotating reference frame accounting solar radiation pressure of the Sun can be expressed as,³²

$\ddot{x}-2\dot{y}=\frac{\partial U}{\partial x},$   (1)

$\ddot{y}+2\dot{x}=\frac{\partial U}{\partial y},$   (2)

$\ddot{z}=\frac{\partial U}{\partial z},$   (3)

Where U is the pseudo-potential of the system and it is expressed as

$U=\frac{\left(x^2+y^2\right)}{2}+\frac{\left(1-\mu \right)q}{r_1}+\frac{\mu }{r_2}.$   (4)

Where q is known as the mass reduction factor, $r_1$  and $r_2$  are the position vectors of the spacecraft from the Sun and the Earth, respectively, and these quantities are defined as,³⁴

$\{\begin{array}{l}q:=1-\frac{F_p}{F_q},\\ r_1:=\sqrt{{\left(x+\mu \right)}^2+y^2+z^2},\\ r_2:=\sqrt{{\left(x+\mu -1\right)}^2+y^2+z^2},\end{array}$   (5)

In Equation (5), $F_p$  and $F_q$  are solar radiation pressure and gravitational attraction forces, respectively. Note that when $q=1$ , the governing equations of motion (1)-(3) reduce to the classical CRTBP.

The Jacobi constant of the motion also exists and is given by

$C=2U\left(x,y,z\right)-\left({\dot{x}}^2+{\dot{y}}^2+{\dot{z}}^2\right).$   (6)

When the kinetic energy is zero, Equation (6) reduces to

$C=2U\left(x,y,z\right),$   (7)

And defines the zero velocity surfaces in the configuration space. These surfaces projected in the rotating xy-plane generate some lines called zero velocity curves.

We construct a third-order analytic approximation using the method of successive approximation (Lindstedt-Poincare method) to compute the halo orbit around the collinear Lagrangian points L₁ and L₂ in the photogravitational Sun-Earth system. The origin is translated at the libration points L₁ and L₂ , and the distance is normalized by taking distance between the Earth to the Lagrangian point as a unit,²⁶ using the new coordinates

$X=\frac{x-\left(1-\mu \mp \gamma \right)}{\gamma },Y=\frac{y}{\gamma },Z=\frac{z}{\gamma },$   (8)

Where X, Y and Z are the new coordinates when origins are shifted at the Lagrangian points, and γ is the distance from the Lagrangian point to the Earth. In Equation (8), the upper (lower) sign corresponds to L₁ (L₂).

Now using the transformation (8), the equations of motion (1)-(3) are expressed as

$\gamma \left(\ddot{X}-2\dot{Y}\right)=\frac{\partial \Psi }{\partial X},$   (9)

$\gamma \left(\ddot{Y}+2\dot{X}\right)=\frac{\partial \Psi }{\partial Y},$   (10)

$\gamma \ddot{Z}=\frac{\partial \Psi }{\partial Z},$   (11)

Where

$\Psi ={\gamma }^2\frac{\left(X^2+Y^2\right)}{2}+\frac{\left(1-\mu \right)q}{R_1}+\frac{\mu }{R_2}.$   (12)

In Equation (12), R₁ and R₂ are given by

$\{\begin{array}{l}{R}_1=\sqrt{{\left(\gamma X+1\pm \gamma \right)}^2+{\left(\gamma Y\right)}^2+{\left(\gamma Z\right)}^2},\\ R_2=\sqrt{{\left(\gamma X\pm \gamma \right)}^2+{\left(\gamma Y\right)}^2+{\left(\gamma Z\right)}^2}.\end{array}$   (13)

The location of the Lagrangian points $L_1$ and $L_2$ from the Earth are computed from the root of the fifth degree polynomial:

${\gamma }^5\pm \left(\mu -3\right){\gamma }^4+\left(3-2\mu \right){\gamma }^3\mp \left[\left(1-q\right)\left(1-\mu \right)\pm \mu \right]{\gamma }^2\pm 2\mu \gamma -\mu =0.$   (14)

In Equation (14) upper and lower signs correspond to $L_1$  and $L_2$  points, respectively. We expand the nonlinear terms, $\frac{\left(1-\mu \right)q}{R_1}+\frac{\mu }{R_2}$  , in Equation (12) using the formula,²⁶

$\frac{1}{\sqrt{{\left(x-A\right)}^2+{\left(y-B\right)}^2+{\left(z-C\right)}^2}}=\frac{1}{D}{{\displaystyle \sum _{m=0}^{\infty }\left(\frac{\rho }{D}\right)}}^m{P}_m\left(\frac{Ax+By+Cz}{D\rho }\right),$   (15)

Where $D^2=A^2+B^2+C^2$ and ${\rho }^2=x^2+y^2+z^2$ , and $P_m\left(\frac{x}{\rho }\right)$ is the $m^{th}$  degree Legendre polynomial of first kind with argument $\frac{x}{\rho }$ . After some algebraic manipulation, the equations of motion (9)-(11) can be written as:^35,36

$\ddot{X}-2\ddot{Y}-\left(1+2c_2\right)X=\frac{\partial }{\partial X}{\displaystyle \sum _{m\ge 3}^{\infty }{c}_m{\rho }^m{P}_m\left(\frac{X}{\rho }\right)},$   (16)

$\ddot{Y}+2\dot{X}+\left(c_2-1\right)Y=\frac{\partial }{\partial Y}{\displaystyle \sum _{m\ge 3}^{\infty }{c}_m{\rho }^m{P}_m\left(\frac{X}{\rho }\right)},$   (17)

$\ddot{Z}+c_{2}Z=\frac{\partial }{\partial Z}{\displaystyle \sum _{m\ge 3}^{\infty }{c}_m{\rho }^m{P}_m\left(\frac{X}{\rho }\right)},$   (18)

Where the left hand side contains the linear terms and the right hand side contains the nonlinear terms. The coefficient $c_m\left(\mu \right)$  is expressed as

$c_m\left(\mu \right)=\frac{1}{{\gamma }^3}\left\{{\left(\pm 1\right)}^m\mu +{\left(-1\right)}^m\frac{q\left(1-\mu \right){\gamma }^{m+1}}{{\left(1\mp \gamma \right)}^{m+1}}\right\},m=0,1,2,\mathrm{...}$   (19)

Where the upper sign is for $L_1$  and the lower one for $L_2$ .

A third-order approximation of Equations (16)-(18) is given by,²⁵

$\ddot{X}-2\dot{Y}-\left(1+2c_2\right)X=\frac{3}{2}{c}_3\left(2X^2-Y^2-Z^2\right)+2c_{4}X\left(2X^2-3Y^2-3Z^2\right),$   (20)

$\ddot{Y}+2\dot{X}+\left(c_2-1\right)Y=-3c_{3}XY-\frac{3}{2}{c}_{4}Y\left(4X^2-Y^2-Z^2\right),$   (21)

$\ddot{Z}+c_{2}Z=-3c_{3}XZ-\frac{3}{2}{c}_{4}Z\left(4X^2-Y^2-Z^2\right).$   (22)

A correction term $\Delta ={\lambda }^2-c_2$ is required for computing the halo orbit which is introduced on the left-hand-side of Equation (22) to make the out-of-plane frequency equals to the in-plane frequency. The new third-order z-equation then becomes:

$\ddot{Z}+{\lambda }^{2}Z=-3c_{3}XZ-\frac{3}{2}{c}_{4}Z\left(4X^2-Y^2-Z^2\right)+\Delta Z.$   (23)

While using the successive approximation procedure, some secular terms arise. To avoid the secular terms, one uses a new independent variable and introduces a frequency connection $\omega$  through $\tau =\omega t$ . The equations of motion (20), (21) & (23) can be then rewritten in terms of new independent variable $\tau$ :

${\omega }^2{X}^{″}-2\omega Y^{\prime }-\left(1+2c_2\right)X=\frac{3}{2}{c}_3\left(2X^2-Y^2-Z^2\right)+2c_{4}X\left(2X^2-3Y^2-3Z^2\right),$   (24)

${\omega }^2{Y}^{″}+2\omega X^{\prime }+\left(c_2-1\right)Y=-3c_{3}XY-\frac{3}{2}{c}_{4}Y\left(4X^2-Y^2-Z^2\right),$   (25)

${\omega }^2{Z}^{″}+{\lambda }^{2}Z=-3c_{3}XZ-\frac{3}{2}{c}_{4}Z\left(4X^2-Y^2-Z^2\right)+\Delta Z,$   (26)

Where $X^{\prime }=\frac{dX}{d\tau },X^{″}=\frac{d^{"2}X}{d{\tau }^2}$  etc.

We assume the solutions of Equations (24)-(26), using the perturbation technique, of the form:

$\{\begin{array}{l}X\left(\tau \right)=\epsilon X_1\left(\tau \right)+{\epsilon }^2{X}_2\left(\tau \right)+{\epsilon }^3{X}_3\left(\tau \right)+\mathrm{...},\\ Y\left(\tau \right)=\epsilon Y_1\left(\tau \right)+{\epsilon }^2{Y}_2\left(\tau \right)+{\epsilon }^3{Y}_3\left(\tau \right)+\mathrm{...},\\ Z\left(\tau \right)=\epsilon Z_1\left(\tau \right)+{\epsilon }^2{Z}_2\left(\tau \right)+{\epsilon }^3{Z}_3\left(\tau \right)+\mathrm{...},\end{array}$   (27)

And

$\omega =1+\epsilon {\omega }_1\left(\tau \right)+{\epsilon }^2{\omega }_2\left(\tau \right)+{\epsilon }^3{\omega }_3\left(\tau \right)+\mathrm{...}+{\epsilon }^i{\omega }_i\left(\tau \right)+\mathrm{...},where{\omega }_i<1.$   (28)

In Equation (28) ∈ is the perturbation parameter. Using Equations (27) & (28) into Equations (24)-(26) and equating the coefficient of the same order of ∈, ∈², and ∈³ from both sides we get the first, second, and third-order equations, respectively.¹

1. First order equations

The first order linearized equations are given by

$\ddot{X}-2\ddot{Y}-\left(1+2c_2\right)X=0,$   (29)

$\ddot{Y}+2\dot{X}+\left(c_2-1\right)Y=0,$   (30)

$\ddot{Z}+{\lambda }^{2}Z=0,$   (31)

Whose periodic solutions are given by:

$\{\begin{array}{l}X\left(t\right)=-A_X\cos \left(\lambda t+\varphi \right),\\ Y\left(t\right)=\kappa A_X\sin \left(\lambda t+\varphi \right),\\ Z\left(t\right)=A_Z\sin \left(\lambda t+\psi \right).\end{array}$   (32)

2. Second order equations

Collecting the coefficients of ${\in }^2$ , we get

${X^{″}}_2-2{Y^{\prime }}_2-\left(1+2c_2\right)X_2=-2{\omega }_1\left({X^{″}}_1-{Y^{\prime }}_1\right)+\frac{3}{2}{c}_3\left(2X_1^2-Y_1^2-Z_1^2\right),$   (33)

${Y^{″}}_2+2{X^{\prime }}_2+\left(c_2-1\right)Y_2=-2{\omega }_1\left({Y^{″}}_1+{X^{\prime }}_1\right)-3c_3{X}_1{Y}_1,$   (34)

${Z^{″}}_2+{\lambda }^2{Z}_2=-2{\omega }_1{Z^{″}}_1-3c_3{X}_1{Z}_1.$   (35)

Now using Equation (32) into (33)-(35), the following equations are obtained

${X^{″}}_2-2n{Y^{\prime }}_2-\left(1+2c_2\right)X_2=2{\omega }_1\lambda A_X\left(\kappa -\lambda \right)\cos {\tau }_1+{\alpha }_1+{\gamma }_1\cos 2{\tau }_1+{\gamma }_2\cos 2{\tau }_2,$   (36)

${Y^{″}}_2+2{X^{\prime }}_2+\left(c_2-1\right)Y_2=2{\omega }_1{A}_X\lambda \left(\kappa \lambda -1\right)\sin {\tau }_1+{\beta }_1\sin 2{\tau }_1,$   (37)

${Z^{″}}_2+{\lambda }^2{Z}_2=2{\omega }_1{A}_Z{\lambda }^2\sin {\tau }_2+{\delta }_1\sin \left({\tau }_1+{\tau }_2\right)+{\delta }_2\sin \left({\tau }_2-{\tau }_1\right),$   (38)

Where ${\tau }_1=\lambda \tau +\varphi ,{\tau }_2=\lambda \tau +\psi$ . Equations (36)-(38) are a set of non-homogenous linear differential equations whose bounded homogenous solution is incorporated from the first-order equations. We need to find only particular solutions of (36)-(38). The secular terms $\sin {\tau }_1$ ,$\cos {\tau }_1$ and $\sin {\tau }_2$ are eliminated by setting ${\omega }_1=0$ . Hence, the solutions of the second-order equations are given by

$\{\begin{array}{l}{X}_2\left(\tau \right)={\rho }_{20}+{\rho }_{21}\cos \left(2{\tau }_1\right)+{\rho }_{22}\cos \left(2{\tau }_2\right),\\ Y_2\left(\tau \right)={\sigma }_{21}\sin \left(2{\tau }_1\right)+{\sigma }_{22}\sin \left(2{\tau }_2\right),\\ Z_2\left(\tau \right)={\kappa }_{21}\sin \left({\tau }_1+{\tau }_2\right)+{\kappa }_{22}\sin \left({\tau }_2-{\tau }_1\right).\end{array}$   (39)

3. Third order equations

Now collecting the coefficients of ${\in }^3$ and setting  ${\omega }_1=0$ , we get

${X^{″}}_3-2{Y^{\prime }}_3-\left(1+2c_2\right)X_3=-2{\omega }_2\left({X^{″}}_1-{Y^{\prime }}_1\right)+3c_3\left(2X_1{X}_2-Y_1{Y}_2-Z_1{Z}_2\right)+2c_{4}X\left(2X_1^2-3Y_1^2-3Z_1^2\right),$   (40)

${Y^{″}}_3+2{X^{\prime }}_3+\left(c_2-1\right)Y_3=-2{\omega }_2\left({Y^{″}}_1+{X^{\prime }}_1\right)-3c_3\left(X_1{Y}_2+X_2{Y}_1\right)-\frac{3}{2}{c}_4{Y}_1\left(4X_1^2-Y_1^2-Z_1^2\right),$   (41)

${Z^{″}}_3+{\lambda }^2{Z}_3=-2{\omega }_2{Z^{″}}_1+\frac{\Delta }{{\in }^2}{Z}_1-3c_3\left(X_2{Z}_1+X_1{Z}_2\right)-\frac{3}{2}{c}_4{Z}_1\left(4X_1^2-Y_1^2-Z_1^2\right).$   (42)

Using Equations (32) and (39) into Equations (40)-(42), we get

${X^{″}}_3-2{Y^{\prime }}_3-\left(1+2c_2\right)X_3=\left[{\upsilon }_1+2{\omega }_2{A}_X\lambda \left(n\kappa -\lambda \right)\right]\cos {\tau }_1+{\gamma }_3\cos {\tau }_1+{\gamma }_4\cos \left(2{\tau }_2+{\tau }_1\right)+{\gamma }_5\cos \left(2{\tau }_2-{\tau }_1\right),$   (43)

${Y^{″}}_3+2{X^{\prime }}_3+\left(c_2-1\right)Y_3=\left[{\nu }_2+2{\omega }_2\lambda A_X\left(\lambda \kappa -1\right)\right]\sin {\tau }_1+{\beta }_3\sin 3{\tau }_1+{\beta }_4\sin \left({\tau }_1+2{\tau }_2\right)+{\beta }_5\sin \left(2{\tau }_2-{\tau }_1\right),$   (44)

${Z^{″}}_3+{\lambda }^2{Z}_3=\left[{\nu }_3+A_Z\left(2{\omega }_2{\lambda }^2+\frac{\Delta }{{\in }^2}\right)\right]\sin {\tau }_2+{\delta }_3\sin 3{\tau }_2+{\delta }_4\sin \left(2{\tau }_1+{\tau }_2\right)+{\delta }_5\sin \left(2{\tau }_1-{\tau }_2\right).$   (45)

The secular terms in the $X_3-Y_3$ equations (43)-(44) and in the $Z_3$  equation (45) cannot be removed by setting a value of ${\omega }_2$ .¹ These terms from Equations. (43)-(45) are removed by adjusting phases of ${\tau }_1$ and ${\tau }_2$  so that $\sin \left(2{\tau }_1-{\tau }_2\right)~\sin {\tau }_2$ which can be achieved by setting the phase constraint relationship

$\varphi =\psi +p\frac{\pi }{2}\text{, }wherep=0,1,2,3.$   (46)

After removing the secular terms from Equation (46), the $Z_3$  solution is bounded when

${\nu }_3+A_Z\left(2{\omega }_2{\lambda }^2+\frac{\Delta }{{\in }^2}\right)+\zeta {\delta }_5=0,\zeta ={\left(-1\right)}^p.$   (47)

The phase constraint (47) reflects the $X_3-Y_3$ equations, each now contains one secular term. The secular terms from both equations are removed by using a single condition from their particular solutions:

$\left({\nu }_1+2{\omega }_2\lambda A_X\left(\kappa -\lambda \right)+\zeta {\gamma }_5\right)-\kappa \left({\nu }_2+2{\omega }_2\lambda A_X\left(\kappa -\lambda n\right)+\zeta {\beta }_5\right)=0.$   (48)

Condition (48) is satisfied if

${\omega }_2=\frac{{\nu }_1-\kappa {\nu }_2+\zeta \left({\gamma }_5-\kappa {\beta }_5\right)}{2\lambda A_X\left[\lambda \left(1+{\kappa }^2\right)-2\kappa \right]}=s_1{A}_X^2+s_2{A}_Z^2,$   (49)

Where similar type of expressions for $s_1$  and $s_2$  can be referred in.¹ Substituting the value of ${\omega }_2$  from Equation (49) into Equation (48), we get

$l_1{A}_X^2+l_2{A}_Z^2+\frac{\Delta }{{\in }^2}=0,$   (50)

Where similar type of expressions for $l_1$  and $l_2$  can be followed from.¹ Equation (50) gives a relationship between the in-plane and the out-of-plane amplitudes. Assuming these constraints, the third-order equations become

${X^{″}}_3-2{Y^{\prime }}_3-\left(1+2c_2\right)X_3=\kappa {\beta }_6\cos {\tau }_1+\left({\gamma }_3+\zeta {\gamma }_4\right)\cos 3{\tau }_1,$   (51)

${Y^{″}}_3+2{X^{\prime }}_3+\left(c_2-1\right)Y_3={\beta }_6\sin {\tau }_1+\left({\beta }_3+\zeta {\beta }_4\right)\sin 3{\tau }_1,$   (52)

${Z^{″}}_3+{\lambda }^2{Z}_3=\{\begin{array}{l}{\left(-1\right)}^{\frac{p}{2}}\left({\delta }_3+{\delta }_4\right)\sin 3{\tau }_1,p=0,2,\\ {\left(-1\right)}^{\frac{p-1}{2}}\left({\delta }_4-{\delta }_3\right)\cos 3{\tau }_1,p=1,3,\end{array}$   (53)

Where ${\beta }_6={\nu }_2+2{\omega }_2\lambda A_X\left(\lambda \kappa -1\right)+\zeta {\beta }_5$ . Thus, the solutions of Equations (51)-(53) are given as

$\begin{array}{l}{X}_3\left(\tau \right)={\rho }_{31}\cos 3{\tau }_1,\\ Y_3\left(\tau \right)={\sigma }_{31}\sin 3{\tau }_1+{\sigma }_{32}\sin {\tau }_1,\\ Z_3\left(\tau \right)=\{\begin{array}{l}{\left(-1\right)}^{\frac{p}{2}}{\kappa }_{31}\sin 3{\tau }_1,p=0,2,\\ {\left(-1\right)}^{\frac{p-1}{2}}{\kappa }_{32}\cos 3{\tau }_1,p=1,3.\end{array}\end{array}$   (54)

The expressions for all the coefficients can be referred to.²⁹

4. Final approximation

Halo orbits of third-order approximations are obtained on removing $\in$ from all solutions of equations by using the transformation $A_X\mapsto A_X/\in$ and $A_Z\mapsto A_Z/\in$ . Then one can use $A_X$  or $A_Z$  as a small parameter. Combining the above computed solutions, the third-order approximate solution is thus given by

$\begin{array}{l}X\left(\tau \right)={\rho }_{20}-A_X\cos {\tau }_1+\left({\rho }_{21}+\zeta {\rho }_{22}\right)\cos 2{\tau }_1+{\rho }_{31}\cos 3{\tau }_1,\\ Y\left(\tau \right)=\left(\kappa A_X+{\sigma }_{32}\right)\sin {\tau }_1+\left({\sigma }_{21}+\zeta {\sigma }_{22}\right)\sin 2{\tau }_1+{\sigma }_{31}\sin 3{\tau }_1,\\ Z\left(\tau \right)=\{\begin{array}{l}{\left(-1\right)}^{\frac{p}{2}}\left[A_Z\sin {\tau }_1+{\kappa }_{21}\sin 2{\tau }_1+{\kappa }_{31}\sin 3{\tau }_1\right],p=0,2\\ {\left(-1\right)}^{\frac{p-1}{2}}\left[A_Z\cos {\tau }_1+{\kappa }_{21}\cos 2{\tau }_1+{\kappa }_{22}+{\kappa }_{32}\cos 3{\tau }_1\right],p=1,3.\end{array}\end{array}$   (55)

Time period (in non-dimensional form) of the halo orbit is expressed as

$T_{halo}=\frac{2\pi }{\lambda \omega },where\omega =1+{\omega }_1+{\omega }_2;{\omega }_1=0.$   (56)

In this section, Newton’s method of differential correction is briefly described for the numerical computation of halo orbit. Assume X denote a column vector containing all of the six state variables of the governing equations of motion, i.e.,

$X={\left[xyz\dot{x}\dot{y}\dot{z}\right]}^T,$   (57)

Where superscript “T” denotes the transpose.

The state transition matrix (STM), $\Phi \left(t,t_0\right)$ , is a $6\times 6$ matrix composed of the partial derivatives of the state:

$\Phi \left(t,t_0\right)=\frac{\partial X\left(t\right)}{\partial X\left(t_0\right)},$   (58)

With initial conditions $\Phi \left(t_0,t_0\right)=I$ . Note that the state transition matrix is called monodromy matrix for the full periodic orbit. The eigenvalues of the monodromy matrix tells about the stability of the halo orbit.

The STM is propagated using the relationship:

$\frac{d\Phi \left(t,t_0\right)}{dt}=A\left(t\right)\Phi \left(t,t_0\right),$   (59)

Where the matrix $A\left(t\right)$ is known as variational matrix and is made of the partial derivatives of the state derivative with respect to the state variables, i.e.,

$A\left(t\right)=\frac{\partial \dot{X}\left(t\right)}{\partial X\left(t\right)}.$   (60)

The variation matrix $A\left(t\right)$ can be partitioned into four sub-matrices:

$A\left(t\right)=\left(\begin{array}{cc}O& I\\ \Upsilon & 2\Omega \end{array}\right),$   (61)

Where

$\begin{array}{l}O=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right],I=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right],\Omega =\left[\begin{array}{ccc}0& 1& 0\\ -1& 0& 0\\ 0& 0& 0\end{array}\right]\\ \Upsilon =\left[\begin{array}{ccc}\frac{\partial \ddot{x}}{\partial x}& \frac{\partial \ddot{x}}{\partial y}& \frac{\partial \ddot{x}}{\partial z}\\ \frac{\partial \ddot{y}}{\partial x}& \frac{\partial \ddot{y}}{\partial y}& \frac{\partial \ddot{y}}{\partial z}\\ \frac{\partial \ddot{z}}{\partial x}& \frac{\partial \ddot{z}}{\partial y}& \frac{\partial \ddot{z}}{\partial z}\end{array}\right]=\left[\begin{array}{ccc}{U}_{XX}& U_{XY}& U_{XZ}\\ U_{YX}& U_{YY}& U_{YZ}\\ U_{ZX}& U_{ZY}& U_{ZZ}\end{array}\right]\end{array}$   (62)

Note that the matrix $\Upsilon$  is a symmetric matrix of second order partial derivatives of U  with respect to x, y, and z evaluated along the orbit. Thus, Equation (59) represents a system of 36 first-order differential equations. These equations, coupled with the equations of motion (1)-(3), are the basic equations that define the dynamical model in the photo-gravitational CRTBP accounting solar radiation pressure. Trajectories are computed by simultaneous numerical integration of the 42 first-order differential equations. It can be easily seen that the governing equations of motion (1)-(3) are symmetric about the xz-plane by using the transformation $y\to -y$  and $t\to -t$ .

Let $X\left(t_0\right)$ be the state of a periodic symmetric orbit at the xz-plane crossing and let $X\left(t_{\raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)$ denotes the state of the orbit half of its orbital period later at the xz-plane. If the orbit is periodic and symmetric about the xz- plane, then

$X\left(t_0\right)={\left[x_{0}0z_{0}0\dot{y}0\right]}^T\Rightarrow X\left(t_{\raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)={\left[x_{\raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}0z_{\raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}0{\dot{y}}_{\raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}0\right]}^T.$   (63)

Assume that $\widehat{X}\left(t_0\right)$ be an initial state of a desirable state. Integrating this state forward in time until the next xz-plane crossing, we obtain the state $\widehat{X}\left(t_{\raisebox{1ex}{$\widehat{T}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)$ :

$\widehat{X}\left(t_{\raisebox{1ex}{$\widehat{T}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)={\left[x0z_{\raisebox{1ex}{$\widehat{T}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\dot{x}}_{\raisebox{1ex}{$\widehat{T}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\dot{y}}_{\raisebox{1ex}{$\widehat{T}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{\dot{z}}_{\raisebox{1ex}{$\widehat{T}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right]}^T.$   (64)

We adjust the initial state of the trajectory in such a way so that the values of ${\dot{x}}_{\raisebox{1ex}{$\widehat{T}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$ and ${\dot{z}}_{\raisebox{1ex}{$\widehat{T}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}$  become zero. Note that by adjusting the initial state, not only the values of and change, but the propagation time, $\raisebox{1ex}{$\widehat{T}$}\!\left/ \!\raisebox{-1ex}{$2$}\right.$ , needed to penetrate the xz -plane also changes. In order to target a proper state $X\left(t_{\raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)$ , one may vary the initial values of x, z and/or $\dot{y}$ . The linearized system of equations relating the final state to the initial state can be written as:

$\delta X\left(t_{\raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)\approx \Phi \left(t_{\raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$2$}\right.},t_0\right)\delta X\left(t_0\right)+\frac{\partial X}{\partial t}\delta \left(\raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right),$   (65)

Where $\delta X\left(t_{\raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}\right)$ denotes the deviation in the final state due to a deviation in the initial state, $\delta X\left(t_0\right)$ , and a corresponding deviation in the orbit’s period, $\delta \left(\raisebox{1ex}{$T$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)$ .

The variation in the locations of $L_1$  and $L_2$  with the mass reduction factor, q are given in Table 1 from the Barycenter. It can be observed that as the value of q decreases, the distance between $L_1(L_2)$   and the Barycenter decreases (increases). Thus, as solar radiation pressure dominates, the location of $L_1$ moves towards the Sun while that of $L_2$ moves away from the Sun.

Halo orbits are computed using the constructed third-order analytic approximate solution as the starting initial guess. Figures 2–5 are generated using the following characteristic properties of ISEE-3 mission:³⁷ mass of the spacecraft = 435kg; solar reflectivity constant, $k=1.2561$ ; spacecraft effective cross sectional area, $A=3.55$ m²; speed of light, $c=2.998\times {10}^8$ m/sec; solar light flux, $S_0=\text{1352}\text{.098 }$  kg/sec² at one astronomical unit from the Sun. We have chosen the out-of-plane amplitude, $A_Z=$ 1,10000km of ISEE-3, for the sake of simplicity, the corresponding value of the in-plane amplitude, $A_X$  is 2,06000km.

Figures 2 & 3 depict the projections of xy, xz, and yz- planes of northern branch of the halo orbit around $L_1$  and $L_2$ , respectively, whereas Figure 4 depicts its three dimensional (3D) state. Similarly, its southern branch can be obtained by changing the sign of z since both branches are mirror images to each other. Jacobi constant of the halo orbit around $L_1$ is $C_{halo}=$ 3.00082686598735 while it is $C_{halo}=$ 3.00082167380548 for $L_2$ . The halo orbit and it’s zero velocity curves around $L_1$  and $L_2$  are shown in Figure 5. It can be observed that the halo orbit lies in the neck and goes around  $L_1(L_2)$ .

Figure 2 Projection of (A) xy-plane, (B) xz-plane, and (C) yz-plane of the halo orbit around L₁.

Figure 3 Projection of (A) xy-plane, (B) xz-plane, and (C) yz-plane of the halo orbit around L₂.

Figure 4 3D state of the halo orbit around L₁ (left) and L₂ (right).

Figure 5 Halo orbit and its zero velocity surfaces around L₁ (left) and L₂ (right).

The effect of solar radiation pressure on the spacecraft’s velocity in a halo orbit around $L_1$  and $L_2$  with q is shown in Figure 6. It is observed that as solar radiation pressure increases, velocity of the spacecraft increases in a halo orbit around $L_1$  while it decreases around $L_2$ . Also maximum (minimum) magnitude of the spacecraft’s velocity are 2.918$\times {10}^{-1}$ (9.823$\times {10}^{-2}$ ), 2.947$\times {10}^{-1}$ (9.985$\times {10}^{-2}$ ), and 3.043$\times {10}^{-1}$  (1.056$\times {10}^{-1}$ ) km/sec for q = 1, 0.998668, and 0.984334, respectively around $L_1$ . For q=1, 0.999668, and 0.999334, maximum (minimum) magnitude of the spacecraft’s velocity are 2.963$\times {10}^{-1}$ (9.942$\times {10}^{-2}$ ), 2.956$\times {10}^{-1}$ (9.903$\times {10}^{-2}$ ), and 2.948$\times {10}^{-1}$ (9.864$\times {10}^{-2}$ )km/sec, respectively around $L_2$ . Maximum (minimum) values of the spacecraft’s velocity are obtained at the xz (xy)-plane crossing time. Time period of the halo orbit and Jacobi constant with q are shown in Tables 2 & 3 for $L_1$  and $L_2$ , respectively. As solar radiation pressure prevails, time period of the halo orbit increases whereas Jacobi constant decreases about both libration points. The effect of solar radiation pressure on shape of the halo orbit around $L_1$  and $L_2$ are shown in Figures 7 & 8, respectively. One more important observation is that as solar radiation pressure dominates, shape of the halo orbit increases and move towards the Sun for $L_1$  while it shrinks towards the Earth around $L_2$ .

Figure 6 Velocity variation of the spacecraft in the halo orbit with q around L₁ (left) and L₂ (right).

Figure 7 Variation in shape of (A) xy-plane, (B) xz-plane, and (C) yz-plane of the halo orbit with q around L₁.

Figure 8 Variation in shape of (A) xy-plane, (B) xz-plane, and (C) yz-plane of the halo orbit with q around L₂.

According to Floquet theory,^26,38,39 stability of the halo orbits is described by the eigenvalues of its monodromy matrix. The monodromy matrix corresponding to the halo orbit around $L_1$  and $L_2$ has six eigenvalues $\left({\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4,{\lambda }_5,{\lambda }_6\right)$ which are given by

$\left(1732.916,0.0005770619,1.0000018,1.0000018,0.9999982,0.9968152\pm 0.0797459i\right),$   (71)

And

$\left(1664.2099,0.00060088575,0.9970228\pm 0.0771079i,1.000000\pm 0.0000024i\right),$   (72)

for $L_1$  and $L_2$ , respectively.

The periodic orbit is stable only if the modulus of all eigenvalues of its monodromy matrix is less than one.²⁸ It can be noted that for these orbits two eigenvalues are non-unity (${\lambda }_1$  and ${\lambda }_2$ ) near both $L_1$ and $L_2$ , and the complex eigenvalues lie on the unit circle. Thus, there exists both stable and unstable halo orbits, and orbits near the halo which remain near the halo for all time around $L_1$  and $L_2$ . The eigenvectors corresponding to stable and unstable eigenvalues direct stable and unstable manifolds of the orbit whereas complex eigenvalues correspond to the center directions of the orbit. The eigenvalues show that halo orbits have and typecharacteristics behavior around $L_1$ and $L_2$ points, respectively.

In this study, a third-order analytic approximate solution using the Lindstedt-Poincare method and Newton’s single step differential correction scheme are used to compute the halo orbit analytically and numerically around the collinear points $L_1$  and $L_2$  in the photogravitational circular restricted three-body problem accounting radiation pressure of the Sun. The effects of solar radiation pressure are studied around both collinear libration points. For q = 1, 0.999934, 0.999668, and 0.999336, the Barycentric distance of $L_1(L_2)$ are 1.481019$\times {10}^8$ (1.511071$\times {10}^8$ ), 1.481008$\times {10}^8$ (1.511082$\times {10}^8$ ), 1.480963$\times {10}^8$  (1.511126$\times {10}^8$ ), and 1.480906$\times {10}^8$  (1.511183$\times {10}^8$ ) kilometers, respectively for $L_1(L_2)$ which shows that as solar radiation pressure dominates, the distance between $L_1(L_2)$ and the Barycenter decreases (increases). The achieved maximum (minimum) velocity (in terms of magnitude) of the spacecraft are 2.918$\times {10}^{-1}$ (9.823$\times {10}^{-2}$ ), 2.947$\times {10}^{-1}$ (9.985$\times {10}^{-2}$ ), and 3.043$\times {10}^{-1}$  (1.056$\times {10}^{-1}$ )km/sec around $L_1$  whereas around $L_2$ are 2.963$\times {10}^{-1}$ (9.942$\times {10}^{-2}$ ), 2.956$\times {10}^{-1}$ (9.903$\times {10}^{-2}$ ), and 2.948$\times {10}^{-1}$ (9.864$\times {10}^{-2}$ )km/sec, respectively for q=1, 0.999668, and 0.999334, respectively. In other words, as solar radiation pressure increases, velocity of the spacecraft increases around $L_1$  point while it decreases around $L_2$ point. It is found that time period of the halo orbit increases around both $L_1$  and $L_2$ points. Further, as solar radiation pressure dominates, shape of the halo orbit around $L_1$ increases and moves towards the Sun while it shrinks around $L_1$  and moves towards the Earth. The eigenvalues of the monodromy matrix depict that the halo orbits have $saddle\times center$ ($L_1$ ) and $saddle\times center\times center$  ($L_2$ ) type of behavior in the photogravitational circular restricted three-body problem for the Sun-Earth system.

None.

Authors declare there is no conflict of interest.

1. Thurman R, Worfolk PA. The geometry of halo orbits in the circular restricted three-body problem. Technical report. CiteSeerX, USA; 1996.
2. Schuerman DW. The restricted three-body problem including radiation pressure. The Astrophysical Journal. 1980;238:337–342.
3. Simmons JFL, McDonald AJC, Brown JC. The restricted 3-body problem with radiation pressure. Celestial mechanics. 1985;35(2):145–187.
4. Ragos O, Zafiropoulos FA. A numerical study of the influence of the Poynting-Robertson effect on the equilibrium points of the photogravitational restricted three-body problem. I. Coplanar case. Astronomy and Astrophysics. 1995;300:579–590.
5. McInnes CR. Solar Sailing: Technology, Dynamics and Mission Applications. Springer-Verlag, Berlin, Germany. 1990.
6. Waters TJ, McInnes CR. Periodic orbits above the ecliptic in the solar-sail restricted three-body problem. Journal of Guidance Control and Dynamics. 2007;30(3):687–693.
7. Kishor R, Kushvah BS. Periodic orbits in the generalized photogravitational Chermnykh-like problem with power-law profile. Astrophysics and Space Science. 2013;344(2):333–346.
8. Verrier P, Waters W, Sieber J. Evolution of the L1 halo family in the radial solar sail circular restricted three-body problem. Celestial Mechanics and Dynamical Astronomy. 2014;120(4):373–400.
9. Srivastava VK, Kumar J, Kushvah BS. The effects of oblateness and solar radiation pressure on halo orbits in the photogravitational Sun-Earth system. Acta Astronautica. 2014;129:389–399.
10. Moulton FR. An Introduction to Celestial Mechanics. Kessinger Publishing, USA. 1914.
11. McCuskey SW. Introduction to Celestial Mechanics. Dover Publications, USA. 1963.
12. Szebehely V. Theory of Orbits, The Restricted Problem of Three Bodies. Academic Press, USA. 1967.
13. Baoyin H, McInnes C. Solar sail halo orbits at the Sun-Earth artificial L1 point. Celestial Mechanics and Dynamical Astronomy. 2006;94(2):155–171.
14. Heiligers J, McInnes C. Novel Solar Sail Mission Concepts for Space Weather Forecasting. 24th AAS/AIAA Space Flight Mechanics Meeting. Santa Fe, New Mexico, USA, 2014. p. 1–20.
15. Farquhar RW, Kamel AA. Quasi-periodic orbits about the translunar libration point. Celestial mechanics. 1973;7(4):458–473.
16. D’Amario LA, Bright LE, Wolf AA. Galileo trajectory design. Space Science Reviews. 1992;60(1–4):23–78.
17. Martin WL. Libration point trajectory design. Numerical Algorithms. 1997;14(1–3):153–164.
18. Gomez G, Koon WS, Martin WL, et al. Connecting orbits and invariant manifolds in the spatial restricted three-body problem. Nonlinearity. 2004;17(5):1571–1606.
19. Xu R, Cui P, Qiao D, et al. Design and optimization of trajectory to Near-Earth asteroid for sample return mission using gravity assists. Advances in Space Research. 2007;40(20):220–225.
20. Pergola P, Geurts K, Casaregola C, et al. Earth–Mars halo to halo low thrust manifold transfers. Celestial Mechanics and Dynamical Astronomy. 2009;105.
21. Landgraf M, Renk F, Vogeleer BD. Mission design and analysis of European astrophysics missions orbiting libration points. Acta Astronautica. 2013;84:49–55.
22. Vaquero M, Howell KC. Design of transfer trajectories between resonant orbits in the Earth–Moon restricted problem. Acta Astronautica. 2014;94(1):302–317.
23. Farquhar RW. The control and Use of Libration-Point Satellites. Department of Aeronautics and Astronautics, Stanford University, USA; 1968. p. 1–138.
24. Zazzera FB, Topputo F, Massari M. Assessment of Mission Design Including Utilization of Libration Points and Weak Stability Boundaries. ESA, ACT, France; 2005. p. 1–173.
25. Richardson DL. Analytic construction of periodic orbits about the collinear points. Celestial Mechanics. 1980;22(2):241–253.
26. Koon WS, Lo MW, Marsden JE. Dynamical Systems: The Three-Body Problem and Space Mission Design. Interdisciplinary Applied Mathematics, Springer, Berlin, Germany; 2011. p. 1–331.
27. Breakwell JV, Brown JV. The ‘Halo’ family of 3-dimensional periodic orbits in the Earth-Moon restricted 3-body problem. Celestial mechanics. 1979;20(4):389–404.
28. Howell KC. Three-dimensional, periodic, ‘halo’ orbits. Celestial mechanics. 1984;32(1):53–71.
29. Tiwary RD, Kushvah BS. Computation of halo orbits in the photogravitational Sun-Earth system with oblateness. Astrophysics and Space Science. 2015;73:357.
30. Dunham DW, Farquhar RW. Libration point missions, 1978–2002. In: Gómez G, et al. editors. Proceedings of the Conference on Libration Point Orbits and Applications, World Scientific, Singapore; 2003.
31. Radzievskii VV. Astronomical Journal USSR. 1950;27(5):250.
32. Eapen RT, Sharma RK. A study of halo orbits at the Sun–Mars L1 Lagrangian point in the photogravitational restricted three-body problem. Astrophysics and Space Science. 2014;352(2):437–441.
33. Bell JL. The impact of solar radiation pressure on Sun-Earth L1 Libration point orbits. MS thesis, Department of Aeronautics and Astronautics, Purdue University, USA. 1991.
34. Srivastava VK, Kumar J, Kushvah BS. Regularization of circular restricted three-body problem accounting radiation pressure and oblateness. Astrophysics and Space Science. 2017;362.
35. Jorba A, Masdemont J. Dynamics in the center manifold of the collinear points of the restricted three body problem. Physica D. 1999;132(1–2):189–213.
36. Bucciarelli S, Ceccaroni M, Celletti A. Annali di matematica pura e applicata 01, Springer, Germany. 2015.
37. Pernicka HJ. The numerical determination of nominal Libration point trajectories and development of a stationkeeping strategy. Department of Aeronautics and Astronautics, Purdue University, USA. 1990.
38. Baig S, McInnes MC. Artificial halo orbits for low-thrust propulsion spacecraft. Celestial Mechanics and Dynamical Astronomy. 2009;104(4):321–335.
39. Gomez G, Mondelo, JM. The dynamics around the collinear equilibrium points of the RTBP. Physica D. 2001;157(4):283–321.