Planetary Satellite Orbiters : Applications for the Moon

Low-altitude, near-polar orbits are very desirable as science orbits for missions to planetary satellites, such as the Earth’s Moon. In this paper, we present an analytical theory with numerical simulations to study the orbital motion of lunar low-altitude artificial satellite. We consider the problem of an artificial satellite perturbed by the nonuniform distribution of the mass of the Moon J2–J5, J7, and C22 . The conditions to get frozen orbits are presented. Using an approach that considers the single-averaged problem, we found families of periodic orbits for the problem of an orbiter travelling around the Moon, where frozen orbits valid for long periods of time are found. A comparison between the models for the zonal and tesseral harmonics coefficients is presented.


Introduction
Frozen orbits around the Moon, natural satellites, or asteroids are of current interest because several space missions have the goal of orbiting around such bodies see, for instance, 1-5 and references therein .The dynamics of orbits around a planetary satellite, taking into account the gravitational attraction of a third-body and the nonuniform distribution of mass J 2 , J 3 , and C 22 of the planetary satellite Europa , is studied in 6-11 .This paper contains a continuation of the work developed in 1, 12, 13 , where the problem of the third-body gravitational perturbation and the non-uniform distribution of mass J 2 , C 22 of the planetary satellite Moon on a spacecraft is studied using a simplified dynamical model.Also in 14, 15 , the problem of critical inclination orbits for lunar artificial satellites including the terms J 2 zonal, C 22 sectorial, and the lunar rotation due to the lunar potential, besides the Keplerian term is considered.In 16 , a study done to control the increasing eccentricity of a polar lunar satellite for different altitudes due to the perturbations of the Earth is presented.The approach presented here is the use of low-thrust propulsion in order to keep the orbital eccentricity of the satellite at low values.A fourth-order analytical theory is presented for the accurate computation of quasiperiodic frozen orbits in 17 .
In our research, we use a simplified dynamical model that considers the effects caused by the nonsphericity of the Moon J 2 -J 5 , J 7 , and C 22 using single-averaged analytical model over the short satellite period.The motion of equations is obtained from the Lagrange planetary equations, and then they are integrated numerically, using the software Maple.Emphasis is given to the case of frozen orbits, defined as orbits which have constant values of eccentricity, inclination, and argument of the periapsis when the single-averaged system is considered.An approach is used for the case of frozen orbits similar to the one used in 1, 2, 5 .The basic idea is to eliminate the terms due to the short time periodic motion of the spacecraft to show smooth curves for the evolution of the mean orbital elements for a long time period.The interesting set of results allow us to find good conditions close to frozen orbits.
In other words see 5, 11, 18, 19 , a better understanding of the physical phenomenon can be obtained and it allows the study of long-term stability of the orbits in the presence of disturbances that cause low changes in the orbital elements.
The paper has four sections.Section 2 is devoted to the terms due to the nonsphericity of the Moon and to the Hamiltonian of the system.Applications taking into account the longperiod disturbing potential using the averaging method are presented in Section 3. Section 4 shows the conclusions.

