Global searches of frozen orbits around an oblate Earth-like planet

A frozen orbit is beneficial for observation owing to its stationary apsidal line. The traditional gravitational field model of frozen orbits only considers the main zonal harmonic terms J2 and limited high-order terms, which cannot meet the stringent demands of all missions. In this study, the gravitational field is expanded to J15 terms and the Hamiltonian canonical form described by the Delaunay variables is used. The zonal harmonic coefficients of the Earth are chosen as the sample. Short-periodic terms are eliminated based on the Hori–Lie transformation. An algorithm is developed to solve all equilibrium points of the Hamiltonian function. A stable frozen orbit with an argument of perigee that equals neither 90◦ nor 270◦ is first reported in this paper. The local stability and topology of the equilibrium points are obtained from their eigenvalues. The bifurcations of the equilibrium points are presented by drawing their global long-term evolution of frozen orbits and their orbital periods. The relationship between the terms of the gravitational field and number of frozen points is addressed to explain why only limited frozen orbits are found in the low-order term case. The analytical results can be applied to other Earth-like planets and asteroids.


Introduction
Considering the gravity of non-spherical perturbations on a satellite, both precession of the orbital plane and apsidal line rotation in the inertial space exist. However, a special case called a frozen orbit occurs where the rotation of the apsidal line stops [1,2].
Frozen orbits are generated by non-spherical gravitational perturbations. Because of the importance of satellite orbit perturbations in the design of a satellite orbit in space, satellite orbit perturbation theory to analytically solve the real orbit has been developed and improved for more than 70 years. Since the 1950s, several well-known scholars, e.g., Brouwer [3] and Kozai [4], have studied orbital evolution under zonal harmonic perturbations. Brouwer used Delaunay variables instead of orbital elements to transform the Lagrangian planetary perturbation equation into its Von Zeipel form and then divided the perturbation terms into secular, long-periodic, and short-periodic [3]. Under the effect of J 4 zonal harmonics, Kozai applied the double-averaging method to analytical solutions of the mean orbital elements [4]. Musen improved Hansen's method to extend solutions of any order under zonal harmonic perturbations [5]. Vinti derived a closed-form analytical solution under the influence of oblateness by replacing conventional circular coordinates with rotating ellipsoidal coordinates in Ref. [6]. Kaula developed a second-order solution under both zonal and tesseral harmonic perturbations using mean orbital elements [7]. Lorell et al. used averaging to develop the second-order long-periodic and secular effects of oblateness [8]. To avoid the difficulties caused by small eccentricity or inclination, Poincaré variables instead of Delaunay variables were used to improve Brouwer's results in Ref. [9]. Hori [10] and Deprit [11] employed the Lie series method in the canonical transform to simplify the canonical perturbation equation following xuming@buaa.edu.cn Nomenclature U potential of spacecraft subject to the attraction of an axisymmetric body F Hamiltonian function µ gravitational parameter of the oblate body R disturbing potential ae mean equatorial radius of the oblate body r distance of the spacecraft from the center of the body ϕ latitude of the spacecraft Jn zonal harmonics coefficients of degree n Cmn, Smn tesseral harmonics coefficients of degree m, n Pn(x) Legendre polynomials Pmn(x) associated Legendre polynomials a semi-major axis e eccentricity i inclination ω argument of perigee Ω right ascension of the ascending node (RAAN) M mean anomaly and its derivative in terms of time λ geographical longitude f true anomaly L, l, G, g, H, h Delaunay variables the approach of Brouwer [3]. Barrar improved Vinti's work [6] to achieve an analytical solution with the same accuracy, yet in a simpler form [12]. The majority of the methods for frozen orbits have been developed along with orbital perturbation theory. The concept of a frozen orbit was proposed in the SEASAT-A mission by Cutting in 1977 [13]. Because its eccentricity and argument of perigee remain constant on average, the frozen orbit has been applied to high-latitude communication missions and remote sensing satellites around the Earth and other celestial bodies. Currently, frozen orbits can be divided into two types. The first type of frozen orbit considers only the gravity field up to the J 2 perturbation term, which has been typically adopted in other studies [14,15]. The inclination angle of this frozen orbit is either 63.43 • or 116.57 • ; this angle is also called the critical inclination. The Russian satellite Molniya [16], which was conducted in 1965, is representative of this type of frozen orbit. The other type of frozen orbit considers higher zonal harmonic terms. The argument of perigee is 90 • or 270 • with small eccentricity (1 × 10 −3 ) in this type of frozen orbit, which is represented by the SEASAT, LandSat [17], and GeoSat [18] conducted by the United States in 1978, 1972, and 1985, respectively. The latter type of frozen orbit is the primary focus of this study.
Numerous scholars have studied frozen orbits using different methods. Coffey used the Lie transformation to derive the normalized Hamiltonian canonical equation at J 9 zonal harmonics and derived families of frozen orbits around an Earth-like planet [19]. The bifurcation in the study by Coffey et al. is similar to that in this study. However, certain equilibrium points were ignored owing to the limitations of the low-order gravitational field. Moreover, the global distribution and evolution of the equilibrium points were not displayed in that study. Aorpimai et al. derived an analytical solution of Earth frozen orbits with arbitrary terms of zonal harmonics by avoiding the singularity of the equation through epicycle elements [20]. In addition, the study of frozen orbits has been extended to other celestial bodies such as the Moon and Mars. Liu et al. proposed a frozen condition and analyzed the differences in frozen orbits around the Earth and Moon based on the different coefficients of the gravitational field model between different celestial bodies [21]. Abad et al. [22] used the same method as San-Juan et al. [23] and applied it to frozen orbits for a lunar orbiter. Liu et al. found four families of Martian frozen orbits, three of which were orbits with small eccentricity, using the Lie transformation [24,25]. Wang et al. proposed several control strategies to realize a Sun-synchronous frozen orbit with arbitrary orbital elements using continuous low-thrust [26]. Lara used a methodology based on the construction of both eccentricity vector diagrams and inclination-eccentricity diagrams of frozen orbits to disclose the qualitative and quantitative differences introduced in the Earth's orbit behavior by truncation to different degrees of the Legendre polynomial expansion of the zonal geopotential [27] . The goal of this study is to explore the additional properties of frozen orbits and identify unreported frozen orbits that could not be obtained in the low-order case. Thus, the gravitational model is extended to the J 15 term described in Delaunay variables, which are used as the main variables to simplify the calculation; orbital elements are adopted to better present the properties of the frozen orbits. Short-periodic terms are eliminated to accelerate the long-term evolution of the frozen orbits through a Hamiltonian canonical transformation based on the Lie series. The equilibrium points of dynamical systems after the Hori-Lie transformation can also be referred to as frozen orbits. Only equilibrium points (or frozen orbits) with the argument of perigee ω equal to either 90 • or 270 • have been widely studied in previous reports under low-order gravitational fields; these points are defined as traditional equilibrium points (or traditional frozen orbits) in this study. The other equilibrium points, whose argument of perigee ω is not equal to 90 • or 270 • , are defined as nontraditional equilibrium points or frozen orbits; they can be considered as the bifurcations of the traditional orbits presented by their global evolutions. The relationship between the terms of the gravitational field and number of frozen points is used to explain why fewer frozen orbits are found in lower-order cases.
2 High-order gravitational field reduced by Hori-Lie transformation The main perturbation of a satellite near an Earthlike planet is caused by the non-central force of the gravitational field. The gravitational potential function is typically written as [28]: J n a e r n P n (sin ϕ) n m=2 a e r n P mn (sin ϕ)(C mn cos mλ + S mn sin mλ) where P n (x) are Legendre polynomials, P mn (x) are the associated Legendre polynomials, a e is the equatorial radius, ϕ is the geocentric latitude, λ is the geographical longitude, and r is the radial distance. The terms involving J n are zonal harmonics. The terms involving C mn and S mn are tesseral harmonics. Then, R is defined as the disturbing potential, i.e., R = U − µ r . Because the measure of a non-recursive orbit is larger than that of a recursive orbit, the non-recursive orbit is mainly considered in this study. For a non-recursive orbit, the effects of the tesseral harmonic terms are primarily offset, on average, by the rotation of the Earth-like planet. In addition, the tesseral harmonic terms do not produce resonance terms with significantly enhanced perturbation effects on low-orbit satellites. Therefore, only zonal harmonics J 2 -J 15 are considered in the disturbing potential. To calculate the frozen orbit, an analytic derivation is required to eliminate the short-periodic terms. Keplerian elements are replaced by Delaunay variables to describe the orbital motion. The relationship between the Delaunay variables and Keplerian elements is as Eq. (2): where a is the semi-major axis, ω is the argument of perigee, Ω is the right ascension of the ascending node, e is the eccentricity, i is the inclination, and M is the mean anomaly. The units of length, time, and mass are meter, second, and kilogram, respectively. Thus, the units of the Delaunay variables L, G, and H are m 2 /s. For simplicity, the units of the Delaunay variables are omitted in this study.
The perturbed equations of motion can be expressed in canonical form:Ẋ where and Hamiltonian function F is described by Based on the relationships: sin ϕ = sin i sin(g + f ) and r = a(1−e 2 ) 1+e cos f , where f is the true anomaly, the potential R can be expanded in the Delaunay variables. Most terms involving f are ignored because they can be eliminated in the following derivation process. Because J 3 -J 15 of an Earth-like planet are of the same order as J 2 2 , the Hamiltonian function F can be expressed as The short-periodic perturbation is eliminated through canonical transformation based on the Lie series, which is also known as the Hori-Lie transformation. After the Hori-Lie transformation, the new Hamiltonian function F * is given by where the symbol {} is the Poisson bracket, F 1s and F 2s are the averages of F 1 (t * ) and F 2 (t * ) about t * , and The detailed analytical expression of the new Hamiltonian function F * is presented in the Electronic Supplementary Material. The variables after this paragraph are all mean variables after applying the Hori-Lie transformation.
After the Hori-Lie transformation, the new Hamiltonian function F * based on the mean elements does not have short-periodic terms. The angular variables l and h do not exist in the new Hamiltonian function F * . Based on Eq. (3), it can be proven that L = const and H = const, where L and H are the mean elements after the canonical transformation. Thus, the order of the new canonical equation is reduced.
Frozen orbits are traditionally defined as those where both the eccentricity and argument of perigee remain constant on average [29]. Thus, the frozen conditions areė = 0,ω = 0,i = 0, which are equal toĠ = 0,ġ = 0, H = const in the Delaunay variables. Because shortperiodic perturbation is not considered in frozen orbits, the variables above are the mean elements after canonical transformation. In this study, the differential equation: is chosen as the main nonlinear dynamical differential equation, and the angular variables l and h are not considered. Based on Eq. (2): H = L · √ 1 − e 2 · cos i, L is assumed to be constant. Then, H is considered an important parameter in the analysis of frozen orbit dynamics.
Based on the nonlinear differential equations in Eq. (7), the Delaunay variables G and g are independent variables. The equilibrium points are calculated using numerical methods. In addition, phase portraits are drawn through the following steps: (1) select the initial points (G 0 , g 0 ) in their domain (H < G 0 < L, 0 < g 0 < 2π) and (2) integrate the differential equation in Eq. (7) using the Runge-Kutta method to obtain the corresponding trajectories. The differences between different phase portraits are changes in the Delaunay variables H and terms considered in the gravitational model.

