The effect of zonal harmonics on dynamical structures in the circular restricted three-body problem near the secondary body

The circular restricted three-body model is widely used for astrodynamical studies in systems where two major bodies are present. However, this model relies on many simplifications, such as point-mass gravity and planar, circular orbits of the bodies, and limiting its accuracy. In an effort to achieve higher-fidelity results while maintaining the autonomous simplicity of the classic model, we employ zonal harmonic perturbations since they are symmetric about the z-axis, thus bearing no time-dependent terms. In this study, we focus on how these perturbations affect the dynamic environment near the secondary body in real systems. Concise, easily implementable equations for gravitational potential, particle motion, and modified Jacobi constant in the perturbed model are presented. These perturbations cause a change in the normalized mean motion, and two different formulations are addressed for assigning this new value. The shifting of collinear equilibrium points in many real systems due to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_2$$\end{document}J2 of each body is reported, and we study how families of common periodic orbits—Lyapunov, vertical, and southern halo—shift and distort when \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_2$$\end{document}J2, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_4$$\end{document}J4, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_6$$\end{document}J6 of the primary and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$J_2$$\end{document}J2 of the secondary body are accounted for in the Jupiter–Europa and Saturn–Enceladus systems. It is found that these families of periodic orbits change shape, position, and energy, which can lead to dramatically different dynamical behavior in some cases. The primary focus is on moons of the outer planets, many of which have very small odd zonal harmonic terms, or no measured value at all, so while the developed equations are meant for any and all zonal harmonic terms, only even terms are considered in the simulations. Early utilization of this refined CR3BP model in mission design will result in a more smooth transition to full ephemeris model.


Introduction
The circular restricted three-body problem (CR3BP) is a useful platform for understanding and designing trajectories in the presence of two large, perturbing bodies (Szebehely 1967;Koon et al. 2006). Widespread use of the CR3BP has given rise to the field of low-energy trajectory design, which can help minimize fuel usage and allow access to a wide variety of orbital geometries (Howell 1983;Koon et al. 2006;Parker et al. 2013;Vaquero Escribano 2013;Bosanac 2016;Restrepo 2018). However, popular solutions obtained with the classical CR3BP are sometimes not dynamically robust, particularly in sensitive regions. The CR3BP makes large assumptions such as the circular motion of bodies and point-mass gravity. The purpose of this paper is to provide the reader easily implementable tools that can increase the accuracy of the CR3BP dynamic model. This goal is accomplished by modifying the CR3BP to account for zonal harmonic perturbations. These are terms that describe non-spherical properties of a body that happen to be symmetric about the z-axis, such as J 2 , which describes the difference between a body's equatorial and polar radii. Because of their z-axis symmetry, these perturbations are autonomous in the CR3BP, making them ideal candidates for further study. In order to make these perturbations easily implementable, derivations of the CR3BP gravitational potential, equations of motion, and modified Jacobi constant which account for zonal harmonic terms are presented, and the effects of these perturbations on collinear equilibrium points and periodic orbits are studied.
Certain zonal harmonics in the CR3BP have previously been considered. One of the immediate effects of accounting for zonal harmonics is a change in the system potential, which can affect the positions of the Lagrange equilibrium points. For a number of threebody environments in the Solar system,  study the locations of all five equilibrium points when J 2 of the primary body is accounted for. Markellos et al. (2004) also consider J 2 of the primary, but study homoclinic and heteroclinic connections between Lyapunov orbits in the Hill problem. Jain and Aggarwal (2015) study non-collinear equilibria accounting for J 2 of the primary, but also account for a dissipative Stokes drag force from the primary body. Carvalho et al. (2014) study the design of near-circular frozen orbits when considering the third-body perturbations of Jupiter along with J 2 , J 3 , and C 22 of Europa. Zotos (2015) categorizes a wide range of starting conditions in a three-body system as being either bounded, escaping, or collisional trajectories, and studies how these categorizations change when accounting for J 2 of the primary. Cinelli et al. (2015) studied the use of polynomial equations to account for zonal harmonics of Europa up to J 4 and a second-order Legendre polynomial expansion to account for Jupiter's third-body perturbations in order to find sunsynchronous orbits about Europa. Singh and Umar (2017) study periodic orbits around L 4 and L 5 when both bodies are oblate in the ER3BP. Cinelli et al. (2019) looked at longlifetime, high-inclination orbits about Europa when considering J 2 and C 22 of Europa via a double-averaged expression of disturbing potential.
There are a number of studies that look at the location and stability of the equilibrium points when J 2 of both the primary and secondary bodies are accounted for , Bhatnagar and Hallan (1979), Elshaboury (1989), Arredondo et al. (2012). In addition to analyses of equilibrium points, some studies look at the effects of zonal harmonics on planar periodic orbits. Certain periodic orbits have been shown to exist when J 2 is taken into account for only the primary body (Sharma 1972a, b, c), for only the secondary body (Mittal et al. 2009), and for both (Abouelmagd et al. 2013). Although most work in this field looks specifically at the J 2 term, some studies also consider the effects of J 4 . Abouelmagd et al. (2015a) verify the existence of equilibrium points, find their locations, and study their stability when both J 2 and J 4 of the primary and secondary bodies are accounted for. These authors also show that the L 4 and L 5 equilibrium points have planar periodic orbits when J 2 and J 4 of the primary body are accounted for Abouelmagd et al. (2015b).
This study aims to bridge results from the previous literature and simultaneously fill in gaps. This is accomplished by developing concise, easily implementable equations of motion which are able to account for any zonal harmonic perturbations due to both the primary and secondary bodies in the CR3BP. The normalized mean motion of a system is also affected by the added perturbations, and two formulations are addressed for determining this new value. Locations of collinear equilibria are studied under these dynamics, and the effects of these terms are studied near the secondary bodies in the Jupiter-Europa and Saturn-Enceladus systems. Lyapunov, vertical, and halo families of periodic orbits are studied near the L 2 points in these systems. The primary effects of the selected zonal harmonics on these families are a combination of an x-axis shift of the family, distortions in shape, and changes in specific energy at a given time period. In this study, only the effects of even zonal harmonic terms are analyzed because the main focus is on the outer planets, many of which have little or no measured value for odd zonal harmonic terms.