Oblateness of the Moon
To analyze the motion of a lunar low-altitude artificial satellite, it is necessary to take into account the Moon's oblateness 1, 14, 20 .Besides that, the Moon is much less flattened than the Earth but it also causes perturbations in space vehicles.Table 1 shows the orders of magnitude for the zonal and sectorial harmonics for both celestial bodies to make a comparison.The C 20 term describes the equatorial bulge of the Moon, often referred to as oblateness.The coefficient C 22 measures the ellipticity of the equator.
The force function, the negative of the total energy as used in physics, is given by F K V − T here, V is the opposite of the potential energy, and T is the kinetic energy.The force function can be written as K μ/r R − T or K μ 2 /2a R. The function R, comprising all terms of V except the central term, is known as the disturbing potential.The space vehicle is a point of mass in a three-dimensional orbit with orbital elements a semimajor axis , e eccentricity , i inclination , g argument of the periapsis , h longitude of the ascending node , and n mean motion given by the third Kepler law n 2 a 3 Gm 0 , where G is the gravitational constant and r is the radius vector of the body m 0 mass of the planetary satellite .The term due to the unperturbed potential is given by Considering the lunar equatorial plane as the reference plane, the disturbing potential can be written in the form 21 where μ is the lunar gravitational constant, R M is the equatorial radius of the Moon R M 1738 km , P n represent the Legendre polynomials, P nm represent the associated Legendre polynomials, the angle φ is the latitude of the orbit with respect to the equator of the Moon, and the angle λ is the longitude measured from the direction of the longest axis of the Moon, where λ λ − λ 22 , since λ is the longitude measured from any fixed direction, and λ 22 is the longitude of the Moon's longest meridian from the same fixed direction.However, λ 22 contains explicitly the time 20 .Using spherical trigonometry, we have sin φ sin i sin f g , where f is the true anomaly of the satellite.
The following assumptions will be used in this work: a the motion of the Moon is uniform librations are neglected , b the lunar equator lies in the ecliptic we neglect the inclination of about 1.5 • of the lunar equator with respect to the ecliptic and the inclination of the lunar orbit with respect to the ecliptic of about 5 • , c the longitude of the lunar longest meridian is equal to the mean longitude of the Earth librations are neglected , and d the mean longitude of the Earth λ ⊕ is equal to λ 22 .
Since the variables Ω and λ ⊕ appear only as a combination of Ω − λ ⊕ , where λ ⊕ n M t const, n M being the lunar mean motion.The degrees of freedom can be reduced by choosing h Ω − λ ⊕ as a new variable.A new term must then be added to the Hamiltonian in order to get ḣ Ω − n M −∂ K n M H /∂H.This term n M H is used here in the same way it is used in 20 .The Hamiltonian is still time dependent through λ ⊕ .Since the longest meridian is always pointing toward the Earth, it is possible to choose a rotating system whose x-axis passes through this meridian.
Regarding the Earth's potential, the dominant coefficient is J 2 .The rest are of higher order 2 .Opposite to the situation of the Earth, the first harmonics of the lunar potential are almost of the same order see Table 1 .This fact complicates the choice of the harmonic where the potential can be truncated, and this fact makes its choice a little arbitrary.The influences of the Earth and of the nonsphericity of the Moon on the stability of lunar satellites were also shown by 22 , but the sectorial harmonics were not considered.
In terms of the orbital elements, the Legendre-associated functions for the zonal up to J 5 and sectorial C 22 terms can be written in the form 20, 21 where we use the shortcuts ξ cos g cos h−c sin g sin h, χ − sin g cos h−c cos g sin h, s sin i, and c cos i.With this, the zonal perturbation due to the oblateness is 14

2.4
Replacing the relation μ n 2 a 3 and using Cayley's tables 23 to express the true anomaly in terms of the mean anomaly, after some algebraic manipulations, we get e 2 e 2 cos 2l .

2.5
For the sectorial perturbation 20, 21 , we get where δ C 22 R 2 M .With some manipulations, we get where the disturbing potential is written in the form R H 20 H 22 .

The Hamiltonian System
The Hamiltonian of the dynamical system associated with the problem of the orbital motion of the satellite around the Moon taking into account the non-uniform distribution of the mass of the planetary satellite J 2 -J 5 , J 7 , and C 22 can be written in the following form: where Now, using We arrange the Hamiltonian K as follows where

2.13
The motion of the spacecraft is studied under the single-averaged analytical model.The average is taken over the mean motion of the satellite.The standard definition for average used in this research is 24 F 1/2π 2π 0 F dl where l is the mean anomaly of the satellite.The single-average method is applied in 2.13 to eliminate the terms of short period, to analyze the effect of the disturbing potential on the orbital elements.Periodic terms will be calculated, replacing the result in the Lagrange planetary equations 25 and numerically integrated using the Maple software.The long-period disturbing potential is given in the appendix.

Applications: The Long-Period Disturbing Potential Using the Averaging Method
With a simplified model for the disturbing potential, it is possible to analyse the orbital motion of the satellite.
To study the lifetimes of low-altitude artificial satellites of the Moon 26 we took into account the zonal J-J 5 and sectorial C 22 , C 31 terms in the disturbing potential.However, 5 shows that the effect of the J 7 term is an order of magnitude larger than the J 2 -J 5 terms see Table 2, model given by 27 , and therefore it should be considered in the potential.In 5 , an analytic model is presented to find frozen orbits for lunar satellites, where the potential only due to J 2 and J 7 zonal harmonics is taken into account.In this paper, we present an approach taking into account the J 2 -J 5 , J 7 , and C 22 terms.