Expansion of traditional frozen orbits in high-order gravitational field
The earliest frozen orbit, which was designed based on the perturbation model of the gravitational field, is accurate for J 2 and J 3 . With the improved terms of the gravitational field model, the order of Eq. (7) increases. Thus, the solutions to the differential equation, which are also called frozen points, increase in number. In traditional frozen orbits, the Delaunay variable g, which corresponds to the argument of perigee ω, is equal to either 90 • or 270 • . In fact, we can make dG dt and dg dt equal zero by adjusting the values of Delaunay variables G and g at a given H and L. Thus, there are equilibrium points whose g values are equal to neither 90 • nor 270 • . However, such frozen orbits have rarely been mentioned in other literature reports because of the limitation of the terms of the gravitational field model adopted in their work. This type of equilibrium point is defined as a nontraditional equilibrium point. Only traditional frozen orbits are discussed in this section; nontraditional frozen orbits are investigated in the next section.

Expansion and distribution of traditional frozen orbits
Based on the analytical expressions of dG dt and dg dt , dG dt = 0 can be satisfied by setting ω to either 90 • or 270 • , which is considered a basic condition of traditional frozen orbits. Based on Eq. (7), the remaining frozen condition dg dt = 0 can be satisfied by changing the value of G when L and H are given. Because of the gravitational field model with higher terms derived herein, more traditional equilibrium points are obtained in this study.
Based on the analysis above, the global distribution of the frozen orbits is determined by Delaunay elements G and H when L is constant. H ranges from zero to L to include the corresponding G that can satisfy the frozen condition. The relationship between the frozen elements G and H after the numerical calculation is displayed in Fig. 1. G and H are transformed into orbital elements e and i after the numerical calculation. The new relationship is displayed in Fig. 2. Every point in Figs. 1 and 2 corresponds to one frozen orbit.
Each H corresponds to multiple equilibrium points and the number of equilibrium points decreases with an increase in H, as indicated in Fig. 1. Based on H = L · √ 1 − e 2 ·cos i. and H = G·cos i, the ranges of i and G can be obtained by 0 < i < arccos H L and H < G < L. With an increase in H, the ranges of i and G decrease gradually; thus, the number of equilibrium points is reduced. There is a maximum of 19 equilibrium points, with at least two points corresponding to each H. Figure 2 indicates that there are typically frozen orbits with an inclination i ranging from 0 • to 90 • when the argument of perigee is equal to 90 • or 270 • . More precisely, there is a maximum of four frozen eccentricity values for a certain inclination when the argument of perigee is 90 • , and there is a maximum of two eccentricity values for each angle of inclination when the argument of perigee is 270 • .
With the exception of the frozen orbits near the critical inclination, the frozen orbits corresponding to the equilibrium points in Fig. 2 can be divided into quasi-circular orbits and large elliptic orbits based on the eccentricity. First, the quasi-circular frozen orbits with eccentricity less than 0.01 correspond to those in the bottom of Fig. 2, which are typically studied in the majority of reports. Secondly, the eccentricity values of large elliptic frozen orbits are greater than 0.5 and these orbits supplement the quasi-circular frozen orbits because of the higher-order terms of the gravitational field model derived in this study; these terms are rarely studied in other reports. Frozen orbits near the critical inclination are the link between the quasi-circular orbits and large elliptic orbits. Thus, it is implied that the high-order terms of the gravitational field cannot be neglected for orbits with large eccentricity [19,30]. This phenomenon reflects the complexity of the gravitational field and proves the importance of the gravitational model with high precision. Because the order of Delaunay variable G in the differential equation increases with the increasing terms of the gravitational field, which contributes significantly to the evolution of the orbit when G is sufficiently small, the influence of high-order terms gradually increases with increasing eccentricity.

