Nonsingular vectorial reformulation of the short-period corrections in Kozai’s oblateness solution

We derive a new analytical solution for the first-order, short-periodic perturbations due to planetary oblateness and systematically compare our results to the classical Brouwer–Lyddane transformation. Our approach is based on the Milankovitch vectorial elements and is free of all the mathematical singularities. Being a non-canonical set, our derivation follows the scheme used by Kozai in his oblateness solution. We adopt the mean longitude as the fast variable and present a compact power-series solution in eccentricity for its short-periodic perturbations that relies on Hansen’s coefficients. We also use a numerical averaging algorithm based on the fast-Fourier transform to further validate our new mean-to-osculating and inverse transformations. This technique constitutes a new approach for deriving short-periodic corrections and exhibits performance that are comparable to other existing and well-established theories, with the advantage that it can be potentially extended to modeling non-conservative orbit perturbations.


Introduction
In the brief span of time after the launch of Sputnik, a whole succession of analyses was devoted to the problem poised by the drag-free motion of an artificial satellite about an fundamentally involves removing all terms that depend on the fast-varying mean anomaly, thus retaining the (secular and long-period) mean-element motion. As noted by Kozai (1959), the long-periodic perturbations of the first order come from the terms of the second order, and the singularity associated with the critical inclination only arises when such long-period effects are retained. Nevertheless, the distinction between secular and mean elements has often been muddled in the literature, and tabulated mean-to-osculating transformations (Gim and Alfriend 2003;Schaub and Junkins 2018) apply Brouwer's full periodic corrections. This approach will inevitably lead to significant errors near the critical inclination, but can be properly amended by neglecting the long-periodic terms (Breakwell and Vagners 1970).
An important application of the mean-to-osculating transformation is in the computation of the nominal osculating orbits from the frozen-orbit conditions determined in mean-elements space (Gurfil and Lara 2013). Frozen orbits correspond to equilibria for the averaged equations of motion, and, in the oblateness model, occur when secular effects due to even zonal harmonics are canceled by the long-period perturbations of the odd harmonics. While modern formulations compute frozen orbits directly from the non-averaged equations, based on the underlying quasi-periodic structure of librations around these mean equilibria, explicit analytical solutions form the starting point for the numerical optimization process. Recovering the short-periodic effects needed for this initialization can be troublesome at the critical inclination, whose precise frozen-orbit location is very sensitive to model truncation (Lara 2018).
Here, we present a new formulation of the mean-to-osculating and inverse conversions for first-order oblateness perturbations based on the Milankovitch elements Scheeres 2013, 2014). We use the direct method of Kozai (1959), as further elucidated by Scheeres (2012), and present an explicit analytical short-period correction in vector form that is valid for all elliptical orbits. We adopt the mean longitude as the fast variable and present a compact power-series solution in eccentricity for its short-periodic perturbations that can be truncated to achieve the necessary accuracy. We establish a numerical averaging approach based on the fast-Fourier transform (Uphoff 1973;Ely 2015) and make detailed comparisons between our vectorial solution and the classical Brouwer-Lyddane (BL) theory. For the latter, we adapt the more streamlined formulas presented in Schaub and Junkins (2018) and Gim and Alfriend (2003).