Frozen Orbits Solutions
The definition of frozen orbit is written in the form dg dt 0, de dt 0, di dt 0 .
Replacing by the potential in A.1 in the Lagrange's planetary equations 25 and solving the equation dg/dt 0, we get an equation that depends on three variables, e, i, and a, that can be represented as a 3-D surface as shown in Figures 1, 2, 3, and 4. Points on this surface correspond to equilibria of the system in equations dg/dt 0 and de/dt 0, that is, frozen orbits.Figures 1 and 3 take into account the coefficients J 2 -J 5 and J 7 , while Figures 2  and 4 take into account the J 2 -J 5 , J 7 , and C 22 terms.Comparing the models given by 27 see Table 2 and 28 see Table 3 , we can analyze that the results are very close and the model given by 28 presents better results to obtain a family of frozen orbits for inclinations above 69 • .In 30 , an analytical theory based on Lie-Deprit transformation is used to locate family of frozen orbits for asteroids and natural satellites.
For inclinations larger than 69 • , we found a family of frozen orbit, as we can see in Figures 1, 2, 3, and 4, mainly considering the model LP165P for the coefficients of the lunar gravity field 28 .We also observed that the C 22 term contributed efficiently to obtain of such orbits.
The curves obtained for constant-a sections are very similar; thus, in order to go into more detail, we define a constant value for the semimajor axis a 1861 km, which indeed  corresponds to low orbits.For the chosen value of a, we plot the curve i versus e see Figures 5 and 6 .Figure 5 does not take into account the C 22 term in the potential, while Figure 6 considers the sectorial term.We observed an increase of the eccentricity for inclinations between 28.65 • and 45.84 • comparing Figures 5 and 6 .For orbits with inclination larger than 69 • degree, we find a family of frozen orbits see Figure 6 .The region where the eccentricity is frozen is for small values of the initial eccentricity, when compared with Figure 5.
Taking into account that Figure 6 represents a family of frozen orbits for inclinations of 90  8 shows frozen orbits for a case where the semimajor axis is 1838 km.We found frozen orbits with initial eccentricity equal to 0.02 that can be applied for missions of long period without the needs for doing orbit corrections for a period of very long time larger than 4000 days .The initial orbit with appropriate characteristics can extend Mathematical Problems in Engineering the time life and reduce maintenance costs.The set of initial conditions for frozen orbits with long lifetime is given by a 1861 km, e 0.02, i 90 • , g 90 • , and h 270 • .
In the literature, in general, frozen polar orbits are found with periapsis g 270 • 3, 31-33 , however, for orbits with inclinations near 90 • , the argument of periapsis moves to g 90 • 2, 31, 34 .These authors take into account only the zonal terms.In 34 , there is an approach that compares the disturbing potential when a model of the gravitational field of high order taking into account the terms J 7 , J 20 , and J 50 is considered.For small inclinations, the three models show similar results for low-altitude orbits see 34 , while for high inclination orbits the J 20 truncation only approximates roughly the J 50 case, giving frozen orbits that, in general,  show lower eccentricities and shifting the zero eccentricity solutions to lower inclinations by several degrees.In this case, the potential truncated in the J 7 term provides a qualitative description only correct in the case of low inclinations.In 35 , a numerical study considering models for the lunar gravity field of high order is applied to near-circular, low-altitude orbits.The behavior of the periapsis altitude is analyzed for some harmonics, for example, J 7 and J 9 .Here, we truncate the potential of the lunar gravity field up to the J 7 term, to analyze the effect caused by this harmonic in particular 5 .When we consider the potential taking into account the terms J 2 -J 5 and C 22 zonal and sectorial , the frozen orbit i 90 • , see Figure 10 librates around the equilibrium point for g 270 • which is in agreement with the literature see 3, 31, 32 .But, when we take into account the potential up to the term, we found frozen orbits that librate around the equilibrium point for g 90 • , which is in agreement with 5 .In 34 , a frozen orbit considering the following initial conditions: a 1861 km, e 0.0036, i 84 • , and g 90 • , is found, and it corresponds to a near-polar orbit.
Then, we develop the disturbing potential up to the J 9 term J 2 -J 9 , C 22 to analyze what happens with the value of g for a low-altitude, near-circular, near-polar orbit.Figure 11 shows the behavior for the diagram e versus g for a polar orbit which librates around the equilibrium point for the value of g 270 • see 33 for more details with respect to the potential truncated up to the J 9 term. .In this case, the initial eccentricity varies from 0.001 to 0.004.The numerical values used for the harmonics up to J 9 is given in Table 3.Thus, we conclude that, because of the characteristics of the lunar gravity field, the potential truncated up to the J 7 produces the effect of the frozen orbit obtained that librates around g 90 • and not g 270 • .
The idea here is to consider a simplified potential to analyze the effect of the harmonics on the orbital elements e and g.Therefore, this study did not realize the full potential integrations, which can be viewed in 33, 34 .Another important factor that we must analyze is the value of the argument of the periapsis, that strongly affects the determination of the frozen orbit.As we previously showed, the values of the periapsis to assure the condition of frozen orbit are g 90 • or g 270 • , equally for the longitude of the ascending node that is h 90 • or h 270 • .Since that, in general, we find frozen orbits for the value of g 90 • .Figures 12-15 show the behavior of the eccentricity as a function of time for different values of the inclinations.Note that Figure 12 presents an orbit with lesser variation of the eccentricity for inclination of i 90 • , considering g 90 • .When we compare Figure 12 with Figure 13, we get that, for the value of g 270 • , the orbit presents great variations of the eccentricity.We also analyze the case for the initial eccentricity of the order of 10 −3 and verify that the behavior of the orbit does not suffer the same effect as occurred for the case of the eccentricity of order 10 −2 .Here, the value of g 90 • or g 270 • has smaller effect on the amplitude of the eccentricity see Figures 14 and 15 .