Evolution period of traditional frozen orbits
To analyze the nonlinear dynamical property of the orbits near the equilibrium points, the canonical equation after the Hori-Lie transformation is linearized using the differential method.
Taylor expansion is performed around the equilibrium points to linearize the differential formulations of the mean orbital elements. The canonical equation is expanded via a Taylor expansion to approximately linearize the dynamical equation. For the sake of convenience, assuming point (G 0 , g 0 ) is the equilibrium point, the differential equation is expanded to a first-order Taylor expansion at (G 0 , g 0 ): The stability of the equilibrium point can be determined by the eigenvalues of the coefficient matrix A. The eigenvalues of the coefficient matrix A can be divided into two classes: one class is a pair of conjugate complex numbers for which the real part is zero, corresponding to the centers; the other class is a pair of real eigenvalues of opposite sign and equal magnitude, corresponding to the saddle points. Because the center is a stable point, the evolution of the orbit around a center is a closed loop; the saddle point is an unstable equilibrium point.
The center point corresponds to a frozen orbit, i.e., a periodic orbit, whereas the blue curves around it correspond to quasi-periodic orbits. The evolution periods of the quasi-periodic orbits around the frozen points (centers) can be obtained by integrating. The orbits around the selected center are indicated in Fig. 3. The points on the red dashed line from Point A to Point B in Fig. 3 are the initial values of the integrals, with g equal to π/2. The x -axis coordinates in Fig. 4 are the values of the Delaunay variables G in these initial values. The relationship between G and the periods is displayed in  As revealed in Fig. 4, the evolution periods of the orbits around the same center do not change with the initial G changing and can be considered the same value. Therefore, the evolution period of the frozen orbit can be defined by the eigenvalue of the equilibrium point. Assume [0 ± βi] are the eigenvalues of the frozen orbits. The evolution period of the frozen orbit can be given by T = 2π β based on the theory of ordinary differential equations that describes the evolution of the frozen orbits. This period, obtained by eigenvalue, represents the evolution periods of the orbits around the frozen orbit. Next, the eigenvalues of all the traditional equilibrium points in Fig. 1 are calculated to determine the stability of each equilibrium point and obtain the period of each center. Figure 2), eccentricity e is a monotonically decreasing function of the Delaunay variable G, i.e., period T increases as the eccentricity e of the center decreases.
An interesting phenomenon is indicated in Fig. 5. By defining the frozen orbits with similar values of eccentricity e as a family of frozen orbits in this chapter, it can be observed that the period of a family of frozen orbits changes slowly and actually remains constant as the Delaunay variable H changes. The period T approaches infinity during the derivative or disappearance of the center corresponding to the frozen orbit, confirming the result in Section 3.