Analytical expressions for the inclusion of zonal harmonics in the CR3BP
Closely following the procedure presented in Roy's Orbital Motion (Roy 2005) for deriving the classical CR3BP equations of motion, this section presents a method for also expressing zonal harmonic perturbations from both the primary and secondary bodies in the CR3BP equations of motion, gravitational potential, and modified Jacobi constant. The derivation starts by including terms for J 2 , J 3 , and J 4 of each body, and works toward general, recursive forms. It is important to note that by including these terms which do not vary with time, it is assumed that the secondary body is orbiting in the equatorial plane of the primary and that the spin axis of each body is aligned. This is a reasonable assumption for many three-body systems. Europa and Enceladus, for example, are each inclined at less than 1 • from their primary's equatorial plane and are also tidally locked.
To begin the derivation for the equations of motion, the normalization of the system and basic inertial frame are defined in Eqs. (1)-(2) where M 1 and M 2 are the masses of the primary and secondary body. The model is illustrated in Fig. 1. The velocities and positions of the bodies remain in thex-ŷ plane, so theẑ andζ vectors remain constant and parallel. Normalization of the CR3BP system First, we define the coordinates of the spacecraft and each body in the inertial frame originating at the barycenter. In this formulation, the variables x, y, and z are reserved for coordinates in the rotating frame, while ξ , η, and ζ are the inertial frame equivalents.
(2) Fig. 1 Layout of the circular restricted three-body problem, showing the relationship between the bodies and the inertial (black) and rotational (red) frames From the recursive form of the equation for gravitational potential with zonal harmonics considered in Eq. (3) (Schaub and Junkins 2013), the accelerations due to both a primary and secondary body can be derived. The results of this derivation for the zonal harmonic terms of J 2 -J 4 are expanded for clarity in Eqs. (4)-(7). To clarify some of the variables used in this section, with respect to the body being considered, G M refers to the gravitational parameter, R is the radius, r is the distance from the body to the spacecraft, J k is the kth zonal harmonic term, and P k ( ζ r ) is the kth Legendre polynomial evaluated with ζ r .
In Eq. (4), the equations of motion for the CR3BP derived without zonal harmonics are shown.ξ Now the expressions are condensed with a set of γ and β variables. In this notation, γ 2 p and γ 2s correspond to the expressions for the effects of J 2 p and J 2s , respectively. The β terms follow the same notations, but while γ is used for x and y terms, β is used for z terms.
Substituting into Eqs. (5), (6), and (7), With these substitutions, the inertial equations of motion are concise: To relate the inertial and rotating frames, we start by presenting relevant vectors in the rotating frame, which can all be seen in Fig. 1: We now describe the inertial coordinates in terms of the rotating frame, which has rotated at an angle of nt, where n is the mean motion of the system: Differentiating one has: ξ =ẋ cos nt − nx sin nt −ẏ sin nt − ny cos nt, η =ẋ sin nt + nx cos nt +ẏ cos nt − ny sin nt, Now (15) is substituted into (13) and the sin and cos terms are collected. Theζ term is temporarily ignored due to its trivial conversion between inertial and rotational frames. See Eq. (17):ξ Equating (19) with (17), By multiplying (20) by cos nt and (21) by sin nt, adding, then separately multiplying (20) by − sin nt and (21) by cos nt, and adding again, the following relation emerges 1 By substituting in values for x 1 and x 2 from (14) and rearranging, we get the final equations of motion for the CR3BP model which includes J 2 , J 3 , and J 4 of each body:

