2.1 Equilibrium Point and Stability Analysis

Let x=(x, y, z, x˙, y˙, z˙) ∈ R6 be the state vector of a massless particle and f:R6 → R6 be

The equilibrium points L _i , i = 1,2,3,4,5 are defined as the zeros of f(x). To determine the linear stability of each equilibrium, one needs to translate the origin to the position of this equilibrium point and linearize the equations of motions around this point. Thus, the linearization over any of these equilibrium points is

where DLi is the derivative of f(x) computed at the equilibrium point L _i .

To determine the linear stability of the equilibrium points (L _i , i = 1,2,3,4 and 5), it is necessary to transfer the origin of the coordinate system to the position of the equilibrium points (x0,y0,z0) and then to linearize the equations of motion around these points, obtaining the results shown below.

where the partial derivatives in (x0,y0,z0) mean that the value is calculated at the equilibrium point being analyzed. Partial derivatives are shown in equations 8 - 13.

In equations 5 - 7, ξ, η and ζ represent the position of the particle with respect to the equilibrium point. Through numerical analysis, we observed that the equilibrium points exist only in the xy plane, regardless of the values assigned to d, μ* and k. Due to the fact that the equilibrium points for the rotating mass dipole model are in the xy plane, the equation 7 is decoupled (it does not depend on ξ and η); therefore, the equation of motion 7 becomes

where ϑ is constant and depends on the values assigned to d, μ* and k. Equation 14 shows that the motion perpendicular to the xy plane is periodic with frequency ω = ϑ. Motion in the z direction is therefore limited with

where c₃ and c₄ are integration constants.

When the motion is in the xy plane, the non-trivial characteristic roots of the equations 5, 6 were obtained in ^{Barbosa Torres dos Santos et al.(2017)} (considering k = 1.). The linearization around L₁ and L₂ provides a pair of real eigenvalues (saddle), corresponding to one-dimensional stable and unstable manifolds, and one pair of imaginary eigenvalues, suggesting a two-dimensional central subspace in the plane xy, which accounts for an oscillatory behavior around the equilibrium point of the linear system ^Howell(1982), ^{Haapala et al.(2015)}. Hence, in general, for L₁ and L₂, the stability type is saddle×center×center for the problem studied here and also for the CRTPB considering 0 < μ ≤ 0.5. The Lyapunov Center Theorem guarantees, for the planar case, the existence of a one-parameter family of periodic orbits emanating from each of the collinear equilibrium points. Thus, for the spatial case, two one-parameter families of periodic orbits around L₁ and L₂ are expected. It was observed that the nature of the eigenvalues of the collinear equilibrium points is not altered when we vary d, μ* and k.

Consider the linearized dynamics around the L₁ equilibrium point. We will adopt the coordinates x'=(ξ; η; u; v), where u and v are the velocities in the x and y direction, respectively, for the physical variables in the linearized planar system. To differentiate, we will use the coordinates x0=(x;y;x˙;y˙) for the physical variables in the nonlinear system and, finally, y = (y1; y2; y3; y4) for the variables in the diagonalized system. We know that if we choose an initial condition anywhere near the equilibrium point the real components of the eigenvalues (stable and unstable) will dominate the particle’s behavior. But instead of specifying any initial condition for the system, we want to find an orbit around the equilibrium point L₁, for example, with some desired behavior, such as a periodic orbit. This becomes easy if we use the diagonalized system (y0) to determine the initial conditions. As we want to minimize the component in the unstable direction of the non-linear path, we must choose the initial conditions that correspond to the harmonic motion of the linear system. Thus, we choose the initial condition in the diagonalized system as y0 = (0; 0; y3; y4), where the non-zero initial values can be complex numbers and are intended to amplify the oscillatory terms. Null terms have the function of nullifying exponential (unstable) terms. In fact, if we want to get real solutions at the x' coordinates, we must consider y3 and y4 as complex conjugates. Transforming these conditions back to the original coordinates of the linear system, from the transformation x0 = Ty0, we find the initial conditions in the linearized system x'0 = (x', y', u, v), where T is the matrix of the eigenvectors of the state transition matrix A. The Jacobian matrix A contains the Hessian pseudo potential derived from the truncated Taylor series expansion over the reference solution.

Due to the fact that these initial conditions were chosen to nullify the unstable and stable eigenvectors, they provide a harmonic movement in the linear system.