Topology
and bifurcation of nontraditional equilibrium points

Evolution and topology of nontraditional equilibrium points
Apart from more traditional frozen orbits, a high-order gravitational field also brings the equilibrium points whose g values are neither 90 • nor 270 • , which are defined as nontraditional frozen orbits in Section 3. Because of the high-order terms of the gravitational field model containing the J 15 term, we calculate the accurate values of nontraditional equilibrium points and discover the global evolution of all equilibrium points containing centers and saddle points with the Delaunay variable H changing based on constant L. When the Delaunay variable g is not equal to 90 • or 270 • , the solution of Eq. (7) is transformed from a solution problem of a onedimensional equation into a two-dimensional equation (compared with that based on g equaling 90 • or 270 • ).
Because of the complexity of the differential equation, this solution requires particular methods to determine all the equilibrium points at any Delaunay variable H.
Variables e x and e y are used to replace the orbital elements e and ω used in several studies. The relationship between them is as follows: e x = e cos ω, e y = e sin ω. As indicated in Ref. [22], Coffey drew the contour map of F * in the space (e x , e y ) to obtain a local phase portrait within the small range of G, i.e., e. We obtain the same result under the same initial condition (the gravity field up to J 9 and J 15 in Fig. 6).
The method described above can be applied to draw a local phase portrait. However, the Hamiltonian function F * varies significantly with changes in G. We attempted to draw the contour map of F * to obtain the global phase portrait (G ∈ [H, L]) and obtained the following result as shown in Fig. 7.
It can be observed that the information for numerous equilibrium points is lost. For whatever the resolution of the counter map of F * selected, the total equilibrium points cannot be displayed entirely. This is because the Hamiltonian function F * changes dramatically with G in the high-order gravitational field. Thus, the contour map of F * can only present a local phase portrait; it cannot be used to draw the global phase portrait. The global phase portrait is drawn as described in the following.
First, select four typical values for Delaunay variable H based on the distribution of the traditional equilibrium   and Table 1 display the boundaries of each part and the corresponding number of equilibrium points, centers, and saddle points; they also display the evolution of the global frozen orbits for global H. Figure 8 indicates the respective phase portrait at the particular Delaunay variable H of each part that can describe the distribution of equilibrium points in the corresponding part. Figure 9 confirms the existence of nontraditional equilibrium points at any value of H. It indicates that the majority of the nontraditional equilibrium points are saddle points, and there are a maximum of four nontraditional centers for any Delaunay variable H value. As can be observed in Fig. 8, the number of traditional equilibrium points decreases as Delaunay variable H increases, accounting for the results in Section 3. In general, the total number of equilibrium points containing both traditional and nontraditional equilibrium points decreases as Delaunay variable H increases. However, the number of equilibrium points does not continue to decrease because of the bifurcation of nontraditional centers. The number of centers and saddle points are denoted as n c and n s , respectively. It can be concluded from Fig. 9 that n c = n s + 2. Thus, there always exists an even number of equilibrium points for any value of H, which can also be derived from the fact that n c = n s + 2. Furthermore, nontraditional equilibrium points are found dually.