Recursive equations of motion, gravitational potential, and Jacobi constant
From (23), there is a clear pattern which allows the full equations to be written in recursive where With this form, the full gravitational potential for the CR3BP with zonal harmonics is clear: Multiplying these three equations byẋ,ẏ, andż, respectively, and adding, results in: Integrating this, an expression is found where C is a constant of integration and V is the magnitude of the particle's velocity in the rotating frame. This C can be thought of as a modified Jacobi constant:

Normalized mean motion in the zonal harmonics perturbed CR3BP
Let us denote the Keplerian mean motion as Following the formula for Keplerian mean motion [Eq. (31)] and plugging in the standard normalization values [Eq.
(1)], we see that in the classical formulation of the CR3BP, the normalized mean motion is equal to 1. However, when zonal harmonic perturbations are considered, this can no longer be assumed and a new formulation must be adopted. There are at least two approaches for coming up with this value. These approaches will be referred to as the "theory-based" and "ephemeris-based" methods, and each has advantages and disadvantages. The ephemeris-based method is used to generate the numerical results for this study. For clarity, a subscript 'n' will be used to denote a normalized, dimensionless value (e.g., T P n = 2π, classically).

Theory-based method (example with J 2p and J 2s )
For the theory-based method, the orbit of the primary-secondary system is treated as circular.
The expression for mean motion is then determined by defining and equating the dynamical and kinematic expressions forr 12 between the two gravitational bodies. To find the dynamic expression forr 12 , set up a system with two gravitational bodies and an inertial origin, define r 1 and r 2 as vectors from the origin to the respective bodies, find the expressions forr 1 and r 2 based on Newton's laws and zonal harmonic potentials [Eq.
(3)], and difference the two to findr 12 . A similar process has been seen in much of the literature regarding the threebody problem with zonal harmonics Markellos et al. 2004;Jain and Aggarwal 2015;Zotos 2015;Singh and Umar 2017;Abouelmagd et al. 2013Abouelmagd et al. , 2015a.
In the unperturbed problem, this would show up as the recognizable: However, when we include J 2 p and J 2s , we find the following relationship: Next, we want to solve for the kinematic expression forr 12 , which comes from a time derivative of the position vector: Finally, by relating the dynamic and kinematic expressions forr 12 , an analytical expression for mean motion is obtained:

Ephemeris-based method
The key to this method is that no attention is paid to specific perturbations which are included in the model. Instead, a time-normalizing constant, t N , is formulated to ensure that G(m 1 + m 2 ) = 1. This t N can then be used to normalize the real, measured value of mean motion.
To start, we recognize that the units of G M are km 3 s 2 . This relationship can then be exploited to find t N : Using t N to normalize the measured value of mean motion, we arrive at the final expression, which in general does not equal 1.

