On mean elements in artificial-satellite theory

The merits of a perturbation theory based on a mean-to-osculating transformation that is purely periodic in the fast angle are investigated. The exact separation of the perturbed Keplerian dynamics into purely short-period effects and long-period mean frequencies is achieved by a non-canonical transformation, which, therefore, cannot be obtained by Hamiltonian methods. For this case, the evolution of the mean elements strictly adheres to the average behavior of the osculating orbit. However, due to the unavoidable truncation of perturbation solutions, the fact that this kind of theory confines in the mean variations the long-period terms of the semimajor axis, how tiny they may be, can have adverse effects in the accuracy of long-term semi-analytic propagations based on it.


Introduction
Mean elements are useful in a variety of aerospace engineering tasks, like maneuver design or fast orbit propagation.They are roughly understood as the result of an averaging process whose aim is to isolate the long trend variations of the osculating orbit, which are modulated with short period fluctuations.While the averaging ensues from a particular transformation of variables it happens that different transformations can be chosen.Therefore, the "mean elements" terminology is loosely used with the meaning of elements that are free from short-period effects, a vague definition that may cause confusion in the implementation and use of available perturbation solutions in the literature -refer to the caveats in [1], pp.616, 690, and 695.
Ideally, the transformation from mean to osculating elements should be pure periodic in the mean anomaly, so that the mean elements comprise all the secular and long-period variations of the osculating orbit.But this is not often the case, and the mean to osculating transformation of different satellite theories may involve non-periodic and long-period terms in addition to the short-period terms.This is, in particular, the case of common analytical solutions in closed form of the eccentricity, in which some long-period terms remain hidden in the mean to osculating transformation due to the implicit dependence of the true anomaly on the mean anomaly and the eccentricity through Kepler's equation [2].Because of that, the mean orbit will depart from the average osculating orbit in long period displacements of small amplitude.
The lack of agreement between mean and average dynamics is not a concern in usual analytical or semi-analytical orbit theories, a case in which the osculating solution is obtained after recovering the periodic effects from the mean to osculating transformation.This is the case of usual solutions computed by canonical perturbation theory [3,4,5,6].However, the disagreement may be unwanted in those cases in which the information provided by the mean elements alone is crucial for the problem at hand, as may happen in some space geodesy applications [7,8] or when high accuracy is needed in preliminary orbit and maneuver design [9,10].Conversely, the reported aim of some non-canonical perturbation methods is to achieve the removal of only the short-period fluctuations from the mean orbit (see [11,12], and references therein).In this way, it is expected that the mean to osculating transformation is free from long-period terms, and, therefore, the resulting mean elements accurately represent the average behavior of the osculating orbit [8].
Taking these facts into account is of primordial importance when comparing canonical and non-canonical solutions [13,14].Moreover, special attention must be paid to the correct computation of the short-period elimination when Hamiltonian simplification techniques are involved in the computation of the canonical solution [15,16].On top on that some confusion is added since canonical perturbation methods have been reported to be particularizations of non-canonical methods for the appropriate choice of the arbitrary integration functions of the slow variables arising in both different methods [17].But this connection between canonical and non-canonical perturbation theories simply recalls the obvious: that canonical transformations are very particular instances of the wider field of transformation theory.Remarkably, the mean to osculating transformation that produces the exact separation between short-and long-period variations has been recognized as noncanonical for perturbed Keplerian motion [8].More precisely, the exact separation of terms that are pure periodic in the mean anomaly only can be guaranteed up to the first order of the small parameter when using Hamiltonian methods [18].
Regarding artificial satellite theory, this issue is mainly relevant in the computation of second order effects of the dominant zonal harmonic of the second degree, hereafter noted J 2 .Indeed, for Earth-like bodies second order effects of J 2 are comparable to those produced by the higher order harmonics of the gravitational potential, and, therefore, are routinely incorporated into operational software [19,5,20,21].
The generalized method of averaging (see [22] p. 168 and ff., for instance) is a suitable, non-canonical option when approaching semi-analytically the solution of non-trivial perturbation models, which may include both Hamiltonian and non-Hamiltonian perturbations [11].Rather, we focus on the simple J 2 perturbation and resort to the versatile method of Lie transforms to compute explicit expressions of the mean variations as well as their concomitant periodic corrections.Beyond its original aim of dealing with perturbed Hamiltonian flows [23], the Lie transforms method generally applies to systems of ordinary differential equations with minor modifications of the original formulation, and is readily implemented by means of efficient recursive algorithms [24,25].Most notably, the method allows for the explicit computation of the inverse, osculating to mean transformation, in addition to the direct, mean to osculating one.The former, we will see, is fundamental for understanding the consequences of choosing a mean to osculating transformation that is pure periodic in the mean anomaly in order to obtain the mean frequencies of the motion.
To clearly illustrate the non-canonical nature of a mean to osculating transformation that is pure periodic in the mean anomaly, we strip the formulation of non-essential complexities stemming from the closed form integration of implicit functions of the mean anomaly.That is, we base our exposition on a J 2 -type potential in which the implicit dependence of the J 2 -problem on the mean anomaly is avoided by resorting to the usual expansions of the elliptic motion [26].Furthermore, we truncate these expansions to the lower orders of the eccentricity to shorten the length of printed expressions.
The toy model is assumed exact in the computation of the mean variations so that the mandatory analytical tests carried out to ensure the correctness of the theory may provide immediate insight without need of making additional expansions.In particular, at each order of the perturbation approach we checked that the successive composition of the mean to osculating and the osculating to mean analytic transformations yields formally the identity up to the truncation order of the theory.We also checked that replacing the transformation equations of the short-period elimination in the original, osculating variation equations, leads to mean variation equations that match term for term the mean frequencies obtained with the perturbation approach up to the truncation order of the perturbation solution.These tests and model also serve us to illustrate the paradox that the inverse of a non-canonical transformation that is pure periodic in the fast angle may bear long-period as well as non-periodic terms.
Moreover, the toy model permits the computation of high orders of the per-turbation approach with relatively little computational effort.We need to compute them to show that, for a mean to osculating transformation that is pure periodic in the mean anomaly, the mean variation of the semimajor axis ceases to vanish at the third order.This is in clear contrast with the canonical transformation approach, in which making the mean anomaly ignorable in the Hamiltonian definitely results in constant mean semimajor axis at any order.Because of that, we may anticipate the inadequacy of a perturbation theory based on a pure periodic mean to osculating transformation for long-term propagation.Indeed, because the more important source of errors of a perturbation solution is related to the truncation of the mean frequencies, the neglected long-period effects in the variation of the semimajor axis unavoidably introduce long-period errors that will become dominant in the long term.Tests carried out to check the accuracy of perturbation solutions based on different choices of the mean to osculating transformation confirm this conjecture.
For the same reasons of simplicity and insight, we constructed the perturbation solution in Keplerian orbital elements, which are singular for circular orbits as well as for equatorial orbits.The implementation of alternative non-singular solutions useful for operational purposes can be approached analogously.
The paper is organized as follows.For completeness, the basic equations of the Lie transforms method are summarized in Section 2.Then, the toy model is formulated in usual orbital elements in Section 3. Also in this section we discuss the two basic instances in which the computation of the m-th order terms of the mean variations can be decoupled from the computation of the mean to osculating transformation of the same order, and illustrate the paradox that the inverse of a pure periodic transformation may yield long period effects of higher order than the first.The perturbation solution based on a mean to osculating transformation that is pure periodic in the mean anomaly is described in Section 4, showing that the appearance of long-period terms in the variation of the mean semimajor axis starts at the third order.Section 5 deals with the case of a perturbation solution based on the pure periodic character of the vectorial generator of the Lie transforms method, which, like in the Hamiltonian case, results in null mean variation of the semimajor axis at any order.Finally, selected propagations that illustrate the main features of each perturbation approach are presented in Section 6.