Bifurcation of equilibrium points
The global evolution of frozen points with Delaunay variable H changing from 0 to 5.29 × 10 10 can also (1) (2)  Fig. 9 Corresponding phase portraits in each part of H from zero to 5.29 × 10 10 : red points indicate centers, green points indicate saddle points, and red stars indicate frozen orbits widely studied in other reports.
be obtained based on Fig. 9. Equilibrium points move, appear, and disappear as H changes, which results in a change in the number of equilibrium points. When H is a small value, the majority of the equilibrium points are located in the bottom region. The distribution of frozen points tends to be uniform as Delaunay variable H increases. Finally, when H is a large value, the majority of the equilibrium points disappear, and there are only two centers in Part (35). In addition, there is a maximum of 38 equilibrium points for the same Delaunay variable H in Parts (2), (5), and (7). Because centers and saddle points appear and disappear at the same time, the evolution of the centers is selected to describe the bifurcation in this study. There are two evolution cases for nontraditional centers. In the first case, either one traditional center changes into a pair of nontraditional centers and one traditional saddle point as H increases or an inverse process occurs. Figure 10 displays this evolution case, which reflects the evolution from Part (1) to Part (2), which is denoted by Parts (1)- (2). The evolutions between Parts (7)- (8) and (23)- (24) are also examples of this case. In the second case, a pair of nontraditional centers and a pair of nontraditional saddle points appear or disappear at the same time. Figure 11 displays this evolutionary case, which reflects the evolution between Parts (4)-(5), The evolution of traditional centers has a primary case. As indicated in Fig. 12, the two traditional centers and two saddle points constitute the basic structure. Thus, the stable range of one center is narrowed and the center and two saddle points approach until changing into one saddle point. Finally, the saddle point and other center approach until they disappear. This process causes two changes in the equilibrium points, where each change corresponds to a decrease in one center and one saddle point. This     3.734 × 10 9 -4.536 ×