Comparison of theory-based and ephemeris-based methods
Following the theory-based method, the calculated n n value is consistent with the exact model being used. However, the analytical process can be laborious as a different formulation is required for any combination of perturbations being considered. Further, there are still assumptions being made in this formulation that do not describe real-world behavior, such as the fact that the full-body problem is not being addressed (i.e., the effects of J 2 p on J 2s and vice versa). Theoretically, it may be that the n n value calculated with the theory-based method converges toward the ephemeris-based method as perturbations are added. An advantage to the ephemeris-based method is that its formulation is independent of the perturbations being included in the model, so the process always remains the same. Since the measured mean motion is a function of all the real perturbations, when we normalize this value, all the perturbations are still being represented. It is possible that solutions produced from this method may translate more easily into an ephemeris model, although this has not been investigated. Along with this, the ephemeris-based method may provide the most flexible way to use real data while still keeping a simple model for numerical studies.
The main drawback of the ephemeris method is that it can overcompensate for perturbations-the normalized mean motion value is based off of all real perturbations, but the simplified dynamical model being used in the study might include very few of them.
To compare the effects of the two methods of obtaining normalized mean motion, Tables 1 and 2 present the shifting of the collinear equilibrium points when the normalized mean motion value is changed according to the theory-based method (Table 1) and the ephemerisbased method (Table 2), but no perturbations are added to the dynamic model itself. The point of this is to isolate the effects of a change in normalized mean motion. For the theorybased case, the normalized mean motion is calculated via Eq. (37)-that is, J 2 p and J 2s are accounted for.
The results show that an increase in n n acts to move the collinear points closer to the origin, and a lower n n moves the points further away. The magnitude of shifting is larger in the ephemeris-based case for all systems except Saturn-Enceladus. An important distinction to be aware of when comparing equilibria shifting due to changes in modeled potential versus changes in normalized mean motion is that the former is a change occurring within the frame, while the latter is a change of the frame itself, so we are observing the results of very different processes. In the remainder of the study, the ephemeris-based method for calculating n n is used to simplify the process when various zonal harmonic terms are being considered.

Effects on dynamic environment
Before studying the dynamical effects of zonal harmonic perturbations in the CR3BP, we present the system parameters that were chosen for this study. Tables 3 and 4 provide the chosen parameters relating to primary and secondary bodies, respectively. In Table 4, the semimajor axis, a, and the mean motion, n, are shown as these terms are necessary for system normalization.
When terms for zonal harmonics are included in the CR3BP, the gravitational potential of the system grows and the normalized mean motion changes. As a result, the collinear equilibrium points, L 1 , L 2 , and L 3 , shift from their nominal locations. Knowing precisely where these equilibria are located is important for low-energy trajectory design, especially when looking to make use of dynamical systems theory to study periodic orbits and their manifolds. Table 5 provides the shifting seen in the collinear equilibria for real three-body systems when J 2 of both the primary and secondary is accounted for, and the ephemeris-based method is used to determine n n [Eq. (40)].
As the system dynamics have changed and the Lagrange points have shifted, we can also expect trajectories to behave differently. One example of this in Fig. 2 shows a set of initial conditions corresponding to an unstable manifold from a northern halo orbit at Enceladus that, under standard CR3BP dynamics, goes on to impact the south pole at a near-tangent angle after several orbits. However, as we see in blue, when the same conditions are propagated with J 2 p , J 4 p , J 6 p , and J 2s accounted for, the trajectory quickly diverges from the nominal path and impacts Enceladus near the equator in less than one orbit. Clearly such differences could necessitate redesigning a targeted trajectory. To further study the effects of these zonal harmonic perturbations on trajectories, we examine shifts and distortions in common families of periodic orbits.