The Lie transforms method for vectorial flows
In this section we present the basic formulation that should allow interested readers to reproduce our computations.Full details in the derivation of the algorithms can be consulted in the original reference [24] as well as in the reformulation in [25] to which we rather adhere.
A Lie transformation x j " x j py k , q, with j and k integers ranging form 1 to l, is defined as the solution of the differential system where is the independent variable and the W j denote the l components of some vectorial generating function, for the initial conditions When such kind of transformation is applied to an analytical function F " F px j , q given by its Taylor series expansion the transformed function F ˚" F px j py k , q, q can be directly obtained as a Taylor series in the new variables by the recursive computation of the coefficients F 0,m from Deprit's triangle in which L denotes the linear, scalar operator and W k,m are the coefficients of the Taylor series expansion of the vectorial generator.Namely, On the other hand, the formulation of a differential system, say in the new variables y j " y j px k , q needs to take the Jacobian of the transformation into account.Namely, whose series expansion in the new variables is Now, the coefficients Φ j,0,n are obtained directly in the new variables from the new recursion Φ j,m,q`1 " Φ j,m`1,q `m ÿ iě0 ˆm i ˙Lj ,i`1 pΦ k,m´i,q q, j, k " 1, 2, . . ., l, (11) which still adheres to the structure of Deprit's triangle but with a different, vectorial operator under the summations.To wit, In a perturbation approach, we search for a generator that transforms the vectorial flow (8) in some desired, simplified form given by ( 10) -commonly with one of the variables removed up to a given order of the small parameter .For this task, the first order of recursion (11) yields where Φj,0,1 " Φ j,1,0 .After choosing the terms Φ j,0,1 in accordance with our simplification criterion, Eq. ( 13) is solved for the W j,1 from the system of partial differential equations obtained replacing L j,1 from Eq. (12).Analogously, at the second order, repeated iterations of Eq. ( 11) yield where Φj,0,2 " Φ j,2,0 `Lj ,1 pΦ k,1,0 q `Lj ,1 pΦ k,0,1 q depends only on terms of the first order, and hence is known.Then, after choosing Φ j,0,2 at our convenience, we solve Eq. ( 14) for W j,2 .
That is, in general, Eq. ( 11) is rearranged as the homological equation in which the terms Φj,0,m are known from previous computations, whereas those Φ j,0,m are chosen in agreement with some desired objective.Therefore, everything becomes known in the homological equation save for the m-components of the vectorial generating function (7), which are then solved from the corresponding partial differential system.The direct transformation is derived from the vectorial generating function W as a particular instance of the scalar recursion (5) for the functions x j " ř mě0 p m {m!q x j,m,0 px k q, in which x j,0,0 " x j and x j,m,0 " 0 for m ě 1.
The inverse transformation y j " y j px k , q is analogously obtained from the vectorial generating function V j " ´Wj , which in turn must be written in the y k variables.Because W is itself a vectorial flow, the reformulation of W " W px k py k , q, q as a power series in the y k variables is done based on Eq. ( 11), which is readily implemented replacing B j,i,0 " W j,i`1 px k q.Once the terms B j,0,i are computed recursively from Eq. ( 11), the inverse transformation is obtained by making V j,i`1 " V j,i`1 py k q " ´Bj,0,i px k q| "0 .In particular, V 1 " V 1 py k q " ´W1 px k q| "0 , from which it follows that the first order terms of direct and inverse transformations are formally opposite -yet they are evaluated in distinct sets of variables.On the other hand, making m " q " 0 in Eq. ( 11) immediately yields B j,0,1 " B j,1,0 " W j,2 , and hence V 2 " V 2 py k q " ´W2 px k q| "0 .
Remark that this trivial equivalence between direct and inverse generating function terms no longer applies beyond the second order of .

