Precise Analytical Computation of Frozen-Eccentricity , Low Earth Orbits in a Tesseral Potential

J 2 -J 3 modelmay not be accurate enough for some applications. Furthermore, lower order truncations of the geopotential are known to fail in describing the behavior of elliptic frozen orbits properly. Inclusion of a higher degree geopotential, which also takes into account the short-period effects of tesseral harmonics, allows for the precise computation of frozen-eccentricity, low Earth orbits that show smaller long-period effects in long-term propagations than those obtained when using the zonal model design.


Introduction
The procedure of designing and maintaining mission orbits for artificial satellites is made of different phases that go from simple to complex.First of all it is essential to have a deep knowledge of the natural dynamics involved, which requires determining the origin and magnitude of the forces that prevent exact Keplerian motion to happen.The main forces are then used to create a model that, while being simple enough for providing analytical insight, must retain the bulk of the real dynamics.The simplified model allows for the selection of the range of trajectories useful for the mission, from which the nominal trajectory will be chosen and refined in a more sophisticated model.Finally, propagations in a real, ephemeris model are required for maneuver planning.
The selection of the simplified model is eased by the common practice of separating perturbations into secular, long-period, and short-period effects [1].Then, long-term behavior is efficiently studied in mean elements, whose evolution equations are obtained by averaging the disturbing function over the fast rotating angles so as to remove shortperiod effects.Ideally, the long-term model would be solved analytically.However, as the analytical solution is rarely attained, one must be satisfied with taking advantage of the fair performances of semianalytical integration [2,3].
However, the perturbation problem may admit Hamiltonian formulation in some cases, basically when the perturbative forces are of gravitational origin.Then, in analogy to the Hamilton-Jacobi method, which allows integrating a problem by finding the canonical transformation to action-angle variables [4], one can resort to classical perturbation theories where canonical transformations are used to construct a set of formal integrals, such that the new Hamiltonian only depends on actions-the momenta conjugate to the angles-up to the required truncation order [5].
Traditionally, perturbation approaches that proceed by averaging all the angles are used to obtain an explicit analytical time solution [6,7].However, averaging is only allowed when the involved angles rotate or "circulate, " as is the case of the mean anomaly.Hence, a full averaging process that removes both short-and long-period effects to obtain the socalled "proper" elements [8] may miss important pendulartype solutions related to the oscillation or "libration" of some angles, thus constraining the validity of the analytical solution to certain regions in the phase space.On the contrary, the averaging may only be applied to the circulating angles in Mathematical Problems in Engineering order to reduce the Hamiltonian to one degree of freedom.Then, the Hamiltonian flow will be free of the constraints that apply to the full averaged problem.Even in the cases in which an explicit closed form solution to the reduced Hamiltonian is not known, the search for particular solutions of the reduced problem, and notably the stationary solutions, may provide important results.
A conspicuous example is the so-called Kozai mechanism [9], by means of which asteroids in high-inclination orbits may experience large orbital eccentricity and inclination variations under third-body perturbations.Kozai's eccentricityvarying orbits stem from an unstable stationary solution of a reduced Hamiltonian, involving only the variation of the eccentricity and the argument of the periapsis (see also [10]).The Kozai mechanism is not only of interest for astronomers [11,12], but also attracts the attention of aerospace engineers in view of its direct application to the case of science missions about natural satellites [13][14][15][16][17][18].
Similarly, orbits with stationary perigee and eccentricity, on average, are found to exist under the influence of the geopotential alone.These types of specialized orbits reduce stationkeeping requirements by using the natural dynamics.Furthermore, because almost constant eccentricity and fixedperigee orbits minimize altitude variations, they are found quite useful as nominal orbits for Earth mapping missions [19].In the usual aerospace engineering lingo, stationary solutions of these kinds are generally termed as frozen orbits and are loosely defined as orbits that, on average, maintain constant certain orbital design parameters, such as the inclination for sun-synchronous orbits, the node rotation rate as in the case of repeat-groundtrack orbits, or the eccentricity and argument of the perigee.The latter, which is the topic of this paper, are termed frozen-eccentricity orbits [20].Historical details on the frozen orbit concept can be found in [21].
In the preliminary design of Earth mapping missions, the frozen-eccentricity orbit is usually computed in a  2 - 3 model [19], although the influence of higher order harmonics may be required in some applications [22,23].Specifically, it must be noted that second-order effects of  2 and other zonal harmonics may introduce radical changes in the general frozen orbits picture when compared to the  2 - 3 case [21].
Because frozen orbits are determined in the mean elements space, whereas nominal orbits are realized by their osculating elements, the computation of precise initial conditions needs the readmission of the short-period effects that were removed in the averaging.This readmission is done by using the transformation from mean to osculating elements.In general, the frozen-eccentricity orbit design is efficiently carried out in a zonal model.However, in some cases the influence of tesserals may worsen the frozen-eccentricity condition considerably, thus making the effort of recovering short-period effects related only to zonals naive [24].In the present paper we extend the definition of frozen-eccentricity orbits to include the short-period effects due to tesserals in the computation of initial conditions for the nominal orbits, in addition to the long-and short-period effects due to zonal harmonics.
The paper is organized as follows.First of all, several classical facts on perturbations of the orbital elements due to noncentralities of the geopotential are recalled in Section 2.Then, in Section 3 we describe the main lines in the computation of the perturbation theory, which consists of the long-period Hamiltonian and the transformation equations from mean to osculating elements, by means of the Lie transforms method [25][26][27].The explicit time dependence introduced in the model by tesseral terms can be avoided by moving to a rotating frame.This change of the reference system requires the concomitant introduction of the Coriolis term into the Hamiltonian.In general, in order to make possible the elimination of short-period angles related to longitude-dependent terms from the Hamiltonian, a preliminary expansion of the true anomaly as a trigonometric series in the mean anomaly, whose coefficients are functions of the eccentricity, is required [28].These kinds of expansions constrain the use of a perturbation theory to the case of small or moderate eccentricities [6].However, in the case of close-Earth satellites, it is well known that the Coriolis term is of a higher order than the Keplerian [29].This fact allows the short-period effects of the geopotential affecting low Earth orbits (LEO), including those related to longitude-dependent terms, to be averaged in closed form of the eccentricity, with a notable reduction in the number of literal terms required by the perturbation theory [7,[30][31][32].Although this averaging is only valid away from tesseral resonances, deep tesseral resonances are not expected for the low-Earth orbits which are the topic of the present study.
The closed-form averaging is based on the computation of three canonical transformations of the Lie type.The standard elimination of the parallax [33] simplifies the original Hamiltonian by reducing the parallactic terms to the form  −2 , and by removing short-period terms related to the argument of the latitude.In the case of LEO orbits, a side effect of the standard elimination of the parallax is the removal of longitude-dependent terms, except for those related to the so-called -dailies [34].After this initial simplification, short-period terms depending on the mean anomaly are removed from the Hamiltonian via Delaunay normalization [35].It follows a simple averaging of the remaining terms that depend on the argument of the node in the rotating frame.As expected from the known fact that the effects of tesseral average out in the long term, except for tesseral resonances, we find that the reduced Hamiltonian is exactly the same as the zonal problem [21].
Numerical experiments related to the computation of frozen-eccentricity orbits are presented in Section 4. Because the Lie transforms are carried out only up to second-order effects of  2 , the computed frozen orbits are affected by residual long-period effects derived from the truncation order of the perturbation theory.In the case of elliptic orbits, we show that short-period terms due to tesserals may worsen the frozen-eccentricity solution to some extent, thus magnifying the residual long-period effects.On the contrary, for a variety of cases tested, the short periods introduced by the tesserals do not seem to play any relevant role for almost circular frozen-eccentricity orbits.Therefore, the classical frozen orbits theory based on a zonal model would be enough for computing accurate initial osculating elements of LEO, frozen-eccentricity, almost circular orbits.
Finally, it must be mentioned that we limit the model to a fifth-degree and order truncation of the geopotential, whereas a ninth-degree model has been suggested to minimize variations of model predictions when compared with long-term observations of real satellites [21].However, the fifth-degree and order model is simple enough to show all the relevant aspects of the frozen orbits picture and hence is considered sufficient for illustrating the theoretical facts discussed in the present research.On the contrary, higher order harmonics should be taken into account in the construction of a frozeneccentricity orbit theory intended for practical applications.

Noncentralities of the Geopotential: Classical Facts
We only focus on the effects of the geopotential and assume that the Earth is in uniform rotation about the -axis, without suffering from precessional or any other effects.Then, the Keplerian motion typical of mass-point attraction is disturbed by the noncentralities of the Earth's gravitational field, which in an Earth-centered, Earth-fixed coordinate system, are derived from the disturbing potential [ where  is the Earth's gravitational parameter,  is the equatorial radius of the Earth, (, , ) are spherical coordinates, with  standing for geocentric radius,  for latitude, and  for longitude,   = − ,0 ,  , , and  , are harmonic coefficients, and  , are associated Legendre polynomials of degree  and order .
For orbital applications, it is customary to formulate the disturbing potential (1) in orbital elements (, , , , Ω, and ), respectively, standing for semimajor axis, eccentricity, inclination, argument of the perigee, longitude of the ascending node, and mean anomaly.The conversion from spherical coordinates to orbital elements is straightforward, based on the transformation from spherical to Cartesian coordinates in the Earth-centered Earth-inertial frame: where (, , ) are usual Cartesian coordinates and  is the Greenwich sidereal time, and the fact that the orbital frame, defined by the position vector r and its derivative with respect to time ṙ , and the Earth-centered, Earth-inertial frame are related by simple rotations: where  =  +  is the argument of the latitude,  is the true anomaly,  = ‖r‖, and  1 and  3 , respectively, stand for standard rotation matrices about the and axes: ) . ( The conversion from spherical coordinates to orbital elements shapes the disturbing potential in (1) as a trigonometric series whose arguments are linear, integer combinations of the angles , , Ω, and , and whose coefficients depend on inclination and eccentricity.If, in addition, the true anomaly is expanded as a function of the mean anomaly, one gets [1] where The formulation of the disturbing potential in orbital elements is efficiently carried out using recursion formulas for the computation of the inclination and eccentricity functions, respectively,  and  in (5), as is explained in full detail, for instance, in the popular book by Kaula [1].Note that, in general, Kaula's  functions cannot be expressed in closed form of the eccentricity and, therefore, the use of ( 5) is limited to the case of small or moderate eccentricities.
In view of ( 6) and ( 7), the effects of the disturbing potential when expressed in orbital elements are generally classified in short-periodic effects, due to terms with  − 2 +  ̸ = 0 which are related to the period of the "fast angle" , long-periodic effects, due to terms with  − 2 +  =  = 0 related to the period of the perigee rotation, and secular effects produced by terms with  = 2 and  =  = 0. Specifically, even zonal harmonics contribute secular terms  ,,, = 0, as shown in (6), while odd zonal harmonics introduce long-period terms.
The period of terms depending on Ω −  is close to one day and may be comparable to the orbital period in thecase of high-altitude orbits.On the contrary, for LEO orbits their period is one order of magnitude higher than the orbital period and the inclusion of higher-order -daily terms in the disturbing potential is occasionally recommended [34,36].Terms with  − 2 +  = 0 and  ̸ = 0 are sometimes classified as medium-period terms.
The above classification of the disturbing function in secular, long-period, and short-period terms is quite useful in studying orbit evolution in the long term.The amplitude of the short-periodic effects is normally small and orbit evolution is conveniently studied by retaining only the secular and long-periodic terms of the disturbing function (5).Except for the case of tesseral resonances  Ṁ ≈  θ with  and  integers, where such combinations of the indices that ( − 2 + ) ≈  makes  ,,, evolve slowly, therefore contributing long-period effects to the disturbing function, the long-term disturbing function is represented by the zonal harmonics ( = 0) contribution to (5), which is further "averaged" over the mean anomaly ( − 2 +  = 0).
Lagrange planetary equations provide the rate of variations for the orbital elements, which require the computation of the partial derivatives of the disturbing function (5) with respect to the orbital elements themselves (see [5], for instance).It is easy to check that the partial derivatives of  with respect to , Ω, and -appearing in the Lagrange planetary equations for , , and -change sines into cosines in (7), and vice versa, thus preventing the appearance of secular terms in the semimajor axis, eccentricity, and inclination, as derived from (6).On the other hand, partial derivatives of  with respect to , , and -appearing in the Lagrange planetary equations for the mean anomaly, the longitude of the node, and the argument of the perigee-do not alter (7) and hence contribute secular terms to the rate of variation of , Ω, and , as derived from (6).When dealing only with long-period terms of the disturbing function, Lagrange planetary equations provide the evolution of a different set of elements, called mean orbital elements, which are sometimes written with primes to enhance the differences with respect to the instantaneous or osculating elements.
While the flow in mean elements is representative enough for long-term orbit evolution, the computation of initial conditions corresponding to particular solutions of the mean elements flow requires the knowledge of the transformation from osculating to mean elements: Besides, the long-term orbital evolution of some specific osculating elements can be efficiently explored by propagating corresponding mean elements, which in turn requires the knowledge of inverse transformation T −1 .As a closed form transformation is not available, one has to be satisfied with knowing some terms in the series expansion of the above transformation, which these days are customarily computed by the Lie transforms method [25][26][27].
The knowledge of the flow in mean elements alone may be enough in some engineering applications.But even in these cases, the simple removal of short-period terms  − 2 +  ̸ = 0 from the geopotential [22,23] misses important second-order effects of the second zonal harmonic, such that the averaging of short periods must be properly accomplished with the inclusion of second-order effects of  2,0 = − 2 .

Long-Term Disturbing Function
As the perturbation problem admits Hamiltonian formulation, we perform the averaging using perturbation theory by Lie transforms [25][26][27], which, in addition, provides the explicit transformation from mean to osculating elements, which is required in the computation of the nominal orbit.Besides, in view of  that always appears as a subtraction of the argument of the node, the explicit time dependence introduced by longitude-dependent terms, that is, tesseral and sectorial harmonics, is avoided by formulating the disturbing function in orbital elements in the Earth-centered Earthfixed frame, which is a rotating frame attached to the Earth.In consequence, the Hamiltonian must be supplemented with the Coriolis term.
Thus, the perturbation problem is materialized by the initial Hamiltonian: where  = √ 1 −  2 ,  =  1/2  −3/2 is the orbit mean motion, and  = cos , and orbital elements must be understood as functions of certain sets of canonical variables, such as polarnodal or Delaunay variables.Bear in mind that the polarnodal set (, , , , Θ, ) is formed by the radial distance , the argument of the latitude , the argument of the node , the radial velocity , the modulus of the angular momentum vector Θ, and the projection of the angular momentum vector along the Earth's rotation axis .The set of Delaunay variables (ℓ, , ℎ, , , ) is formed by the mean anomaly ℓ, the argument of the perigee , the argument of the node ℎ, the Delaunay action  = √ , the modulus of the angular momentum vector  = , and the component of the angular momentum vector along the Earth's rotation axis  = .Remember also that the dynamical relation 0 ≤ || ≤  ≤  holds.Besides, since we are using the rotating frame  = ℎ = Ω − .We limit our model to a fifth-degree and order truncation of the geopotential , which incorporates the fundamental qualitative dynamics of the real problem.This model is sufficient to illustrate the importance of taking shortperiod effects due to tesseral harmonics into account in the transformation equations from mean to osculating elements.The consideration of higher degrees and orders does not introduce relevant variations in our approach, except for the increase of the computational burden derived from the rapid growth in the size of the series used.
Given that we limit the theory only to the case of LEO orbits, deep tesseral resonances are not expected to happen and, therefore, the short-period effects to be removed by averaging are those related to the mean anomaly as well as to the argument of the node in the rotating frame.In addition, the ratio θ / ∼ √ 2 =  for LEO orbits, such that it is convenient to arrange the Hamiltonian in the following perturbative order: where is the Keplerian term, is the Coriolis term, and H 2,0 carries the contribution of  2,0 : where  = √ 1 −  2 = sin .Because the contribution of other harmonics of the geopotential is of the second order of  2,0 , which is O( 4 ), there are no terms of the order of  3 in the original Hamiltonian; hence, we choose where, instead of using ( 5), the disturbing potential  = (, , , −, Θ, ) is expressed in polar-nodal variables using (2) and (3), and taking into account that  = Ω −  and cos  = /Θ.Averaging is then made by Lie transforms up to the fourth order, based on the so-called "Lie triangle" inductive algorithm; a full account of the method can be consulted in the original reference [26].For the sake of making the analytical manipulation more efficient, the perturbation theory is split into three Lie transforms, each of which is of a different nature, as sketched out in the following.

Elimination of the Parallax.
First of all, we carry out a preparatory simplification (, , , , Θ, ) called the elimination of the parallax [33].Originally designed for the simplification of zonal models, this transformation only removes part of the short-period terms, those related to the argument of the latitude, while retaining the contribution of parallactic terms of the form 1/ 2 .This transformation was very soon extended to consider also tesseral effects, although the efficiency of the "extended parallax transformation" [37] may be questioned, due to the scarce simplification of tesseral terms and the large size of the generating function.However, in those cases in which the Coriolis term may be considered as a higher order than the Keplerian, the homological equation is released from  (28).
its dependence on the argument of the node; therefore, the standard elimination of the parallax can be applied, providing a notable simplification of tesseral terms [7,[30][31][32].
After applying the standard elimination of the parallax to the perturbed Hamiltonian (10), we find where, neglecting terms above the fourth order, we get and the eccentricity polynomials  , are Finally, the inclination polynomials  , ,  ,, , and  ,, are given in Tables 1, 2, and 3, respectively.It is worth noting that the orbital elements are now functions of the prime variables.(22).We abbreviate  ≡ 1 ± .The Lie transforms method also provides the generating function of the transformation T 1 , from which the explicit transformation equations  =  (  ,   ,   ,   , Θ  ,   ) ,  ∈ (, , , , Θ, ) (25) are obtained by means of the straightforward application of the Lie triangle (see [26], for details).

Delaunay Normalization.
First, the simplified Hamiltonian ( 16) is expressed in Delaunay elements.Then, the elimination of remaining short-period terms related to the mean anomaly is carried out by means of a new transformation: which, up to the truncation order, normalizes the Hamiltonian by introducing the Delaunay action as a formal integral of the perturbation problem [35].The Delaunay normalization is carried out in closed form of the eccentricity based on the standard differential relations between the true and mean anomaly  2 dℓ =  2 d, derived from the preservation of the angular momentum in the Kepler problem.
After the normalization, we obtain a new Hamiltonian: where, neglecting terms of orders higher than 4, we get D 1 = P 1 , D 2 = P 2 , D 3 = P 3 , and  (22).We abbreviate  ≡ 1 ± .where  , and  ,, are the same as before, and the inclination coefficients   , are given in the right column of Table 1.Note that now the orbital elements and the different functions of them are expressed in the double-prime Delaunay variables.
Again, the Lie transforms method provides the generating function of the transformation T 2 , from which the transformation T 2 is itself obtained by means of explicit equations.

Elimination of the Node.
The Hamiltonian is now free from short-period effects related to the mean anomaly; however, it still depends on the argument of the node in the rotating frame.Therefore, we carry out one more Lie transform (−,   , ℎ  ,   ,   ,   ) in order to remove ℎ  from the Hamiltonian, and thus introducing   as a new formal integral of the doubleaveraged problem.
The final Hamiltonian is where, up to the fourth order, K  = D  for  = 0, . . ., 3, and where  , ,  ,,0 , and the inclination polynomials   , are the same as before, and the orbital elements, as well as the different functions of them, are expressed now in the tripleprime Delaunay variables.
Except for tesseral resonances, the effects of tesseral perturbations average out in the long-term.Therefore, it is no surprise to find that the long term Hamiltonian ( 30) is exactly the same as the one provided in the appendix of [21] for a Mathematical Problems in Engineering zonal model, up to the truncation degree of the geopotential used in this paper.This agreement is an additional, important check on the validity of the transformations carried out in this paper.
The long-term flow is derived from Hamilton equations of the final Hamiltonian (30).Because of the averaging, ( 30) is of the form where   =   and   are constant.Therefore, the time evolution of the reduced phase space (  ,   ) decouples from the motion of ℓ  and ℎ  .Thus, we first solve the one degree of freedom system for the initial conditions   0 ,   0 to get Then, the evolution of ℓ  and ℎ  is obtained by the following quadratures: Although the analytical solution (34) may be computed in some cases, it would require the use of special functions, which do not provide the desired insight into the long-term evolution.On the other hand, particular solutions of the reduced flow may be quite useful in specific satellite missions.This is exactly the case of the stationary solutions of (33), corresponding to the frozen-eccentricity orbits which are the topic of this research.
Finally, we note that the reduced Hamiltonian (30) can be written in terms of the mean orbital elements as K = K(  ,   ;   ,   ) = K(, ; , ), where  = (  ,   ),  =   , and the two formal integrals of the reduced flow   and   have been replaced by  =  2 / and  = (, ).The latter is given by where   = arccos(  /  ) stands for the inclination of circular orbits, a case in which  = .

Frozen Orbits Computation: Numerical Experiments
Because the reduced Hamiltonian is of just one-degree-offreedom in the modulus of the angular momentum and the argument of the perigee, then for each pair of the formal integrals (  ,   ) the reduced flow is comprised of closed curves and equilibria.Note that the cylindrical map (  ,   ) introduces singularities in the phase space representation, which should be properly discussed in the invariants of the problem [38].For engineering applications, where rectilinear trajectories are not of practical interest, the reduced flow can also be studied in a sphere [39].Alternatively, insightful information can be obtained directly in orbital variables by means of global inclinationeccentricity diagrams (, ) and local eccentricity vector diagrams ( cos ,  sin ), which give a thorough picture of the solution space for a given semimajor axis  0 .Indeed, after the straightforward reformulation of the relevant equations in orbital elements and in view of the frozen-eccentricity orbits which appear in the "meridian"  = ±/2, (, )-diagrams of frozen-eccentricity orbits can be obtained by the simple, inexpensive evaluation of the contour level 0 = ġ (, ;  = ±/2;  =  0 ) of (33) in the variation range 0 ≤  < 1, 0 ≤  <  [40].
A sample inclination-eccentricity diagram of frozeneccentricity, direct inclination orbits is provided in Figure 1 for a mean semimajor axis  0 = 8000 km, where each point of the respective curves represents a frozen-eccentricity orbit for different values of the (averaged) energy and   (or ).An almost symmetric plot can be obtained analogously for the case of retrograde, frozen-eccentricity orbits.Note the agreement of Figure 1 with results in [21] except for quantitative variations due to the different truncation of the model and the different value selected for the mean semimajor axis.Thus, in Figure 1 we observe two different lines of frozeneccentricity orbits.The first starts from an equatorial almost circular orbit and is comprised of low-eccentricity orbits with the perigee fixed, on average, at 90 degrees, which exist for increasing inclination.Close to the critical inclination  = arccos√1/5 (due to the influence of higher order harmonics, the value of the critical inclination is slightly different from the value  = arccos√1/5, which corresponds to the  2 problem), the eccentricity of the frozen-eccentricity orbits increases with almost constant inclination.A new line starts from a frozen-eccentricity orbit close to the critical inclination, where  = 0.023,  = 63.45 deg., and  = 270 deg.(the large red point in the bottom plot of Figure 1).One branch of this line is comprised of frozen-eccentricity orbits whose eccentricity increases with almost constant inclination, whereas the other branch is comprised of loweccentricity orbits, which exist for increasing inclinations.
For the latter, the eccentricity of the orbits reduces with inclination until finding an exact circular orbit at  = 64.3533deg.Then, for increasing inclinations, eccentricity becomes greater and the argument of the perigee changes to 90 deg.We do not attempt to find other frozen-eccentricity orbits that could exist with very high eccentricities.
On the other hand, local eccentricity vector diagrams can be constructed for each pair of the formal integrals (  ,   ), or, equivalently, (, ), by means of simple contour plots of the reduced Hamiltonian (30).Eccentricity vector diagrams cannot show a clear representation of the reduced flow for orbits with high eccentricities.In spite of that, these kinds of diagrams succeed in presenting valuable information for usual aerospace engineering missions requiring low or moderate eccentricities, and hence their use is quite common in mission design and analysis procedures (see [42,43], for instance).Applications in this paper are limited to  < 0.2 and, therefore, corresponding eccentricity vector diagrams clearly show the frozen orbits' behavior.
Thus, Figure 2 shows a collection of eccentricity vector diagrams that illustrate the typical bifurcation sequence of frozen-eccentricity orbits in the vicinity of the critical inclination.The mean semimajor axis, a dynamical parameter of the reduced flow, is fixed at 8 000 km, corresponding to circular orbits' altitude of about 1 600 km over the Earth's surface-about 400 km below the upper limit for LEO.The variation of the other dynamical parameter of the reduced flow  produces the well-known sequence of bifurcations of frozen orbits in the vicinity of the critical inclination.Further details on the influence of different truncations of a zonal geopotential in the frozen-eccentricity orbits picture can be consulted in [21].While local plots are useful in illustrating qualitative dynamical aspects of the problem and for estimating mean values of frozen-eccentricity orbits of interest, the computation of nominal orbits requires precise knowledge of corresponding mean elements.This is done by computing stationary solutions of the Hamilton equations of the reduced long-term Hamiltonian K = K(, ; , ) in (33).Therefore, we must solve the system Stationary solutions are conveniently found by means of standard root-finding procedures that can be fed by approximate values of fixed points on eccentricity vector diagrams, which are directly estimated from the figures.Thus, for instance, in reference to the right-bottom plot of Figure 2 for the elliptic type solutions, corresponding to stable frozeneccentricity orbits, and for the hyperbolic-type solution, corresponding to an unstable frozen-eccentricity orbit.Even after the improvements, stationary solutions of (37) are just mean values.Therefore, direct propagation of these values in the nonaveraged problem will provide orbits that are only approximately frozen.This is illustrated in Figure 3 for the stable, elliptic frozen-eccentricity orbit, where we see the long-period oscillations of the argument of the perigee with an amplitude of about 30 deg and a period of about one hundred years.Correspondingly, the osculating eccentricity is affected by nonnegligible oscillations between the nominal value and  ≈ 0.128.
Classically, the nominal orbit is computed by recovering the short-period effects due to zonal harmonics, which are the only relevant perturbations in the long-term, assuming that there are no resonances.The readmission of zonal harmonics' short-period effects in the computation of osculating elements, as derived from the generating function of each Lie transform, indeed improves the frozen orbit condition, as shown in Figure 4, where we note that now the amplitude of the long-term oscillations in the eccentricity is almost negligible when compared with the amplitude of the inevitable short-period oscillations.The amplitude of the long-term oscillations of the perigee is also notably reduced to less than 2 deg.
As expected, the amplitude of the long-period oscillations can be further reduced when the computation of the osculating elements of the nominal orbit includes the readmission of short-period effects caused by tesseral harmonics.Thus, Figure 5 shows that now the amplitude of the oscillations in the argument of the perigee is further reduced to about one degree, while the long-period oscillations of the eccentricity are masked by the short periods, whereas the remaining longperiod effects should be attributed to the truncation order of the perturbation theory, as illustrated, for instance, in [44].
The improvements in the frozen orbit computation are further illustrated in the eccentricity vector diagram of  Figure 6, where we note that the readmission of shortperiod effects due to tesseral harmonics centers the osculating elements propagation over the predicted mean value.
The benefits of recovering short-periodic effects due to tesseral perturbations are not so evident in the case of the other two stationary solutions.Thus, as shown in the left plot of Figure 7, either taking or neglecting tesseral effects in the computation of the osculating orbit provides very good results, with slight improvements when recovering the tesseral effects (red plot) that provide an orbit better centered over the mean frozen value and with more regular shape.In the case of the unstable stationary solution (right plot  of Figure 7), the osculating orbit derails over the unstable manifold in both cases, such that the corrections added by the tesserals simply change the unstable branch that the orbit follows.

Conclusions
The classical engineering definition of frozen-eccentricity orbits as those orbits where secular effects due to even zonal harmonics are canceled by the long-period contribution of odd zonal harmonics, which has been more precisely defined in the past as stationary solutions of an averaged zonal model, is now amended to include the short-period effects from tesseral harmonics.Stationary solutions of the averaged model always carry residual long-period effects derived from the truncation order of the perturbation theory.But these undesired long-period effects can be minimized by recovering the short-period effects due to the tesseral harmonics in the computation of the osculating elements of the nominal frozen-eccentricity orbit.This theoretical result has been demonstrated for a five-degree and order truncation geopotential, which shares the main characteristics of the full model.To this effect, we have constructed a perturbation theory that considers up to second-order  2 effects.In the case of LEO orbits, the perturbation theory can be constructed in closed form of the eccentricity, thus notably reducing the number of terms required in the formal series, and, therefore, making their evaluation more efficient.The perturbation theory, of course, is not limited to the computation of frozeneccentricity orbits and can be applied to the fast and efficient semianalytical integration of LEO orbits.Thus, long-term evolution could be used to explain some observable drift in the eccentricity/argument of perigee of debris orbiting in special conditions.
We approached only the mathematical problem of the computation of frozen-eccentricity orbits under the noncentralities of the geopotential.The real, engineering problem is more complex, because second-order effects due to the nonfixed spin axis, lunisolar attraction, solar radiation pressure, or drag may destroy the frozen orbit condition, thus making it necessary to solve an optimization problem.Nevertheless, the proposed solution of this research will be an excellent starting solution in the optimization process.

Figure 1 :
Figure 1: (a) Inclination-eccentricity diagram of frozen-eccentricity orbits of the 5 × 5 tesseral problem, for a mean semimajor axis of 8000 km.(b) Detail in the vicinity of the critical inclination.Full lines are used for  = 90 deg., whereas dashed lines indicate  = 270 deg.

Figure 2 :
Figure 2: Bifurcation sequence of frozen-eccentricity orbits of the 5 × 5 tesseral problem, for a mean semimajor axis  = 8000 km.From left to right and from top to bottom, the inclination of the circular (nonfrozen) orbits is  = 63.33,63.43, 63.456, 63.46, 63.50, and 63.61 deg.

Figure 3 :
Figure 3: Long-term evolution of the stable, elliptic frozeneccentricity orbit in a 5 × 5 tesseral model when the osculating elements are taken the same as mean elements.

Figure 4 :
Figure 4: Long-term evolution of the stable, elliptic frozeneccentricity orbit in a 5×5 tesseral model when short-periodic effects from zonal harmonics are included in the computation of osculating elements.

Figure 5 :
Figure 5: Long-term evolution of the stable, elliptic frozeneccentricity orbit in a 5×5 tesseral model when short-periodic effects from zonal and tesseral harmonics are included in the computation of osculating elements.

Figure 6 :
Figure 6: Eccentricity vector evolution of the orbits in Figure 5 (red points) and Figure 4 (blue points).Contour levels are mean elements.

Figure 7 :
Figure 7: Eccentricity vector evolution of the frozen-eccentricity orbits when incorporating short-period effects from zonals (blue points) and also from tesserals (red points) superimposed to the eccentricity vector diagram in mean elements.Left: almost circular frozen-eccentricity orbit.Right: elliptic, unstable frozen-eccentricity orbit.