Now that we have the initial conditions for the linear system, we want to find a periodic planar orbit in the nonlinear system.

We note that the potential function for the system studied here depends only on the distances that a spacecraft are from the primary bodies, that is, it has symmetry with respect to the x-axis. Taking advantage of the fact that the planar orbits are symmetrical with respect to the x -axis, the initial state vector takes the form x0 = x0 0, 0, y˙0 T. These symmetries were used to find symmetric periodic orbits. This was done by determining the initial conditions, on the x-axis, where the initial velocity is perpendicular to this axis (y˙) and then the integration is done until the path returns by crossing the x-axis with the speed orientation y˙f opposite to the initial condition. This orbit can be used as an initial guess to use Newton’s method, where the target state is quoted above; that is, that the orbit returns to x-axis with normal velocity. The equations of motion and the state transition matrix are incorporated numerically until the trajectory crosses the x-axis again. The final desired condition has the following form: xf = xf 0, 0, y˙f T.

4.1 Planar Orbits

Figure 3 shows a family of planar orbits around L₁ with μ*=10-5 and d = 0. The orbits obtained do not intersect the asteroid although, as seen in Figure 3, as the amplitude increases along the family, the orbits expand from the vicinity of the equilibrium point towards the surface of the secondary body (black asterisk).

Figure 3: Planar orbits around of the equilibrium point L₁ considering μ*=10-5 and d = 0. The color figure can be viewed online. 

In Figure 3, the red orbits indicate where bifurcations occur, that is, when one of the stability indices s₁ or s₂ reaches the critical value 2. Note in Figure 3 that the maximum position reached by the infinitesimal mass body in the x component, when the second bifurcation occurs, is greater than the position of the secondary body.

Although many bifurcations exist in dynamical systems, only two types of bifurcation are of particular interest for the focus of this work; the pitchfork and period-multiplying bifurcations.

A family of periodic orbits undergoes a pitchfork bifurcation when the stability of the periodic orbit changes as a parameter evolves, which in our case is the energy constant. During this type of local bifurcation, a pair of eigenvalues (non-trivial) of the monodromy matrix pass through the critical values λ1=1/λ1 (or λ2=1/λ2) = + 1 of the unit circle. Consequently, the stability index passes through s₁ (or s₂) = 2 ^{Bosanac(2016)}. In addition, the stability of the periodic orbits changes along a family, and an additional family of a similar period is formed. This new family of orbits has the same stability as the members of the original family before the bifurcation arose. On the other hand, a period-doubling bifurcation is identified when a pair of non-trivial eigenvalues (λ1,2 and 1/λ1,2, where λ1,2 means λ1 or λ2), passes through λ1,2=1/λ1,2=-1 of the unit circle. Therefore, it represents a critical value of the stability index, such that s_1,2 = -2 ^{Bosanac(2016)}.

When building the families of planar orbits, with d = 0 and μ* = 10^-5, we observe that the stability index (s₂) reaches the critical value three times for planar orbits around L₁ and L₂, as seen in Figures 4. In both figures, the horizontal axis displays the minimum x value along the orbits. The stability index s1 does not reach the critical value for this μ* and d.

For μ*=10-5, the equilibrium point L₁ is located at x = 0.981278, y = 0 and z = 0, while L₂ is at position x =1.01892, y = 0 and z = 0. The orbits with smaller amplitudes are close to the equilibrium point (right side of Figures 4) and the first bifurcation occurs for a small amplitude orbit (x ≈ 0.97583 for L₁ and x ≈ 1.01577 for L₂). As we continue the planar family, the stability index s₂ shown in Figure 4 continues to increase, reaches a maximum, decreases and reaches the value 2 again, where another bifurcation occurs, with x ≈ 0.99408 for L₁ and x ≈ 1.0056 for L₂. As we continue the families of planar orbits around L₁ and L₂, the stability index decreases, reache a minimum, increases and again and reaches the critical value 2, where another bifurcation occurs. After the third bifurcation, the stability index further increases and we did not detect additional bifurcations given that, as the orbits are very close to the center of mass of the secondary body, our Newton method looses track of planar orbits, converging to a completely different family of orbits.

Figure 4: Stability index (s₂) around L₁ (red) and L₂ (green) considering d = 0 and μ*=10-5. The color figure can be viewed online. 