Relationship between number of harmonic terms and frozen orbits
As mentioned in Section 4, the order of Eq. (7) increases with an increase in the terms of the gravitational field model. The zonal harmonic terms J n influence the number of equilibrium points. That is, both H and J n influence the distribution of equilibrium points at a fixed L. Figure 8 is redrawn in the gravity field up to J 16 . The number of equilibrium points under J 16 is four more than that in J 15 when H is less than 4.3 × 10 10 . Furthermore, the evolution boundaries under J 16 shift left marginally compared with those in J 15 . Next, the numerical relationship between the terms of the gravitational field and frozen points is studied. Phase portraits with different zonal harmonic terms of the gravitational field from J 3 to J 16 are drawn to facilitate the study. The relationship between the number of equilibrium points and terms of the gravitational field J n is displayed in Fig. 14. Figure 15 displays the corresponding phase portrait at H = 3 × 10 10 and L = 5.29096 × 10 10 for each term J n , which displays the corresponding distribution of frozen points. The terms of gravitational field J n indicate that the gravitational field model contains from J 2 to J n zonal harmonic terms.
As indicated in Figs. 14 and 15, the number of equilibrium points typically increases as the zonal harmonic terms of the gravitational field J n increase. Within the range of J n from J 3 to J 16 , the number of equilibrium points increases with increasing J n . This phenomenon is consistent with the hypothesis of Sections 3 and 4; the higher-order terms of the gravitational field result in a greater number of equilibrium points. Note that there are only two centers that are traditional and have no saddle points at the zonal harmonic terms of the gravitational field from J 3 to J 7 under the given Delaunay variables L and H, thus proving the limitations of the low-order gravitational field adopted in other reports. In addition, the number of centers is always two more than the number of saddle points with any order higher than J 3 .
As indicated in Fig. 14, the number of equilibrium points can continue to increase as J n increases (even higher than J 15 ). All the newly increasing equilibrium points are essentially large frozen eccentricities introduced by the J 16 term. More equilibrium points with large eccentricity in the higher-order gravity field can be expected to be discovered. However, the high-order terms have minimal effect on a frozen orbit with small eccentricity. The reference orbits in the majority of practical missions have small eccentricities. It is for this reason that a limited gravity field can complete these practical missions.

Conclusions
Frozen orbits have been widely used in numerous practical missions such as scientific measurements and observation because of their constant mean semi-major axis, semimajor axis, eccentricity, inclination, and the argument of perigee. The evolution and distribution of frozen orbits have not been investigated owing to the limited terms of the gravitational field model adopted in other studies. In this study, the gravitational field model was extended to J 15 zonal harmonic terms to observe these undiscovered frozen points and analyze the distribution and evolution of the frozen orbits. The Lagrangian planetary perturbation equation was transformed into a canonical form described in Delaunay variables, and the averaging process was achieved through the Hori-Lie transformation. Families of traditional frozen orbits whose argument of perigee was either 90 • or 270 • were extended into three families: a family of quasicircular orbits, family of large elliptic orbits, and family near critical inclination. The period of the traditional frozen orbits decreased with increasing eccentricity e and approached infinity rapidly when the corresponding frozen point appeared or disappeared. In addition, nontraditional frozen points were systematically reported in this paper. The distribution and evolution of global frozen orbits were investigated and the regularity of evolution was summarized. Finally, a comparison of phase portraits with different zonal harmonic terms from J 3 to J 16 indicated the limitations of the low-order gravitational field that have not been frequently reported. The results in the high-order gravitational field are the extension and complement of those in the low-order case.
This study focused on the distribution and evolution of global frozen orbits and investigated previously unreported properties. The results can be applied to any orbital height by adjusting L and H. More frozen points and frozen orbits provide more orbital options for designing missions with frozen orbits. In addition, the conclusions can be expanded to frozen orbital missions involving other oblate primary celestial bodies. (1) (2)  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.
The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.