Analytical averaging
The basic idea in orbit-averaging methods is to obtain approximate equations for the system evolution that contain only slowly changing variables by exploiting the presence of a small dimensionless parameter that characterizes the size of the perturbation. The tacit assumption is that the perturbing forces are sufficiently weak so that these approximate mean equations of motion can be used to describe the secular and long-period orbital evolution. The perturbation equations in celestial mechanics, relating the time variation of the orbit parameters to the perturbing accelerations, in Gauss or Lagrange form, are nonlinear, nonautonomous, firstorder differential equations: (Alfriend et al. 2009;Scheeres 2012) in which g(x, t) is assumed to be T -periodic in time t. Equation (1) is trivially solved when = 0, yielding the integrals (Keplerian elements) in the unperturbed problem. The method of averaging consists of replacing Eq. (1) by the averaged autonomous system (Sanders et al. 2007 where the average is performed over time, and it is understood that x in the integrand is to be regarded as a constant during the averaging process. The basis for this approximation is the averaging principle, which states that in the general, non-resonant case, the short-period terms removed by averaging cause only small oscillations that are superimposed on the long-term solution described by the averaged system. Comparison between numerical integrations and the mean solution will in general show a divergence between the two as a result of an inconsistent choice of initial conditions, as depicted in Fig. 1. This offset can be understood by decomposing the osculating elements into mean and short-period components (Kozai 1959;Scheeres 2012): ( Differentiating, we can obtain an approximate equation governing the short-period dynamics: in which the mean state x has replaced the corresponding osculating elements x in the dynamical equations and where the small parameter has been omitted without lack of generality. From Eqs. (3) and (4), the short-periodic perturbations can be obtained as (Kozai 1959) x Accordingly, it can be seen that x sp = 0, which is also the implicit assumption in the averaging process. The interpretation of this result is that given an initial condition for a state x 0 = x(t 0 ), the mean equations have to be initialized at a value provided by Eqs. (3) and (5) as to have the averaged dynamics track the true evolution more closely. The removal of time, or analogously of the mean anomaly, requires computing the quadrature of functions depending implicitly on this variable through the true anomaly. The time averaging is performed over a periodic motion having a period much shorter than the time that characterizes the evolution of the dynamical system; this periodicity necessarily implies that averaging is taken for elliptical orbits. Given a quantity g(x, M) representing the right-hand side of the equations of motion, defined as a function of the dimensionless time variable M (the mean anomaly) in addition to the other orbital elements given by x, the average, Eq. (2), can be redefined as Fig. 1 Schematic showing a comparison between the numerical and mean solutions for an arbitrary orbital element, when starting from the same initial condition. The mean solution changes approximately linearly over one orbit and is referred to as the averaged variation. The short-period oscillations are the fluctuations that happen per orbit in the real motion "True" Solution (True) "Averaged" Solution Short-Period "Mean" Solution where the orbital elements x are held constant in the integration. Although the average is defined with respect to mean anomaly, it is often more easily calculated by means of the true or eccentric anomaly, f and E, respectively, using the differential relationships in which a and b = a √ 1 − e 2 are the semi-major and semi-minor axes, respectively, and e is the eccentricity, yielding the equivalent forms for averaging: Note that r can be expressed in terms of f and E as Here, H is the specific angular momentum and μ is the gravitational parameter. The Milankovitch elements consist of the two fundamental vectorial integrals of motion, namely, the eccentricity vector e and angular momentum vector H. These vectors can be parameterized in terms of the Keplerian elements relative to an inertial frame: where i is the inclination, Ω is the right ascension of the ascending node, and ω is the argument of periapsis. Because of the orthogonality constraint (i.e., H · e = 0), we need a sixth scalar element to fully define an orbit (Roy and Moran 1973). Adopting the mean longitude l = ω + Ω + M, the non-averaged equations of motion for an arbitrary disturbing acceleration a d can be stated in 'dyadic' form 2 as (Battin 1999;Rosengren and Scheeres 2014): where n 2 = μ/a 3 , p = H 2 /μ, andθ = ĥ ·r. Note that Eq. (13c) consists of terms that collect the contributions due to the three components in the radial, transverse, and normal directions of the disturbing acceleration and is valid for elliptical orbits in which e < 1. The position and velocity vectors, r and v, may be expressed as whereê = e/e,ê ⊥ = ĥ ·ê, andĥ = H/H . The Gauss equations have time-varying terms multiplying the accelerations involving the true anomaly, and thus they must each be averaged separately. The mean evolution of these elements can be computed asė The approximate short-period equations of motion for each element can then be formulated by subtracting Eq. (16) from Eq. (13), while holding all orbital elements but f constant. Kozai (1959), in his solution employing the classical elements a, e, i, Ω, ω, and M, obtained the non-averaged, mean, and short-period equations of motion using the Lagrange planetary equations. We note that a more general procedure has been outlined herein.
Particularly, for e and H, following Eq. (5), we have so that Special care, however, must be taken in the case of the mean longitude due to the presence of the mean motion appearing without any factor in Eq. (13c) (Kozai 1959). Expanding the osculating n into a first-order Taylor series about the mean elements, and following Eq. (4), we can write the approximate equation governing the short-period dynamics aṡ where Accordingly, and As a result, the short-periodic perturbations of l are given by

Numerical averaging
The conversion between mean and osculating elements can be obtained numerically, as previously shown by Walter (1967), Uphoff (1973), and Ely (2015). Here, we will exploit a numerical implementation of the near-identity transformation in Eq. (3) between averaged and osculating Milankovitch elements to validate our subsequent analytical developments. The mapping is obtained by numerically solving The boundary condition in Eq. (25) guarantees that the oscillations of x sp (t) are unbiased with respect to x. The formal solution is then (Sanders et al. 2007) where c k (x) are the Fourier coefficients of g (x, t), and i is the imaginary unit.
In this study, we numerically approximate Eq. (26) by truncating the series and by using the fast-Fourier transform (FFT) algorithm as discussed in Ely (2015), under the assumption that g is analytic on the continuation of nt in a non-vanishing complex strip. This assumption guarantees that the Fourier coefficients exhibit exponential decrease as a function of k. It results in a numerical averaging scheme that we call "FFT transformation" from now on. This scheme yields a first-order approximation of the motion, similarly to the Brouwer-Lyddane or Milankovitch ones. For the mean orbit longitude, we incorporate the additional terms in Eq. (19), resulting from the Taylor series of n, into Eq. (25), before performing the numerical quadrature.

Analytical short-period correction for oblateness perturbations
The quadrupolar (i.e., J 2 -truncated) disturbing function arising from an oblate planet can be stated in a general vector expression as (Scheeres 2012) where J 2 is the second zonal harmonic coefficient, R is the mean equatorial radius of the planet, andr = r/r from Eq. (14). Note that the Earth's spin axisp is assumed to be fixed in inertial space, and, as such, is aligned withẑ. The perturbing acceleration is then given by ∂R/∂ r as Accordingly, following Eq. (13), the perturbation equations can be stated aṡ The averaged equations of motion for the eccentricity and angular momentum vectors are given by Ward (1962) and Rosengren and Scheeres (2013). Averaging Eq. (29c) directly requires computing the quadrature of various dyadics and higher rank tensors of the dynamical variables, the details of which we omit. The mean equations can be stated aṡ where the bar operator is omitted from the elements because there is no ambiguity in what follows, i.e., all variables are averaged variables. Note that Eq. (30c) can be seen as the sum of the classical secular precession rates arising from planetary oblateness,l =Ṁ +ω +Ω, whereṀ Following the procedure outlined in Sect. 2.1 and detailed in Appendix B, the shortperiodic perturbations in the eccentricity vector, angular momentum vector, and mean longitude can be stated as where Roman numerals I 1 , I 2 , I I 11 , . . . designate various functions of true anomaly, I 2 , I I 12 , . . . represent the difference of these trigonometric expressions from their averaged values, I 1 , . . . are indefinite integrals of the previous core functions, and I 2 , . . ., represent differences between the doubly-integrated expressions, their averaged values, and the previous core averages. While somewhat cumbersome, the adopted notation mirrors the derivation of Appendix B and is otherwise systematic and methodical. The needed results pertaining to the solution, Eq. (32), are given by where all intermediate terms are given in Appendix B. They are obtained by also using the auxiliary formulas given in Appendix A. Thus, following Eq. (6), given the initial osculating state (e 0 , H 0 , l 0 ), the mean equations of motion, Eqs. (30), have to be initialized as

Validation of the Milankovitch scheme
In this section, we consider an extended grid of initial conditions and test the developed Milankovitch formulation against both BL and the fully numerical transformation of Sect. 2.2. For Brouwer-Lyddane, we have verified that the more streamlined formulas presented in Schaub and Junkins (2018) have been correctly transcribed according to their original sources, excepting the missing sin(2ω) factor in the long-period terms of M, ω, and Ω, which has only been noted in the recent erratum of the latest edition of this widely used monograph. 3 Furthermore, rather than using Lyddane's adhoc modification, Gim and Alfriend (2003) developed a new theory based on Brouwer's generating function that uses equinoctial elements. While still invalid at the critical inclination, Gim and Alfriend (2003) concluded from various numerical simulations that their method produces reasonable results within 0.25 • of this small divisor. Being inherently rooted in Brouwer's theory, it was expected at the outset that Lyddane (1963) and Gim and Alfriend (2003) would yield the same overall degree of accuracy. Nevertheless, for the sake of completeness, and because, superficially, it is not apparent that the formulas of Gim and Alfriend (2003) are mathematically equivalent to those systematized in Schaub and Junkins (2018), 4 we have also extended our numerical campaign to include a comparison between these as well (see Appendix C). Figure 2 shows a numerical confirmation of the validity of the Milankovitch formulation for satellites of the Sun-synchronous and Molniya type. The procedure was to first convert the initial osculating orbit into its corresponding mean elements using the developed formulas. These initial osculating and mean states were then propagated according to their dynamics described by Eqs. (29)  were converted into their corresponding mean elements using the developed Milankovitch formulation, and each set was propagated according to Eq. (29) (denoted "simulated osculating" dynamics in the legend) and Eq. (30) (denoted "simulated mean" dynamics), respectively. The osculating (cyan) and mean (yellow) trajectories were recovered (purple circles and red diamonds) at equal time steps using the aforementioned transformations, taking the respective simulated dynamics as input generally more accurate than, a simple Cowell integration in Cartesian space. The time histories of these evolutions at various subintervals were subsequently used as input to the respective osculating-to-mean and mean-to-osculating transformations in order to recover the aforementioned simulated trajectories. Projecting both the simulated and recovered osculating evolutions into the radial, alongtrack, and cross-track frame, we use the norm root mean square (RMS) of the positional error as a means to assess the accuracy of the transformation. To keep the presentation brief, we do not consider the velocity errors or other metrics. Figure 3 shows the results of this process for Brouwer-Lyddane, Milankovitch, and the FFT transformations, respectively, for Sun-synchronous-like orbits and nearly critically inclined, semi-synchronous ones. The dependence of the resulting errors on the choice of orbit orientation angles is also highlighted by Fig. 3.
For completeness and further validation, Fig. 4 compares the aforementioned Brouwer-Lyddane transformation which includes the long-period terms, with a BL mean-to-osculating implementation that omits them. These test cases are simulated in the trusted semi-analytical propagator, STELA, using a J 2 -only model. While slight discrepancies with STELA are to be expected due to different platforms, mean-element equations, and astronomical constants used in the respective simulations, the results are in good quantitative agreement. Importantly, the long-period terms in the full Brouwer theory do not cause appreciable changes in the recovered solutions throughout the whole of orbital phase space excepting a narrow band centered around the critical inclination. Figure 5 shows error maps corresponding to two different initial eccentricities for 500 × 500 grids of initial inclinations and semi-major axes. Each grid point, together with (M, ω, Ω) = (0, 0, 0), was used to form the full osculating element state vector and propagated for five orbital periods. The same procedure outlined in the previous paragraph was used to characterize the accuracy of each transformation, where the colorbar in Fig. 5 corresponds to the norm RMS of the difference between the recovered and simulated positions over the timescale of the propagations (five orbital periods). Note that prescribed limits were imposed   (i, a) values, for initial eccentricities of 0.01 (left) and 0.2 (right), and where the initial mean anomaly, perigee, and node angles were all set to zero. The colorbar represents the norm RMS of the difference between the recovered and simulated positions over five orbital periods, according to each formulation. The colorbar limit was set to the maximum error found between the Milankovitch and FFT schemes, as the Brouwer-Lyddane formulation becomes singular near the critical inclination. Grid points leading to unphysical errors or to values that exceed this limit are represented in white.  , and where the initial mean anomaly, perigee, and node angles were all set to zero. The colorbar represents the norm RMS of the difference between the recovered and simulated positions over five orbital periods, according to each formulation. The colorbar limit of each map was set to the maximum error found in the Milankovitch scheme to provide a better contrast. Grid points leading to unphysical errors or to values that exceed this limit are represented in white. on the colorbar of each map in this numerical campaign so as to more clearly highlight the differences among the various transformations. Figure 6 shows error maps in the semi-major axis-eccentricity plane for four different initial inclinations. For all three transformations, the error is highest at the boundary of allowable orbits, above which the perigee altitude would equal the Earth's radius.
The remaining slice of the action-like element space is given in Fig. 7, which shows how the errors manifest in the (i, e) plane for two values of initial semi-major axes (representing LEO satellites at 800 km altitude and GEO birds).
From this simulation campaign, we can conclude that the Milankovitch method consistently agrees with the classical Brouwer-Lyddane solution with each of the approaches  (i, a) values, for initial semi-major axes of R + 800 km (left) and a GEO (right), and where the initial mean anomaly, perigee, and node angles were all set to zero. The colorbar represents the norm RMS of the difference between the recovered and simulated positions over five orbital periods, according to each formulation. The colorbar limit of each map was set to the maximum error found in the Milankovitch scheme to provide a better contrast. Grid points leading to unphysical errors or to values that exceed this limit are represented in white. showing a slightly better accuracy in specific orbital regimes. Differences in positional error residuals remain in the order of a few kilometers or less.

Discussion
Our Milankovitch formulation is closed-form in both the eccentricity and angular momentum vectors. We have provided a general series solution for the mean longitude based on Hansen coefficients, which can be taken to any desired order of accuracy 5 , and shown that our scheme performs in agreement with standard BL theories. We note that our formulation bypasses the critical inclination because we do not consider the long-periodic terms that arise from a second-order perturbation treatment. On this account, we also omitted these in Brouwer's theory in Fig. 4 to emphasize that the long-periodic terms are negligible away from the critical inclination in the time-scale of interest.
Being a non-canonical set of elements, our derivation followed the approach used by Kozai (1959), as further elucidated by Scheeres (2012). We note, however, that, like Gim and Alfriend (2003), we could have merely adopted Brouwer's generating function and computed the short-periodic corrections using Poisson-bracket operations. Nevertheless, we chose to present an independent derivation that, as a positive outcome, results to be applicable to other perturbations for which a convenient generating function is not always available, see Shen et al. (2019). Having provided a general Kozai-like scheme, not rooted in canonical perturbation theory, our approach can be potentially extended to non-conservative forces, such as solar radiation pressure and atmospheric drag, in addition to treating lunisolar third-body gravity and other predominant perturbations. Future work will be devoted to the inclusion of long-periodic terms in our transformation and to the application of this method to modeling non-conservative perturbations.

Conclusions
We have developed a new mean-to-osculating and inverse transformation based on the Milankovitch elements with the mean longitude as the fast variable, which is valid for all eccentricity values smaller than one. An extensive numerical campaign was used to validate the vectorial transformation over orbital phase-space grids tailored to the relative distribution of cataloged Earth satellites and debris.
Acknowledgements A deep debt of gratitude is extended to Claudio Bombardelli for help in clarifying a subtle issue related to the mean motion in Kozai's solution, which has led to a substantial improvement in the Milankovitch formulation. Special thanks are also due to K. Terry Alfriend and Daniel J. Scheeres for insightful and motivating conversations, and to Carlos Yanez for running independent simulations using the reference CNES tool, STELA. We are also grateful to the anonymous referees for helping to improve the presentation.

Funding Open access funding provided by Alma Mater Studiorum -Università di Bologna within the CRUI-CARE Agreement.
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. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Auxiliary Formulas
We require the double quadrature of f , sin m f , cos m f , where m is an integer. The traditional approach has resorted to more or less elaborate expansions of various functions into sums of periodic terms, which depend chiefly on Taylor's series and Fourier's theorem, mainly because the integrals of these functions cannot be obtained conveniently in any other way (Brumberg 1995). While detailed explicit tabulations for these quantities may be found in Cayley's and Newcomb's classical tables, these results may be deduced from the Hansen coefficients (Hughes 1981). Namely, where the coefficients S n,m k and C n,m k for integers k, n, and m are power series in the eccentricity, starting with degree |m − k|, and are related to Hansen coefficients X n,m k by C n,m For k = 0, exact analytical expressions exist for the zero-order Hansen coefficients X n,m 0 for all values of n and m; in particular, we recall a special result for the zero-order Hansen coefficients, derived by Kozai (1962a): For k = 0, the analytical expressions for them do not terminate and, consequently, the series have to be truncated at some particular order in the eccentricity.
Hansen's coefficients can be expressed as series invoking Bessel and hypergeometric functions from which several formulae of recurrence can be derived that greatly facilitate their calculation (Challe and Laclaverie 1969;Giacaglia 1976): A useful series representation to seed the recursion, which dates back to Hansen, is given by (Proulx and McClain 1988) where the functions L i (x) are the generalized Laguerre polynomial defined as and In the preceding expressions, η = √ 1 − e 2 and β = (1 − η)/e. Also needed in our derivation is the classical equation of the center, which is given by where J k (x) are Bessel functions of the first kind.

B Derivation of the short-periodic perturbations in the Milankovitch elements
The procedure involves first writing down an approximate differential equation that governs the short-period dynamics for each variable (e, H, l), according to Eq. (4). The first indefinite time integration of Eq. (5) requires the quadrature of dyadics and higher rank tensors of the dynamical variable r. In particular, we can expose all of the fast-variable terms of Eq. (13) and can complete the needed core integrals by changing the dependent variable from time to true anomaly or mean anomaly: Specifically, we have For the vector term, we have Similarly, for the dyadic terms, we find where In the same way, we have for the first triadic terms: (1 + e cos f ) cos 3 fêêê + cos 2 f sin f (êêê ⊥ +êê ⊥ê +ê ⊥êê ) + cos f sin 2 f (êê ⊥ê⊥ +ê ⊥êê⊥ +ê ⊥ê⊥ê ) + sin 3 fê ⊥ê⊥ê⊥ d f = μ H 3 I I I 111êêê + I I I 112 (êêê ⊥ +êê ⊥ê +ê ⊥êê ) + I I I 122 (êê ⊥ê⊥ +ê ⊥êê⊥ +ê ⊥ê⊥ê ) + I I I 222ê⊥ê⊥ê⊥ .

Fig. 10
Error maps in the inclination-eccentricity plane using the Brouwer-Lyddane (top) and Gim-Alfriend (bottom) transformations, respectively. Each panel samples an equidistant grid of initial osculating (i, a) values, for initial semi-major axes of R + 800 km (left) and a GEO (right), and where the initial mean anomaly, perigee, and node angles were all set to zero. The colorbar represents the norm RMS of the difference between the recovered and simulated positions over five orbital periods, according to each formulation. The colorbar limit of each map was set to the maximum error found in the Milankovitch scheme (cf. Fig. 7) to provide a better contrast. Grid points leading to unphysical errors or to values that exceed this limit are represented in white. From left to right, the (maximum, mean) errors (in km) were (1. 2992, 0.2601) and (94.1968, 1.9108) for Brouwer-Lyddane, and (2.5030, 0.2355) and (80.6873, 1.8132) for Gim-Alfriend

Fig. 11
Error maps in the mean anomaly-periapsis angle plane (left) and periapsis-node angles plane (right) using the Brouwer-Lyddane (top) and Gim-Alfriend (bottom) transformations, respectively. Each panel samples an equidistant grid of initial osculating (M, ω) or (ω, Ω) values, for an initial semi-major axis of R + 800 km, eccentricity of 0.1, and inclination of 98 • . The remaining orbital element needed to form the full state vector was set to zero in each case. The colorbar represents the norm RMS of the difference between the recovered and simulated positions over five orbital periods, according to each formulation. The colorbar limit was set to the maximum error found across all schemes as there were no singular or outlier cases for these maps. From left to right, the (maximum, mean) errors (in km) were (1.3241, 0.4100) and (1.0332, 0.4704) for Brouwer-Lyddane, and (1.0173, 0.2060) and (1.3206, 0.4133) for Gim-Alfriend