Figures 5 (a), (b) and (c) provide information about the types of the bifurcations that occur along the family of planar orbits. For μ*=10-5 and d = 0, analyzing the path of the characteristic multipliers in Figures 5 (a) and (b), we find that the first bifurcation is a supercritical pitchfork bifurcation, while the second one corresponds to a subcritical pitchfork case. This suggests that new families of periodic orbits appear in those regions when the bifurcation occurs ^{Feng et al.(2016)}. In fact, after the first bifurcation (low amplitude periodic orbit), it is possible to detect halo orbits, while after the second bifurcation the family of axial orbits appears Grebow(2006). Unlike the planar Lyapunov orbits, halo and axial orbits are three-dimensional.

Figure 5: (a) Behavior of the characteristic multipliers at the first pithckfork bifurcation around L₁ and L₂. (b) Behavior of the characteristic multipliers at the second pithckfork bifurcation around L₁ and L₂. (c) Behavior of the characteristic multipliers that leads to the period-doubling bifurcation around L₁ and L₂. In these cases we consider d = 0 and μ* = 10^-5. The color figure can be viewed online. 

Figure 5 (c) shows the behavior of the eigenvalues at the third bifurcation. The characteristic multipliers start in the imaginary plane and move until they collide on the negative real axis and start to reach only real values on the negative axis. Thus, the eigenvalues indicate a period-doubling bifurcation.

Figure 6 and 7 provides information about the stability index (considering the values of s₂), around L₁ and L₂, respectively, when we increase the dipole dimension from 0 meters, that is, the body is modeled as a mass point, up to the dimension of 2000 meters. In this analysis we consider the constant mass ratio in the value of μ* = 10-5. When we consider the dipole as a point mass body (d = 0), it is possible to observe three bifurcations (the red curve passes through the critical value three times). We can observe that as we increase the dimension of the secondary, the second bifurcation points in the planar orbits around L₁ and L₁ cease to exist because the trajectories collide with the secondary body. This is because, as the dimension of the dipole varies and the planar orbits approach the secondary body, our Newton method looses track of the planar orbits, converging to a completely different family of the orbits. Note that, the larger the dipole size, the smaller the planar orbit family found.

Figure 6: Planar orbit stability index around L₁ for different values of d. The color figure can be viewed online. 

Figure 7: Planar orbit stability index around L₂ for different values of d. The color figure can be viewed online. 

4.2 Influence of the Mass Parameter and the Size of the Dipole on the Planar Orbits

Now, we investigate how the planar orbits evolve as a function of the dipole size and mass ratio in canonical units. With the normalization factor being D = 12000 meters, the dipole sizes used in our study were d = 0, 500, 1000, 1500 and 2000 meters.

Figures 8 provide information about the stability index s₁ of the planar orbits around L₁, respectively, as a function of d and μ*. In both figures, the color code accounts for the size of the dipole (d).

Figure 8: Stability index (s₁) of the planar orbits around L₁ for different values of d and μ*. The color figure can be viewed online. 

First, we investigate the solutions as d varies and μ* is kept constant. Note that, in Figures 8, in general, when the size of the dipole increases, the planar orbits become more unstable. This means that the larger the secondary body, the more unstable the planar orbits are.

If we consider d = 0, which corresponds to the CRTBP, we observe that as μ* increases, the orbits become increasingly unstable. On the other hand, when the elongated form of the secondary body is taken into account, s₁ becomes smaller as μ* increases, and it only increases again after μ* = 10^-1. This information is important for space missions, since a high value in the stability index (s _i ) indicates a divergent mode that moves the spacecraft away from the vicinity of the orbit quickly. In general, the stability index is directly related to the space vehicle’s orbital maintenance costs and inversely related to the transfer costs. This same analysis was performed around the L₂ equilibrium point, where we found similar results.

Next, we analyze the period of the planar orbits in terms of d and μ* around L₁. As shown in Figures 9, for low amplitude, as d increases, with μ* kept constant, the period of the planar orbits decreases. This is because the mass distribution of the secondary body allows part of the mass of the asteroid to be closer to the negligible mass particle, causing the gravitational attraction to become larger, thus increasing the acceleration in the vicinity of the secondary body and decreasing the orbital period. On the other hand, when the x-amplitude is large, the results can be inverted, as shown in Figure 9. In general, when the amplitude of the orbit increases, the orbital period becomes longer. Similar results were found in the vicinity of L₂.