Effect on periodic orbits
Periodic orbits are inherently sensitive to initial conditions, so we can expect initial conditions computed in the classical CR3BP to not correspond to periodic orbits in the modified system. However, families can still exist, though they may shift and undergo changes in shape. Figures  3, 4, and 6 depict families of L 2 Lyapunov, vertical, and southern halo orbits at both Europa and Enceladus, with and without the inclusion of zonal harmonics. In each figure, 10 orbits from each model are shown. The 10 orbits correspond to matching time periods-thus, the largest blue and red orbits in each figure have the same period. The zonal harmonic terms included here are J 2 p , J 4 p , J 6 p , and J 2s . It is clear that the effects of the chosen zonal harmonic terms have a more dramatic relative effect at Enceladus. A common theme among the newly computed families is an x-axis shift corresponding to the magnitude of the shift of the L 2 equilibrium point. Less intuitive is the change in shape exhibited by some families. For the case of Lyapunov orbits, we see a shortening of the family, or a reduction in y-axis amplitude at Enceladus, and the opposite at Europa. The magnitude of this effect is relatively greater at Enceladus. Changes to the vertical orbit family are not easily seen. The tall and skinny nature of vertical orbits makes them difficult for figures of this nature hoping to highlight differences between models, so the vertical orbits in Fig. 4 are shown with unequal axis units, hence the distorted shapes of Europa and Enceladus. Other than the apparent equilibrium point shift, it is not obvious what shape changes are occurring. One way this can be studied is by observing the change in z-axis amplitude between the two models for a given time period. This relationship is shown for Europa and Enceladus in Fig. 5. At Europa, we see that the effect of the zonal harmonics is to increase the amplitude of all vertical orbits in this range of time periods, but especially those that have particularly short or long periods. At Enceladus, we see that the zonal harmonics generate shorter vertical orbits at smaller periods and taller vertical orbits at larger periods. When comparing normalized distances between these two systems, it is important to remember that the normalized distance of 1 × 10 −3 is equal to roughly 671 km at Europa and 238 km at Enceladus (Fig. 6).
Looking to the southern Halo families at Europa and Enceladus, a different trend emerges. Again by comparing orbits of the same time period, the included zonal harmonic terms decrease both the y-and z-axis ranges at Europa, while both are enlarged at Enceladus. An unexpected result is that unlike how vertical and Lyapunov families shift along the x-axis in Fig. 3 Lyapunov orbits computed with and without zonal harmonics accounted for, shown at Europa and Enceladus. In each system, 10 orbits are shown which have the same 10 time periods, spaced evenly from 3.09 to 5.51 at Europa, and 3.11 to 4.29 at Enceladus the same direction as does the Lagrange point, the southern halo families shift in the opposite direction.
To observe the effect of zonal harmonic perturbations on the energy of these periodic orbit families, we can study Figs. 7,8,and 9 (Lyapunov,vertical,and southern halo,respectively). The x-axis values are time periods of periodic orbits since this value is monotonic across these families, making it a convenient parameter with which to scan a family. The y-axis of these figures is approximately a difference of specific energy between the two models. To determine the y-axis value, Jacobi constants are computed across a family in each model. Since the equations for Jacobi constant in the classical and perturbed models are fundamentally different, it is not valuable to compare specific values between the two. This can be addressed by differencing a Jacobi value with the Jacobi constant of L 2 in either model. This J C between each periodic orbit and L 2 can be roughly translated as a specific energy, or −v 2 [Eq. (30)], which can be compared between models as long as the points of comparison are near each other in three-dimensional physical space. (Clearly from Figs. 3, 4, and 6, the families are relatively closer to each other at Europa than Enceladus, so larger errors may be present at the latter.) We subtract the J C values computed for the classical CR3BP from those computed for the perturbed model. This is the relationship shown in the figures, and it describes approximate changes in energy to a family caused by the addition of zonal harmonics. At a specific time period on the x-axis, a negative value on the y-axis implies that the Fig. 4 Vertical orbits computed with and without zonal harmonics accounted for, shown at Europa and Enceladus. In each system, 10 orbits are shown which have the same 10 time periods, spaced evenly from 3.22 to 6.18 at Europa, and 3.24 to 6.0 at Enceladus. Note-for the sake of more clearly showing the relative shape of the orbit families, axes shown are not equal in distance Southern Halo orbits computed with and without zonal harmonics accounted for, shown at Europa and Enceladus. In each system, 10 orbits are shown which have the same 10 time periods, space evenly from 1.81 to 3.114 at Europa, and 2.30 to 3.089 at Enceladus zonal harmonic perturbations have decreased the Jacobi constant, and therefore increased the specific energy. The figures reveal a complex relationship. For Lyapunov orbits with common time periods within the given bounds, the included zonal harmonics strictly raise the energy at Europa and decrease the energy at Enceladus. In each case, the magnitude of change decreases as time period is increased. For vertical orbits, energy strictly increases at Europa, but the sign of energy change is dependent on time period at Enceladus. Here, most orbits lose energy, but at the large end of the time period spectrum, orbits start to gain energy. Interestingly, we see from Fig. 5 that around T p = 5.6 at Enceladus, there is a point where zonal harmonics do not cause a change in z-axis amplitude, yet from Fig. 8, we see that at this point, specific energy is still changing. Looking at the southern halo results, the two systems have inverse results. At Europa, short-period halo orbits gain energy, but switch to losing energy at an increasing rate as time period grows. At Enceladus, short-period halo orbits lose energy, but quickly switch to gaining energy at an increasing rate.