Conclusions
In this paper, the dynamics of orbits around a planetary satellite Moon , taking into account the non-uniform distribution of mass J 2 -J 9 and C 22 , was studied.Several analyses were performed, using the long-period disturbing potential obtained through the single-averaged model.
Considering the long-period disturbing potential, we found a family of frozen orbits for inclinations larger than 69 • considering the model LP165P for the coefficients of the lunar gravity field.It was observed that the C 22 term contributed efficiently to obtain such orbits.
Low-altitude, near-polar orbits are very desirable for scientific missions to study natural satellites, such as the Moon.Due to the characteristics of the lunar gravity field, the potential truncated up to the J 7 produces the effect of the frozen orbit that librates around g 90 • and not g 270 • .Considering the long-period disturbing potential using the singleaveraged method, we found near-circular, near-polar frozen orbits where this choice reduces the maintenance cost considerably.
In general, we conclude that, for the case of the Moon, it is necessary to consider a gravity model of high order, including the zonal and sectorial terms.Because the question of the coefficients are almost of the same order, this fact complicates the choice of the harmonic where the potential can be truncated and fact makes its choice a little arbitrary.
To study the lifetime of a lunar artificial satellite considering an analytical model, the zonal J 2 -J 9 and sectorial C 22 terms of the harmonic coefficients should be taken into account.We observe that the frozen orbits could be strongly affected by terms of the potential containing the coefficient due to the equatorial ellipticity of the Moon C 22 and by terms containing the argument of the periapsis and the longitude of the ascending node.

Appendix
The long-period disturbing potential up to the J 7 term is k1 75 16a where the values for the zonal J 2 -J 7 and tesseral C 22 harmonics coefficients taking into account the two models are given in Table 2 27 and Table 3 28 .

Figure 2 :
Figure 2: e versus i-rad versus a-km.Akim, 27 , J 2 -J 5 , J 7 and C 22 .Green color is the surface when g 3π/2, and the red is obtained for g π/2.

Table 1 :
Orders of magnitude for J 2 and C 22 14 .
•, an example to illustrate the behavior of the frozen orbits for a case where the semimajor Figure 3: e versus i-rad versus a-km.LP165P, Konopliv et al. 28 , J 2 -J 5 , J 7 .Green color is the surface when g 3π/2, and the red is obtained for g π/2.LP165P, Konopliv et al. 28 , J 2 -J 5 , J 7 and C 22 .Green color is the surface when g 3π/2, and the red is obtained for g π/2.
Figure 4: e versus i-rad versus a-km.