Figure 9: Period of planar orbits around L₁ for different values of d and μ*. The color figure can be viewed online. 

Considering the family with d = 0, when μ* increases the period of the orbits remains similar, except when μ* = 10^-1. Conversely, when the elongation of the secondary body is considered, in general, for a given value of d the larger the mass ratio, the longer the orbital period.

Finally, we analyze the energy of the system in terms of d and μ*, as shown in Figure 10.

Figure 10: Jacobi constant of planar orbits around L₁ for different values of d and μ*. The color figure can be viewed online. 

We find that when d or μ* increase, the energy required to orbit a given equilibrium point decreases. That is, the more elongated the secondary body and the larger the value of μ*, the less energy is needed to orbit a given equilibrium point. This also means that, as the size of the dipole increases or as the mass ratio of the system increases, the bifurcations occur at lower energies. The same analysis performed for L₁ can be done for L₂.

4.3 Computing Halo Orbits

Halo orbits are a three-dimensional branch of planar orbits that appear when the planar orbit stability index reaches the critical value s₂ = 2. Figure 11 shows a family of halo orbits around L₁ with μ* = 10^-5 and d = 0. The orbits are in three-dimensional space and as the amplitude increases along the family, the halo orbits expand from the vicinity of the equilibrium point towards the surface of the secondary (black asterisk).

Figure 11: Halo orbits around the L₁ equilibrium point considering μ*=10-5 and d = 0. The color figure can be viewed online. 

To find the initial conditions of the halo orbit, we keep the coordinate x₀ fixed and search for z0*, y˙0* and T/2* such that x˙*(T/2*), z˙*(T/2*) and y*(T/2*) are all null. Then, to find the halo orbit, we use as initial guess the position x₀, velocity (y˙) and period (T) of the planar orbit when the stability index s2=2. Knowing these initial conditions, all that remains is to determine the initial guess of the position on the z axis, such that we can find the halo orbit. Because the halo and planar orbit are similar (when s2=2), the position on the z axis of the halo orbit must have a very small value (almost planar orbit). Thus, in this work, the value of z₀ = 0.0001 canonical units was used as the initial guess for the position on the axis z. A Newton method for this problem is

with x = (z,y˙, T/2) and x0 = (z0,y˙0, T0/2). Here (z0, y˙0, T0/2) is the initial guess of the halo orbit.

The differential is

where ϕi,j are elements of the monodromy matrix, g:U⊂R6→R6 is the vector field of the restricted synchronous three-body problem, z(T/2)=ϕ3(x0,0,z,0,y,0,T/2) and y˙(T/2) = ϕ5(x0,0,z,0,y,0,T/). With this information, we expect that if x0 is close enough to the halo orbit, then xn → x* as n→∞. g is given by equation 18.

All the information we need to start Newton’s method is shown above.

From the cylinder theorem, it was possible to find a halo orbit family. Thus, having found a halo orbit and noticing that it has exactly two unit eigenvalues, we can use that as a starting point to move along the cylinder. We use the initial conditions from the previous halo orbit as a starting point to find the next halo orbit at a slightly larger value of x (x coordinate closer to the secondary asteroid). If we find another halo orbit here, we iterate through the process. In this way it was possible to calculate a halo orbit family. The x coordinate step to determine each halo orbit was x = 0.00002.

4.4 Halo Orbits

Figures 12 and 13 illustrate how the halo orbits appear at the tangent bifurcations of the planar orbits around L₁ and L₂ when μ*=10-5 and d = 0. As we built the family of halo orbits, we observed that the amplitude of the orbits increases as the halo orbits move towards the secondary.

Figure 12: Stability index (s₂) of the planar and halo orbit families around L₁ considering d = 0 and μ* = 10^-5. The color figure can be viewed online. 

Figure 13: Stability index (s₂) of the planar and halo orbit families around L₂ considering d = 0 and μ* = 10^-5. The color figure can be viewed online. 

For the conditions considered here, the halo orbit appears at x ≈ 0.98418 for L₁ and at x ≈ 1.01575 for L₂. Figure 14 shows the path of the characteristic multipliers over the unit circle to the halo orbit around L₁ and L₂. Initially, the characteristic multipliers move in the direction shown by the purple arrows until they collide with the negative real axis, configuring a periodic doubling bifurcation. After moving subtly along the real negative axis, the characteristic multipliers return, moving in the direction of the red arrows, colliding again at -1 and then assuming imaginary values, configuring another periodic doubling bifurcation.