Conclusion
The equations of motion, gravitational potential, and modified Jacobi constant are all provided for the inclusion of zonal harmonics of each body in the circular restricted three-body problem. Specific equations of motion are provided for zonal harmonics through J 4 , and recursive forms are provided so the reader can easily determine necessary equations for any choice of zonal harmonics.
The effect of J 2 of each body is studied by determining how the collinear equilibrium points shift along the x-axis for a variety of real moons in the Solar system. Determining these shifts is a non-trivial process that involves the understanding of two core pieces-the Fig. 7 Approximate change in specific energy due to zonal harmonic perturbations for families of L 2 Lyapunov orbits at Jupiter-Europa (left) and Saturn-Enceladus (right) Fig. 8 Approximate change in specific energy due to zonal harmonic perturbations for families of L 2 vertical orbits at Jupiter-Europa (left) and Saturn-Enceladus (right) Fig. 9 Approximate change in specific energy due to zonal harmonic perturbations for families of L 2 southern halo orbits at Jupiter-Europa (left) and Saturn-Enceladus (right) effects of the change in system potential and the effects of a change in normalized mean motion. By accounting for J 2 and any other zonal harmonics, the potential of the system is raised, and the result is that collinear equilibrium points shift away from the barycenter (when not accounting for changes in n n ). The effect that the normalized mean motion has on the collinear equilibrium points is dependent on its magnitude. If the perturbed normalized mean motion is larger than 1, its effect is to shift the collinear points toward the barycenter. However, if this value is less than 1, its effect is to shift the collinear points away from the barycenter. If the theory-based approach is used for calculating normalized mean motion, and only zonal harmonic terms are included in the model, then n n must be greater than 1. The ephemeris-based method is used in the generating of periodic orbits in this study.
To help study the real effects that these perturbations might have on mission design to the outer Solar system, families of Lyapunov, vertical, and southern halo periodic orbits are studied near the L 2 point in the Jupiter-Europa and Saturn-Enceladus systems. Families are located both in the classical model and in the model perturbed by J 2 p , J 4 p , J 6 p , and J 2s . Odd zonal harmonic terms are ignored, since their values are very small or unmeasured at Jupiter and Saturn. These families are then compared against each other by viewing orbits with 10 common time periods, which make the relative shifting clear. It appears that L 2 Lyapunov and vertical orbits generally shift in the same direction as the L 2 equilibrium point when these zonal harmonic perturbations are considered, while the southern halo family appears to shift in the opposite direction as the equilibrium point. The families also undergo changes in shape and specific energy that are not consistent as a function of time period. The changes to shape can distort the periodic orbits by hundreds of kilometers.
The relationships found in this paper make it clear that the inclusion of zonal harmonic perturbations has the potential to greatly change early mission designs which are based on certain geometries or time periods, although the magnitudes of change caused by zonal harmonics are system-dependent. This result is useful for saving time in iterative mission design when solutions come to reckon with a full ephemeris model which accounts for such perturbations. The early inclusion of zonal harmonic terms can smooth the planning process and lead to a much more realistic solution.