A J 2 -type perturbed Keplerian model
The J 2 nonspherical term exerts the dominant gravitatonial perturbation of the Keplerian motion on Earth orbiting artificial satellites.This is the reason why this perturbation model is known as the "main problem" in the theory of artificial satellites [3].For this problem, the only disturbance of the two body potential V " ´µ{r is due to the term where r " ap1 ´e2 q{p1 `e cos f q is the distance from the center of mass of the attracting body, a, e, I, Ω, ω, and M are usual Keplerian elements, denoting semimajor axis, eccentricity, inclination, right ascension of the ascending node, argument of the periapsis, and mean anomaly, respectively, whereas the true anomaly f is an implicit function of M and involves the solution of the Kepler equation.For a given body, the disturbing potential ( 18) is particularized by the physical parameters µ, R C , and J 2 , standing for the gravitational parameter, the equatorial radius, and the oblateness coefficient, respectively.The fact that the mean anomaly M , which is the fast angle to be removed in the computation of the mean variations, appears implicitly in Eq. ( 18) as a function of the true anomaly f , introduces additional complications in the solution of the integrals involved in the perturbation approach if the closed form is wanted.Therefore, in order to focus on the relevant facts of the solution in mean elements we rely on the usual expansions of the elliptic motion [26], which, besides, are truncated to the lower orders of the eccentricity.In this way we not only sidestep many of the difficulties arising in the computation of the closed form solution, but also shorten higher order expressions stemming from the perturbation solution to a manageable extent.
Therefore, we replace the disturbing potential of the main problem in Eq. ( 18), by the substitute that depends explicitly on the mean anomaly, in which we abbreviate s " sin I.
Then, the variations of the osculating orbital variables are cast in the form of Eq. ( 8), where " J 2 , and dx j {dt denote the time variations of a, e, I, Ω, ω, and M , for j " 1, . . .6, respectively.In particular, ´7" e 2 ps 2 ´2q `s2 ‰ cosp3M `2ωq where n " pµ{a 3 q 1{2 denotes the Keplerian mean motion, and we also abbreviate c " cos I. Terms Φ j,m,0 vanish for m ě 2. The appearance of the eccentricity in denominators of the first order terms of the variations of the argument of the periapsis and the mean anomaly given in Eqs. ( 26) and ( 27), respectively, may cause trouble in the integration of low eccentricity orbits.However, because the toy model is chosen just for illustrative purposes this is not a concern, and we simply will take care in our numerical experiments of choosing initial conditions far away enough from problematic cases.Certainly, a perturbation solution intended for operational purposes should rather be implemented in some set of non-singular variables.

Perturbed Keplerian motion. The homological equation
Remarkably, for perturbed Keplerian motion the homological equation of the Lie transforms method can be solved by indefinite integration.More precisely, for the m-th term, Eq. ( 15) turns into where Bn{Ba " ´3n{p2aq, as follows form the definition of the mean motion of a Keplerian flow.An attentive reader will have noticed the similarities between Eqs. ( 28)-( 29) and Eqs. ( 55)-(58) in [27].
It is important to note that Φ j,0,m should be chosen in such a way that it cancels the average of the known, tilde terms over the mean anomaly, namely, x Φj,0,m y M , in this way preventing the appearance of secular terms in the solution of the homological equation for j " 1, . . . 5. On the other hand, the choice of such Φ 6,0,m that cancels the terms x Φ6,0,m ´3 2 pn{aqW 1,m y M , thus avoiding the undesired secular terms, relies on the previous computation of W 1,m .More precisely, we only need to know xW 1,m y M " C 1,m pa, e, I, Ω, ω, ´q, which can be chosen in advance in two relevant cases due to the freedom provided by the arbitrary integration function that arises in the integration of Eq. ( 28).In particular, imposing W 1,m to be pure periodic in the mean anomaly makes C 1 " 0 trivially, and hence we can choose Φ 6,0,m " x Φ6,0,m y M like in the other cases.Alternatively, imposing the mean to osculating transformation of the semimajor axis to be pure periodic in the mean anomaly, allows for the direct computation of C 1,m from previous orders of the perturbation solution.
For the last case, the mean to osculating transformation in Eq. ( 16) of the semimajor axis a " a 1 `J2 a 0,1 `1 2 J 2 2 a 0,2 `. . . is recursively computed from Eq. ( 5) as follows.At the first order a 0,1 " L 1 paq " W 1,1 , as follows from the definition of the scalar operator in Eq. ( 6).Then, and hence Φ 6,0,1 " xΦ 6,1,0 y M .At the second order, a 0,2 " a 1,1 `L1 pa 0,1 q, where a 1,1 " W 1,2 , as follows from Eq. ( 5).Hence W 1,2 " a 0,2 ´L1 pa 0,1 q, where, for a pure periodic mean to osculating transformation xa 0,2 y M " 0. Therefore, the computation of C 1,2 " xW 1,2 y M yields that only relies on first order terms.Then, Φ 6,0,2 " x Φ6,0,2 y M ´3 2 pn{aqC 1,2 .Analogously, at the third order, Eq. ( 5) leads to the sequence from which W 1,3 " a 0,3 ´L1 pa 0,2 q ´L2 pa 0,1 q ´L1 pa 1,1 q.Therefore, the requirement that xa 0,3 y M " 0 yields and the long-period terms x Φ6,0,3 y M ´3 2 pn{aqC 1,3 that should be cancelled by Φ 6,0,3 are computed with the only knowledge of second order terms.And so on.In a paradox, the choice of arbitrary functions that make the mean to osculating transformation pure periodic may prevent the inverse transformation from having the same nature.Using vectors we denote x " x 1 ` x 0,1 px 1 q`1 2 2 x 0,2 px 1 q`Op 3 q the mean-to-osculating transformation, and x 1 " x ` x 1 0,1 pxq `1 2 2 x 1 0,2 pxq Òp 3 q its inverse.Neglecting terms Op 3 q and higher, their sequential composition yields x " where the term x 0,1 still needs to be expanded.That is, in which x 1 0,1 pxq and x 0,1 pxq cancel each other out because they are just opposite.Hence, where x 0,2 is pure periodic in M by our choice of the arbitrary integration functions.However, the product x 0,1 Bx 0,1 {Bx may give rise -and it certainly does for the J 2 problem-to constant and long-period terms in spite of each factor is pure periodic in the mean anomaly.Details on the construction of perturbation solutions based in these two noteworthy cases are provided in the two following Sections.

Theory 1
In this Section we obtain the perturbation solution based on a mean to osculating transformation that is pure periodic in the mean anomaly, which we label as "Theory 1".
At the first order n " 0 and q " 0 in Eq. ( 11), and hence Φj,0,1 " Φ j,1,0 .Besides, C 1,1 " 0, as follows from Eq. ( 30), and we choose Φ j,0,1 " x Φj,0,1 y M .We obtain Φ 1,0,1 " Φ 2,0,1 " Φ 3,0,1 " 0, and Then, the trivial integration of Eqs. ( 28)-( 29) for m " 1 provides the first order terms W j,1 of the vectorial generating function, each of which will depend on an arbitrary integration function C j,1 " C j,1 pa, e, I, Ω, ω, ´q save for C 1,1 , whose value has already been fixed.The mean to osculating transformation is of the form of Eq. ( 16), and is computed using recursion (5) in which F m,0 " 0 for m ě 1 and F is successively replaced by each orbital element x j P pa, e, I, Ω, ω, M q.In these equations, the arbitrariness of the functions C j,1 is suppressed by imposing that the transformation from mean to osculating elements be pure periodic in the mean anomaly, which yields C j,1 " 0.Then, we obtain ! 6re 2 p3s 2 ´4q `3s 2 ´2s sin M `3re 2 ps 2 ´2q `s2 s sinpM `2ωq ´6eps 2 ´1q ˆsinp2M `2ωq ´7re 2 ps 2 ´2q `s2 s sinp3M `2ωq The reader will forgive the abuse of notation with the aim of reducing the number of printed expressions, for the components of the generator W j,m depend on the osculating orbital variables whereas the periodic corrections are functions of the mean variables.This change of the osculating variables by the mean ones in the corrections, which also affects the mean variations, should be carried out at the end of the whole procedure in order to run the perturbation theory in a semianalytical way.On the other hand, the first order terms of the inverse, osculating to mean transformation, are formally opposite to those in Eqs. ( 37)-( 42), but they are evaluated in the original, osculating variables.At the second order, the solution of Eqs. ( 29)-( 28) needs the previous computation of the tilde terms Φj,0,2 .Repeated application of Eq. ( 11) yields whose evaluation only involves straightforward operations.Then, the second order terms of the mean variations of a, e, I, Ω, and ω, are chosen by averaging, to get 9 16 It follows the solution of Eqs. ( 28) and ( 29) for m " 2 to obtain the components of the vectorial generator W j,2 , which again introduce arbitrary integration functions that are particularized for obtaining a pure periodic mean to osculating transformation.In addition to the function C 1,2 in Eq. ( 49), for j ě 2 we readily obtain C j,2 " 0.
For instance, for the semimajor axis we obtain 51) where δ i,j denotes the Kronecker delta and the inclination polynomials P i,j,k are presented in Table 1.That a 0,2 is pure periodic in the mean anomaly results from the subindex j ě 1.On the contrary, the inverse, osculating to mean correction to the semimajor axis takes the form 3 16 (52) with the new inclination polynomials P 1 i,j,k of Table 2, and the long-period terms 3 16 with coefficients P 1 i,0,k in the same Table, which clearly show the non pure periodic nature of the osculating to mean transformation announced in Eq. ( 35).Indeed, it is immediate to see that xa 1 0,2 y M " a 1 0,2,long .The same features are indeed observed in the direct and inverse corrections to the other orbital variables.
This apparent inconsistency in the osculating to mean transformation does not corrupt the solution, which, as desired is made of the mean elements resulting from the numerical integration of the mean frequencies, and the analytic pure periodic corrections.Still, this kind of solution may be not the best one when the theory is intended for accurate long-term orbit propagation.To show that we compute the third-order terms of the mean variations.Recalling that Φ j,k,0 " 0 for k ě 2, from repeated iterations of recursion (11) we obtain Φ j,0,3 " L j,1 pΦ j,0,2 q `Φj,1,2 Φ j,1,2 " L j,2 pΦ j,0,1 q `Lj ,1 pΦ j,1,1 q `Φj,2,1 Φ j,2,1 " L j,3 pΦ j,0,0 q `2L j,2 pΦ j,1,0 q `Φj,3,0 ,  52) and (53).
Finally, it is worth mentioning that the formulation of Theory 1 in a different set of variables, as for instance, non-singular ones based on the semi-equinoctial elements κ " e cos ω, σ " e sin ω, while feasible, will loose the pure periodic character of the mean to osculating transformation as soon as we reach the second order.The reason is the same stated above that the product of two pure periodic terms yields, in general, non-periodic terms.Thus, for instance, up to second order terms κ " pe 1 ` e 0,1 `1 2 2 e 0,2 q cospω 1 ` ω 0,1 `1 2 2 ω 0,2 q, which, after expansion and truncation to Op 3 q yields κ " e 1 cos ω 1 ` `e0,1 cos where the products e 0,1 ω 0,1 and ω 2 0,1 will naturally spring long-period terms that will remain after reformulation of the orbital elements into semi-equinoctial variables.Therefore, Theory 1 should be recomputed in the new variables from scratch if we want to preserve the pure periodic character in the fast variable of the mean to osculating transformation.

Theory 2
Other case in which the successive order of the mean variations can be computed based only on previous orders of the vectorial generator is when we prescribe that the vectorial generating function be pure periodic in the mean anomaly.That is, the arbitrary integration functions C j,m " xW j,m y M vanish at any order of the perturbation theory.We label this case as Theory 2.
To the first order, both approaches match and the results in Section 4 apply also to this case.At the second order, for j ă 6 the mean frequencies Φ j,0,2 are the same as those given in Eqs. ( 44)-( 48), but the integration function C 1,2 in Eq. ( 49) must be replaced by C 1,2 " 0. Hence, 198η 3 p43s 4 ´24s 2 `8q `572η 2 s 2 ps 2 ´2q ´135ηp45s 4 ´24s 2 `8q ´1044s 2 ps 2 ´1q `6s 2 " 66η 3 p3s 2 ´2q `11η 2 p3s 2 ´4q whose structure is analogous to Eq. ( 50) but with different coefficients.The second order terms of the vectorial generator W j,2 are consequently the same as before for j ă 6, but change for W 6,2 .Since they all are involved in the computation of each transformation equation, now all of them are affected of long-period terms, which at this order are exactly the same in the direct and inverse transformations.
In particular, now a 0,2 " a 0,2,short `a0,2,long , where a 0,2,short is the same as Eq. ( 51), which only comprises short-period terms, and the long-period terms a 0,2,long are exactly one half of those reported in Eq. ( 53).
Analogously, now where the short-period terms a 1 0,2,short are those obtained from Eq. ( 52) as a 1 0,2 á1 0,2,long , whereas a 1 0,2,long " a 0,2,long .The same happens to the other variables, in which the long-period terms appearing in the osculating to mean transformations of the previous approach are now distributed in equal parts between the direct and inverse transformation.
Contrary to Theory 1, the current choice of a pure short-periodic vectorial generator makes that now the variation of the mean semimajor axis vanishes at each order.Besides, since there are not special requisites related to the periodic character of the short-period elimination, Theory 2 can be reformulated in a different set of variables without need of approaching the Lie transformations from scratch, as it also happens in the Hamiltonian case [28,29].

Semi-analytical propagations. A test case
The aim of this section is just to illustrate the different behavior between the two non-canonical perturbation theories described in the previous sections.The actual relevance of the differences observed between both kinds of solutions regarding the implementation of operational software falls out of the scope of the current study.Therefore, we are satisfied with running our tests for a favorable case in the orbital variables used, which we borrow from [19].In particular, we limit the reported tests to the case of an elliptical orbit with eccentricity e " 0.2, inclination I " 20 ˝, and semimajor axis a " 9500 km.The gravity field parameters used in the propagations are µ " 398600.4415km 3 {s 2 , R C " 6378.1363km, and J 2 " 0.001082634.
First of all, we recall that a semi-analytic propagation is made of three parts.The osculating to mean transformation is applied first to the initial conditions in order to get the initial state in mean variables, as well as, more importantly, to initialize the mean frequencies of the motion.Next, the mean variations are integrated numerically for the mean initial conditions, which is achieved with very long steps.Finally, osculating elements are obtained at each evaluation step from the mean-toosculating transformation.Because the more important source of errors stems from the mean elements propagation, usual perturbation theories include a higher order of the mean variations than the order of the periodic corrections [3,4,29,11].On the other hand, the initialization process is crucial to avoiding the abnormal growth of errors in the along-track direction, so the osculating to mean transformation is commonly patched including higher order information related to the semimajor axis conversion [30,31,32,33].In this fashion, we denote as m-th order a theory that comprises the pm `1q-th order effects of the mean frequencies and the m-th order effects of the transformation equations, but in which the osculating to mean transformation is patched with the pm `1q-th order effects of the semimajor axis.The usual calibration of the solution using the energy equation is not used in this case to strictly adhere to the general formulation of a non-Hamiltonian problem.
The first order effects of both distinct perturbation approaches yield identical terms.However, in our notation of a 1st-order theory they disagree in a 1 0,2 , cf.Eqs. ( 52)-(53) and Eq. ( 58), as well as in Φ 6,0,2 , cf.Eqs. ( 50) and (56).These small differences in the formulation produce small but significant changes in the time-history of the errors resulting from each perturbation theory, which are computed with respect to a reference orbit obtained from the numerical integration of the osculating variations provided by Eqs. ( 20)- (27).
The time-histories of the root sum square (RSS) of the position errors for a three-day semi-analytic propagation of the test case with each Theory are shown Fig. 1, where we observe that the RSS errors remain of the same order and present the same behavior.This kind of representation is useful in illustrating the accuracy of each theory for orbit prediction purposes, but does not illustrate relevant differences between the two distinct perturbation approaches.More precisely, Theory 1 yields mean elements that strictly adhere to the average behavior of the osculating orbit, as expected, which is not the case of Theory 2. This is illustrated in Fig. 2 for the semimajor axis, which shows that the propagation errors are of the same order and behave almost identically, with periodic oscillations of the same amplitude.However, the average of the errors, depicted by horizontal black lines in each plot of Fig. 2, nearly vanishes for Theory 1, averaging to about 1 cm in this example, whereas it is affected of a clear shift from the zero average in the case of Theory 2, of about 3 m in this particular propagation.Differences in the time-histories of the errors of the other orbital elements obtained with each theory, as well as in their respective averages, are of minor nature, and are not presented.In particular, with both theories the errors of the eccentricity average to " 10 ´6 for the test case, and to values below one tenth of arc second for the orbital angles.Taking additional terms clearly improves the accuracy of each perturbation theory, whose position errors reduce from tens of meters to several centimeters when propagating the test case with 2nd-order theories, in agreement with a OpJ 2 q improvement.This is shown in Fig. 3, where we also note the lower amplitude of the RSS errors resulting from Theory 1.Nevertheless, in view of the RSS errors remain in both cases with values of the expected order of the truncation of the perturbation theory, we consider that both theories are of comparable accuracy for the test case along the 3-day interval.In particular, in both propagations the errors are Op10 ´8q relative to distance.Correspondingly, the propagation errors of the semimajor axis reduce their amplitude from the meter to the centimeter level, as depicted in Fig. 4, also in agreement with the expected improvement provided by an additional order of the perturbation theories.As shown in the right plot of Fig. 4, the shift from the zeroaverage in Theory 2 remains, but now reduced to less than 1 cm due to the general improvements of the perturbation solution.The semimajor axis errors stemming from the semi-analytical propagation using Theory 1 fall now below the mm level.Regarding the errors of the other orbital elements, they remain small as well as similar in both cases, like it happened to the lower order solutions.On the other hand, as already pointed out in the text, the apparent superiority of Theory 1 must be better qualified.For the truncation of the non-vanishing mean variation of the semimajor axis, additional errors mainly affecting the mean motion, and hence the mean anomaly propagation turn this perturbation theory into a less accurate tool for long-term propagation.This is illustrated for the second order versions of the perturbation theories in Fig. 5, in which the propagation interval of the test case is extended from 3 days to 3 weeks.Now, rather than RSS errors, we display the projections of the position error in the intrinsic (along-track, radial, and cross-track) directions, which better illustrate the nature of the resulting errors.Thus, we clearly see that along-track errors resulting from Theory 1, notably deteriorate passed one week and a half, whereas those of Theory 2 only worsen slightly (top plots of Fig. 5).This is the expected consequence of the truncation error of the mean motion in the propagation of the mean anomaly, which does not affect Theory 2 as far as the mean variation of the semimajor axis vanishes at any order.Due to the elliptic character of the test orbit, this error in the mean anomaly propagation is also observed in the radial errors, yet with less adverse effects (center plots of Fig. 5), but barely affects the errors in the cross-track direction (bottom plots of Fig. 5).
Still, there is no doubt that Theory 1 is certainly correct.To check that, we patched both second order theories with the fourth-order terms of the mean variation of the semimajor axis, which changes nothing in Theory 2 for the vanishing of the mean semimajor axis in this approach, but clearly refines Theory 1.These additional corrections bring both theories again to comparable accuracy for orbit propagation purposes, as shown in Fig. 6, in which the plots in the right column are the same as corresponding ones in the right column of Fig. 5 but represented in different scales to ease comparison.It is worth to mention that in this last case, the 2nd-order patched Theory 1 provides errors of the semimajor axis that average to just a few hundredths of mm, whereas the average of the semimajor axis error of Theory 2 remains in the cm level in spite of the expanded propagation interval.As before, the errors of the other orbital elements presented by each distinct theory are analogous and very small.On the other hand, the amplitude of the errors remain similar to the previous case because the patched theories still rely only on 2nd-order direct corrections.Finally, we must mention that, beyond the illustration purposes of this Section, both perturbation solutions of the toy model have been also obtained in nonsingular variables -recomputed for Theory 1, to preserve the pure periodic character of the mean to osculating transformation, and simply reformulated in the case of Theory 2-in this way allowing us to carry out additional tests for less elliptical orbits that, while non-exhaustive, confirm the reported characteristics of each distinct perturbation approach.

Conclusions
For perturbed Keplerian motion, a non-canonical perturbation theory based on a mean to osculating transformation that is pure periodic in the fast angle has the advantage of providing mean elements that strictly adhere to the average evolution of the osculating orbit.On the other hand, extending this kind of perturbation theory to higher orders shows that the mean variation of the semimajor axis only vanishes up to second order of J 2 effects.Because the accuracy of a perturbation solution for a given time is essentially related to the truncation order of the mean frequencies, and on account of the Lyapunov instability which is inherent to Keplerian motion, this later fact turns the theory based on the pure periodic mean to osculating transformation less accurate for long-term propagation than classical perturbation approaches due to increasing errors in the along-track direction.Therefore, there is not an always-best perturbation approach, and the choice of the proper kind of perturbation theory must be done depending on the user's particular needs.
For didactic purposes the exposition avoids the use of implicit functions of the mean anomaly and relies on a toy model in common orbital elements, which has been derived from the J 2 problem by means of usual expansions of the elliptic motion.Computing an analogous theory in closed form requires to confront the solution of non-trivial integrals.Still, most of the solutions of these integrals are already reported in the literature or can be approached with known techniques of integral calculus.The computation of such kind of perturbation solution is under development, and will be published elsewhere when fully completed and tested.

Figure 2 :
Figure 2: Test case.Semimajor axis errors of the 1st order of Theory 1 (left) and 2 (right).

Figure 4 :
Figure 4: Test case.Semimajor axis errors of the 1st order of Theory 1 (left) and 2 (right).

Figure 5 :
Figure 5: Intrinsic errors of a 3-week semi-analytical propagation of the test case with 2nd-order Theory 1 (left) and 2 (right)

Figure 6 :
Figure 6: Intrinsic errors of a 3-week semi-analytical propagation of the test case with 2nd-order patched Theory 1 (left) and 2 (right)