Figure 14: Behavior of the of characteristic multipliers at the period doubling bifurcation. The color figure can be viewed online. 

Figures 15 and 16 provide information about the s₁ stability index as a function of d and μ*. Note that the smaller the amplitudes of the halo orbits, the larger the value of the stability index s₁, when considering fixed d and μ*. As the amplitude of the halo orbit increases, the stability index decreases. If we set d = 0, we still detect stable halo orbits for small values of μ*. These orbits were also found by several authors using the restricted three-body problem and are called near rectilinear halo orbits (NRHO) ^Howell(1982), ^{Zimovan-Spreen et al.(2020)}. NRHOs are defined as the subset of the halo orbit family with stability indexes around s _i ± 2 and with no stability index considerably greater in magnitude than the others.

Figure 15: Stability index (s₁) of halo orbits around L₁ as a function of d and μ*. The color figure can be viewed online. 

Figure 16: Stability index (s₁) of halo orbits around L₂ as a function of d and μ*. The color figure can be viewed online. 

Figures 17 and 18 provide information about the stability index s₁ as the size of the dipole increases from 0 to 2000 meters and the mass ratio is kept constant at μ* = 10^-5. The influence of the dimension of the secondary body on the stability of the halo orbits is clear in the plots. Note that, as the size of the secondary increases, the values of s₁ become larger in the vicinity of the equilibrium point L₁ and L₂.

Figure 17: Halo orbit stability index around L₁ for different values of d. The color figure can be viewed online. 

Figure 18: Halo orbit stability index around L₂ for different values of d. The color figure can be viewed online. 

Note that it is unlikely to detect NRHOs around L₁ when we take into account the elongated shape of the secondary body and assume small values of μ*. On the other hand, there are several NRHOs around L₂. In this work, we found NRHOs up to d = 1500 meters, as shown in Figure 18.

However, as shown in ^Howell(1982), the stability index also depends on the mass ratio of the system. Considering d = 0 and increasing μ*, the stability index s₁ increases. We did not detect any NRHO for values of μ* ≥ 10^-1 and d = 0. However, we find NRHO for μ* ≥ 10^-1 and d = 0 around L₂. These results are similar to those obtained by ^Howell(1982). However, taking into account the elongation of the secondary and assuming large values of μ* (μ* ≥ 10^-1), it is possible to find family members of stable halo orbits around L₁ and L₂. Thus, in the model used in this article, stable periodic orbits in the vicinity of irregular bodies exist, even when the secondary has a non-spherical shape. This agrees with the results obtained by ^{Chappaz Howell(2015)}, who found stable orbits around L₁ and L₂ taking into account the elongated shape of the secondary body and considering μ = 0.4 with the triaxial ellipsoid model.

Now we analyze how the period of the halo orbits around L₁ and L₂ is affected by d and μ*. As d increases and μ* is kept constant, the periods of the halo orbits decrease, as shown in Figures 19 and 20.

Figure 19: Period of the halo orbits around L₁ as a function of d and μ*. The color figure can be viewed online. 

Figure 20: Period of the halo orbits around L₂ as a function of d and μ*. The color figure can be viewed online. 

This is because the gravitational attraction is stronger near the particle, due to the mass distribution of the secondary body, causing the acceleration to increase and the orbital period to decrease. As the amplitude of the halo orbit increases, its orbital period becomes shorter.

Considering the elongated shape of the asteroid, but keeping d constant and increasing μ*, we notice that the period of the halo orbits become longer. This is because, as μ* increases, the equilibrium point move away from the secondary body, thus the halo orbits are further away from the secondary body, which causes the gravitational acceleration to decrease, and thus the orbital period of the particle along the orbit to increase.

Figures 21 and 22 provide information on the behavior of the Jacobi constant of the halo orbits as a function of d and μ*. Note that when d or μ* increase, the range of value of the Jacobi constant also increases. This is important information in terms of the application to space missions. Note that the larger the mass ratio of the system, or the longer the secondary body, the less energy is needed for the halo orbits to branch from the planar orbits.

Figure 21: Jacobi constant of the halo orbits around L₁ with respect to d and μ*. The color figure can be viewed online. 

Figure 22: Jacobi constant of the halo orbits around L₂ with respect to d and μ*. The color figure can be viewed online.