# On mean elements in artificial satellite theory††thanks: Submitted to
Celestial Mechanics and Dynamical Astronomy

Martin Lara

(February 13, 2023)

###### Abstract

The merits of a perturbation theory based on a mean to osculating
transformation that is pure periodic in the fast angle are investigated. The
exact separation of the purely short-period effects of the perturbed Keplerian
dynamics from the long-period mean frequencies is achieved by a non-canonical
transformation, which, therefore, cannot be computed 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 inescapable
truncation of perturbation solutions, the fact that this theory confines the
long-period variations of the semimajor axis into the mean variation
equations, how tiny they may be, can have adverse effects in the accuracy of
long-term semi-analytic propagations based on it.

##  1 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 non-canonical 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
perturbation 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.

##  2 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}(y_{k},\epsilon)$, with $j$ and $k$ integers
ranging form $1$ to $l$, is defined as the solution of the differential system

$\frac{\mathrm{d}{x}_{j}}{\mathrm{d}\epsilon}=W_{j}(x_{k},\epsilon),\qquad
j,k=1,2,\dots,l,$ (1)

where $\epsilon$ is the independent variable and the $W_{j}$ denote the $l$
components of some vectorial generating function, for the initial conditions

$y_{j}=x_{j}(y_{k},0).$ (2)

When such kind of transformation is applied to an analytical function
$F=F(x_{j},\epsilon)$ given by its Taylor series expansion

$F=\sum_{m\geq
0}\frac{\epsilon^{m}}{m!}\,\left.\frac{\partial^{m}F}{\partial\epsilon^{m}}\right|_{\epsilon=0}=\sum_{m\geq
0}\frac{\epsilon^{m}}{m!}\,F_{m,0}(x_{j}),$ (3)

the transformed function $F^{*}=F(x_{j}(y_{k},\epsilon),\epsilon)$ can be
directly obtained as a Taylor series in the new variables

$F^{*}=\sum_{m\geq
0}\frac{\epsilon^{m}}{m!}\,\left.\frac{\mathrm{d}^{m}F}{\mathrm{d}\epsilon^{m}}\right|_{\epsilon=0}=\sum_{m\geq
0}\frac{\epsilon^{m}}{m!}\,F_{0,m}(y_{k}),$ (4)

by the recursive computation of the coefficients $F_{0,m}$ from Deprit’s
triangle

$F_{m,q+1}=F_{m+1,q}+\sum_{i=0}^{m}{m\choose{i}}\mathcal{L}_{i+1}(F_{m-i,q}),$
(5)

in which $\mathcal{L}$ denotes the linear, scalar operator

$\mathcal{L}_{m}(\psi)=\sum_{k=1}^{l}\frac{\partial\psi}{\partial{x}_{k}}{W}_{k,m},$
(6)

and ${W}_{k,m}$ are the coefficients of the Taylor series expansion of the
vectorial generator. Namely,

$W_{j}=\sum_{m\geq 0}\frac{\epsilon^{m}}{m!}\,W_{j,m+1}(x_{k}),\qquad
j,k=1,2,\dots,l.$ (7)

On the other hand, the formulation of a differential system, say

$\frac{\mathrm{d}x_{j}}{\mathrm{d}t}=X_{j}(x_{k},\epsilon)\equiv\sum_{m\geq
0}\frac{\epsilon^{m}}{m!}\Phi_{j,m,0}(x_{k}),\quad\Phi_{j,m,0}=\left.\frac{\partial^{m}{X}_{j}}{\partial\epsilon^{m}}\right|_{\epsilon=0},$
(8)

in the new variables $y_{j}=y_{j}(x_{k},\epsilon)$ needs to take the Jacobian
of the transformation into account. Namely,

$\frac{\mathrm{d}y_{j}}{\mathrm{d}t}=\sum_{k=1}^{l}\frac{\partial{y}_{j}}{\partial{x}_{k}}X_{j}(x_{k},\epsilon)=Y_{j}(x_{k},\epsilon),$
(9)

whose series expansion in the new variables is

$\frac{\mathrm{d}y_{j}}{\mathrm{d}t}=Y_{j}(x_{k}(y_{i},\epsilon),\epsilon)\equiv\sum_{m\geq
0}\frac{\epsilon^{m}}{m!}\Phi_{j,0,m}(y_{k}),\quad\Phi_{j,0,m}=\left.\frac{\mathrm{d}^{m}{Y}_{j}}{\mathrm{d}\epsilon^{m}}\right|_{\epsilon=0}.$
(10)

Now, the coefficients $\Phi_{j,0,n}$ are obtained directly in the new
variables from the new recursion

$\Phi_{j,m,q+1}=\Phi_{j,m+1,q}+\sum_{i\geq
0}^{m}{m\choose{i}}\mathcal{L}^{*}_{j,i+1}(\Phi_{k,m-i,q}),\quad
j,k=1,2,\dots,l,$ (11)

which still adheres to the structure of Deprit’s triangle but with a
different, vectorial operator under the summations. To wit,

$\mathcal{L}^{*}_{j,i}(\Phi_{k})=\sum_{k=1}^{l}\left(\frac{\partial\Phi_{j}}{\partial{x}_{k}}W_{k,i}-\frac{\partial{W}_{j,i}}{\partial{x}_{k}}\Phi_{k}\right),\quad
j,k=1,2,\dots,l.$ (12)

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
$\epsilon$. For this task, the first order of recursion (11) yields

$\mathcal{L}^{*}_{j,1}(\Phi_{k,0,0})=\Phi_{j,0,1}-\tilde{\Phi}_{j,0,1},$ (13)

where $\tilde{\Phi}_{j,0,1}=\Phi_{j,1,0}$. After choosing the terms
$\Phi_{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 $\mathcal{L}^{*}_{j,1}$ from Eq. (12). Analogously, at the
second order, repeated iterations of Eq. (11) yield

$\mathcal{L}^{*}_{j,2}(\Phi_{k,0,0})=\Phi_{j,0,2}-\tilde{\Phi}_{j,0,2},$ (14)

where
$\tilde{\Phi}_{j,0,2}=\Phi_{j,2,0}+\mathcal{L}^{*}_{j,1}(\Phi_{k,1,0})+\mathcal{L}^{*}_{j,1}(\Phi_{k,0,1})$
depends only on terms of the first order, and hence is known. Then, after
choosing $\Phi_{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

$\mathcal{L}^{*}_{j,m}(\Phi_{k,0,0})=\Phi_{j,0,m}-\tilde{\Phi}_{j,0,m},\qquad
j=1,\dots,l,$ (15)

in which the terms $\tilde{\Phi}_{j,0,m}$ are known from previous
computations, whereas those $\Phi_{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

$x_{j}=x_{j}(y_{k},\epsilon)\equiv\sum_{m\geq
0}\frac{\epsilon^{m}}{m!}x_{j,0,m}(y_{k}),\qquad j,k=1,\dots,l,$ (16)

is derived from the vectorial generating function $W$ as a particular instance
of the scalar recursion (5) for the functions $x_{j}=\sum_{m\geq
0}(\epsilon^{m}/m!)\,x_{j,m,0}(x_{k})$, in which $x_{j,0,0}=x_{j}$ and
$x_{j,m,0}=0$ for $m\geq 1$.

The inverse transformation $y_{j}=y_{j}(x_{k},\epsilon)$ is analogously
obtained from the vectorial generating function $V_{j}=-W_{j}$, which in turn
must be written in the $y_{k}$ variables. Because $W$ is itself a vectorial
flow, the reformulation of $W=W(x_{k}(y_{k},\epsilon),\epsilon)$ 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}(x_{k})$. Once the terms $B_{j,0,i}$
are computed recursively from Eq. (11), the inverse transformation

$V_{j}=\sum_{m\geq 0}\frac{\epsilon^{m}}{m!}\,V_{j,m+1}(y_{k}),$ (17)

is obtained by making $V_{j,i+1}=V_{j,i+1}(y_{k})\equiv-
B_{j,0,i}(x_{k})|_{\epsilon=0}$. In particular, $V_{1}=V_{1}(y_{k})\equiv-
W_{1}(x_{k})|_{\epsilon=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}(y_{k})\equiv-W_{2}(x_{k})|_{\epsilon=0}$. Remark that this
trivial equivalence between direct and inverse generating function terms no
longer applies beyond the second order of $\epsilon$.

##  3 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 $\mathcal{V}=-\mu/r$ is due to the term

$\mathcal{D}=-\frac{\mu}{r}J_{2}\frac{R_{\Earth}^{2}}{r^{2}}\frac{1}{4}\left[2-3\sin^{2}I+3\sin^{2}I\cos
2(f+\omega)\right],$ (18)

where $r=a(1-e^{2})/(1+e\cos{f})$ is the distance from the center of mass of
the attracting body, $a$, $e$, $I$, $\Omega$, $\omega$, 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 $\mu$,
$R_{\Earth}$, 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

$\displaystyle\mathcal{T}=$
$\displaystyle\;\frac{\mu}{2a}J_{2}\frac{R_{\Earth}^{2}}{a^{2}}\frac{1}{4}\Big{[}2(3s^{2}-2)(1+3e\cos{M})+3es^{2}\cos(M+2\omega)$
$\displaystyle-6s^{2}\cos(2M+2\omega)-21es^{2}\cos(3M+2\omega)\Big{]},$ (19)

that depends explicitly on the mean anomaly, in which we abbreviate
$s\equiv\sin{I}$. Then, the variations of the osculating orbital variables are
cast in the form of Eq. (8), where $\epsilon=J_{2}$, and
$\mathrm{d}x_{j}/\mathrm{d}t$ denote the time variations of $a$, $e$, $I$,
$\Omega$, $\omega$, and $M$, for $j=1,\dots 6$, respectively. In particular,

$\displaystyle\Phi_{j,0,0}=$ $\displaystyle\;0,\qquad j<6,$ (20)
$\displaystyle\Phi_{6,0,0}=$ $\displaystyle\;n$ (21)
$\displaystyle\Phi_{1,1,0}=$
$\displaystyle\;an\frac{R_{\Earth}^{2}}{a^{2}}\frac{3}{4}\big{[}(6s^{2}-4)e\sin{M}+es^{2}\sin(M+2\omega)-4s^{2}\sin(2M+2\omega)$
$\displaystyle-21es^{2}\sin(3M+2\omega)\big{]},$ (22)
$\displaystyle\Phi_{2,1,0}=$
$\displaystyle\;n\frac{R_{\Earth}^{2}}{a^{2}}\frac{\eta}{e}\frac{3}{8}\big{[}e\eta(6s^{2}-4)\sin{M}+e(\eta-2)s^{2}\sin(M+2\omega)$
$\displaystyle-4(\eta-1)s^{2}\sin(2M+2\omega)-7e(3\eta-2)s^{2}\sin(3M+2\omega)\big{]},$
(23) $\displaystyle\Phi_{3,1,0}=$
$\displaystyle\;n\frac{R_{\Earth}^{2}}{a^{2}}\frac{cs}{\eta}\frac{3}{4}\big{[}e\sin(M+2\omega)-2\sin(2M+2\omega)-7e\sin(3M+2\omega)\big{]},$
(24) $\displaystyle\Phi_{4,1,0}=$
$\displaystyle-n\frac{R_{\Earth}^{2}}{a^{2}}\frac{c}{\eta}\frac{3}{4}\big{[}2+6e\cos{M}+e\cos(M+2\omega)-2\cos(2M+2\omega)$
$\displaystyle-7e\cos(3M+2\omega)\big{]},$ (25) $\displaystyle\Phi_{5,1,0}=$
$\displaystyle-n\frac{R_{\Earth}^{2}}{a^{2}}\frac{1}{e\eta}\frac{3}{8}\Big{\\{}4e(s^{2}-1)\big{[}1-\cos(2M+2\omega)\big{]}$
$\displaystyle+\big{[}e^{2}(6s^{2}-8)+6s^{2}-4\big{]}\cos{M}+\big{[}e^{2}(s^{2}-2)+s^{2}\big{]}\cos(M+2\omega)$
$\displaystyle-7\big{[}e^{2}(s^{2}-2)+s^{2}\big{]}\cos(3M+2\omega)\Big{\\}},$
(26) $\displaystyle\Phi_{6,1,0}=$
$\displaystyle-n\frac{R_{\Earth}^{2}}{a^{2}}\frac{1}{e}\frac{3}{8}\Big{\\{}4e\big{[}(3s^{2}-2-3s^{2}\cos(2M+2\omega)\big{]}+(7e^{2}-1)$
$\displaystyle\times\big{[}(6s^{2}-4)\cos{M}+s^{2}\cos(M+2\omega)-7s^{2}\cos(3M+2\omega)\big{]}\Big{\\}}.$
(27)

where $n=(\mu/a^{3})^{1/2}$ denotes the Keplerian mean motion, and we also
abbreviate $c\equiv\cos{I}$. Terms $\Phi_{j,m,0}$ vanish for $m\geq 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.

###  3.1 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

$\displaystyle W_{j,m}=$
$\displaystyle\;\frac{1}{n}\int\big{(}\tilde{\Phi}_{j,0,m}-\Phi_{j,0,m}\big{)}\mathrm{d}M,\qquad
j=1,\dots 5,$ (28) $\displaystyle W_{6,m}=$
$\displaystyle\;\frac{1}{n}\int\Big{(}\tilde{\Phi}_{6,0,m}-\Phi_{6,0,m}+\frac{\partial
n}{\partial{a}}W_{1,m}\Big{)}\mathrm{d}M,$ (29)

where $\partial{n}/\partial{a}=-3n/(2a)$, 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 $\Phi_{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, $\langle\tilde{\Phi}_{j,0,m}\rangle_{M}$, in this way preventing the
appearance of secular terms in the solution of the homological equation for
$j=1,\dots 5$. On the other hand, the choice of such $\Phi_{6,0,m}$ that
cancels the terms
$\langle\tilde{\Phi}_{6,0,m}-\frac{3}{2}(n/a)W_{1,m}\rangle_{M}$, thus
avoiding the undesired secular terms, relies on the previous computation of
$W_{1,m}$. More precisely, we only need to know
$\langle{W}_{1,m}\rangle_{M}=C_{1,m}(a,e,I,\Omega,\omega,-)$, 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
$\Phi_{6,0,m}=\langle\tilde{\Phi}_{6,0,m}\rangle_{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^{\prime}+J_{2}a_{0,1}+\frac{1}{2}J_{2}^{2}a_{0,2}+\dots$
is recursively computed from Eq. (5) as follows. At the first order
$a_{0,1}=\mathcal{L}_{1}(a)=W_{1,1}$, as follows from the definition of the
scalar operator in Eq. (6). Then,

$C_{1,1}=0,$ (30)

and hence $\Phi_{6,0,1}=\langle\Phi_{6,1,0}\rangle_{M}$. At the second order,
$a_{0,2}=a_{1,1}+\mathcal{L}_{1}(a_{0,1})$, where $a_{1,1}=W_{1,2}$, as
follows from Eq. (5). Hence $W_{1,2}=a_{0,2}-\mathcal{L}_{1}(a_{0,1})$, where,
for a pure periodic mean to osculating transformation
$\langle{a}_{0,2}\rangle_{M}=0$. Therefore, the computation of
$C_{1,2}=\langle{W}_{1,2}\rangle_{M}$ yields

$C_{1,2}=-\langle\mathcal{L}_{1}(a_{0,1})\rangle_{M}=-\sum_{k=1}^{6}\Big{\langle}\frac{\partial{a}_{0,1}}{\partial{x}_{k}}W_{k,1}\Big{\rangle}_{M},$
(31)

that only relies on first order terms. Then,
$\Phi_{6,0,2}=\langle\tilde{\Phi}_{6,0,2}\rangle_{M}-\frac{3}{2}(n/a)C_{1,2}$.
Analogously, at the third order, Eq. (5) leads to the sequence

$\begin{array}[]{rl}a_{2,1}=&W_{1,3}\\\
a_{1,2}=&\mathcal{L}_{2}(a_{0,1})+\mathcal{L}_{1}(a_{1,1})+a_{2,1}\\\
a_{0,3}=&\mathcal{L}_{1}(a_{0,2})+a_{1,2}\end{array}$

from which
$W_{1,3}=a_{0,3}-\mathcal{L}_{1}(a_{0,2})-\mathcal{L}_{2}(a_{0,1})-\mathcal{L}_{1}(a_{1,1})$.
Therefore, the requirement that $\langle{a}_{0,3}\rangle_{M}=0$ yields

$C_{1,3}=-\langle\mathcal{L}_{1}(a_{0,2})\rangle_{M}-\langle\mathcal{L}_{2}(a_{0,1})\rangle_{M}-\langle\mathcal{L}_{1}(a_{1,1})\rangle_{M},$
(32)

and the long-period terms
$\langle\tilde{\Phi}_{6,0,3}\rangle_{M}-\frac{3}{2}(n/a)C_{1,3}$ that should
be cancelled by $\Phi_{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
$\mbox{\boldmath$x$}=\mbox{\boldmath$x$}^{\prime}+\epsilon\mbox{\boldmath$x$}_{0,1}(\mbox{\boldmath$x$}^{\prime})+\frac{1}{2}\epsilon^{2}\mbox{\boldmath$x$}_{0,2}(\mbox{\boldmath$x$}^{\prime})+\mathcal{O}(\epsilon^{3})$
the mean-to-osculating transformation, and
$\mbox{\boldmath$x$}^{\prime}=\mbox{\boldmath$x$}+\epsilon\mbox{\boldmath$x$}^{\prime}_{0,1}(\mbox{\boldmath$x$})+\frac{1}{2}\epsilon^{2}\mbox{\boldmath$x$}_{0,2}^{\prime}(\mbox{\boldmath$x$})+\mathcal{O}(\epsilon^{3})$
its inverse. Neglecting terms $\mathcal{O}(\epsilon^{3})$ and higher, their
sequential composition yields

$\mbox{\boldmath$x$}=\Big{[}\mbox{\boldmath$x$}+\epsilon\mbox{\boldmath$x$}^{\prime}_{0,1}(\mbox{\boldmath$x$})+\frac{1}{2}\epsilon^{2}\mbox{\boldmath$x$}_{0,2}^{\prime}(\mbox{\boldmath$x$})\Big{]}+\epsilon\mbox{\boldmath$x$}_{0,1}(\mbox{\boldmath$x$}+\epsilon\mbox{\boldmath$x$}^{\prime}_{0,1}(\mbox{\boldmath$x$}))+\frac{1}{2}\epsilon^{2}\mbox{\boldmath$x$}_{0,2}(\mbox{\boldmath$x$}),$
(33)

where the term $\epsilon\mbox{\boldmath$x$}_{0,1}$ still needs to be expanded.
That is,

$\mbox{\boldmath$0$}=\epsilon\left[\mbox{\boldmath$x$}^{\prime}_{0,1}(\mbox{\boldmath$x$})+\mbox{\boldmath$x$}_{0,1}(\mbox{\boldmath$x$})\right]+\frac{1}{2}\epsilon^{2}\left[\mbox{\boldmath$x$}_{0,2}^{\prime}(\mbox{\boldmath$x$})+2\mbox{\boldmath$x$}^{\prime}_{0,1}(\mbox{\boldmath$x$})\frac{\partial\mbox{\boldmath$x$}_{0,1}(\mbox{\boldmath$x$})}{\partial\mbox{\boldmath$x$}}+\mbox{\boldmath$x$}_{0,2}(\mbox{\boldmath$x$})\right],$
(34)

in which $\mbox{\boldmath$x$}^{\prime}_{0,1}(\mbox{\boldmath$x$})$ and
$\mbox{\boldmath$x$}_{0,1}(\mbox{\boldmath$x$})$ cancel each other out because
they are just opposite. Hence,

$\mbox{\boldmath$x$}_{0,2}^{\prime}=-2\mbox{\boldmath$x$}_{0,1}(\mbox{\boldmath$x$})\frac{\partial\mbox{\boldmath$x$}_{0,1}(\mbox{\boldmath$x$})}{\partial\mbox{\boldmath$x$}}-\mbox{\boldmath$x$}_{0,2}(\mbox{\boldmath$x$}),$
(35)

where $\mbox{\boldmath$x$}_{0,2}$ is pure periodic in $M$ by our choice of the
arbitrary integration functions. However, the product
$\mbox{\boldmath$x$}_{0,1}\partial\mbox{\boldmath$x$}_{0,1}/\partial\mbox{\boldmath$x$}$
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.

##  4 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
$\tilde{\Phi}_{j,0,1}=\Phi_{j,1,0}$. Besides, $C_{1,1}=0$, as follows from Eq.
(30), and we choose $\Phi_{j,0,1}=\langle\tilde{\Phi}_{j,0,1}\rangle_{M}$. We
obtain $\Phi_{1,0,1}=\Phi_{2,0,1}=\Phi_{3,0,1}=0$, and

$\Phi_{4,0,1}=-n\frac{R_{\Earth}^{2}}{a^{2}}\frac{3c}{2\eta},\quad\Phi_{5,0,1}=-c\Phi_{4,0,1},\quad\Phi_{6,0,1}=n\frac{R_{\Earth}^{2}}{a^{2}}\frac{3}{2}(3c^{2}-1).$
(36)

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}\equiv{C}_{j,1}(a,e,I,\Omega,\omega,-)$ 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\geq 1$ and $F$ is
successively replaced by each orbital element
$x_{j}\in(a,e,I,\Omega,\omega,M)$. 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

$\displaystyle W_{1,1}=a_{0,1}=$
$\displaystyle-a\frac{R_{\Earth}^{2}}{a^{2}}\frac{3}{4}\big{[}e(6s^{2}-4)\cos{M}+es^{2}\cos(M+2\omega)-2s^{2}$
$\displaystyle\times\cos(2M+2\omega)-7es^{2}\cos(3M+2\omega)\big{]},$ (37)
$\displaystyle W_{2,1}=e_{0,1}=$
$\displaystyle-\frac{R_{\Earth}^{2}}{a^{2}}\frac{1}{8}\Big{[}6\eta^{2}(3s^{2}-2)\cos{M}+3\eta(\eta-2)s^{2}\cos(M+2\omega)$
$\displaystyle+\frac{6e\eta
s^{2}}{1+\eta}\cos(2M+2\omega)-7\eta(3\eta-2)s^{2}\cos(3M+2\omega)\Big{]},$
(38) $\displaystyle W_{3,1}=I_{0,1}=$
$\displaystyle-\frac{R_{\Earth}^{2}}{a^{2}}\frac{cs}{4\eta}\big{[}3e\cos(M+2\omega)-3\cos(2M+2\omega)$
$\displaystyle-7e\cos(3M+2\omega)\big{]},$ (39) $\displaystyle
W_{4,1}=\Omega_{0,1}=$
$\displaystyle-\frac{R_{\Earth}^{2}}{a^{2}}\frac{c}{4\eta}\big{[}18e\sin{M}+3e\sin(M+2\omega)-3\sin(2M+2\omega)$
$\displaystyle-7e\sin(3M+2\omega)\big{]},$ (40) $\displaystyle
W_{5,1}=\omega_{0,1}=$
$\displaystyle-\frac{R_{\Earth}^{2}}{a^{2}}\frac{1}{e\eta}\frac{1}{8}\Big{\\{}6[e^{2}(3s^{2}-4)+3s^{2}-2]\sin{M}$
$\displaystyle+3[e^{2}(s^{2}-2)+s^{2}]\sin(M+2\omega)-6e(s^{2}-1)$
$\displaystyle\times\sin(2M+2\omega)-7[e^{2}(s^{2}-2)+s^{2}]\sin(3M+2\omega)\Big{\\}}$
(41) $\displaystyle W_{6,1}=M_{0,1}=$
$\displaystyle-\frac{R_{\Earth}^{2}}{a^{2}}\frac{1}{8e}\Big{\\{}(4e^{2}-1)\big{[}6(3s^{2}-2)\sin{M}+3s^{2}\sin(M+2\omega)$
$\displaystyle-7s^{2}\sin(3M+2\omega)\big{]}-9es^{2}\sin(2M+2\omega)\Big{\\}}.$
(42)

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 semi-analytical 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 $\tilde{\Phi}_{j,0,2}$. Repeated application of
Eq. (11) yields

$\tilde{\Phi}_{j,0,2}=\Phi_{j,2,0}+\sum_{k=1}^{6}\left[\frac{\partial(\Phi_{j,0,1}+\Phi_{j,1,0})}{\partial{x}_{k}}W_{k,1}-\frac{\partial{W}_{j,1}}{\partial{x}_{k}}(\Phi_{k,0,1}+\Phi_{k,1,0})\right],$
(43)

whose evaluation only involves straightforward operations. Then, the second
order terms of the mean variations of $a$, $e$, $I$, $\Omega$, and $\omega$,
are chosen by averaging, to get

$\displaystyle\Phi_{1,0,2}=$ $\displaystyle\;0,$ (44)
$\displaystyle\Phi_{2,0,2}=$
$\displaystyle-n\frac{R_{\Earth}^{4}}{a^{4}}\frac{9}{16}\frac{s^{2}e}{1+\eta}\left[11\eta^{2}(3s^{2}-2)+3\eta(13s^{2}-10)+12(s^{2}-1)\right]$
$\displaystyle\times\sin 2\omega,$ (45) $\displaystyle\Phi_{3,0,2}=$
$\displaystyle\;n\frac{R_{\Earth}^{4}}{a^{4}}\frac{9}{16}\frac{1}{\eta^{2}}\frac{cse^{2}}{1+\eta}\left[11\eta^{2}(3s^{2}-2)+3\eta(13s^{2}-10)+12(s^{2}-1)\right]$
$\displaystyle\times\sin 2\omega,$ (46) $\displaystyle\Phi_{4,0,2}=$
$\displaystyle
n\frac{R_{\Earth}^{4}}{a^{4}}\frac{3c}{16\eta^{2}}\big{\\{}33\eta^{3}(43s^{2}-12)+104\eta^{2}(s^{2}-1)-81\eta(15s^{2}-4)$
$\displaystyle-116(2s^{2}-1)+6\left[11\eta^{2}(3s^{2}-1)+3\eta(13s^{2}-5)+6(2s^{2}-1)\right]$
$\displaystyle\times(\eta-1)\cos 2\omega\big{\\}},$ (47)
$\displaystyle\Phi_{5,0,2}=$
$\displaystyle\;n\frac{R_{\Earth}^{4}}{a^{4}}\frac{3}{32\eta^{2}}\big{\\{}33\eta^{3}(43s^{4}-86s^{2}+16)+52\eta^{2}(3s^{4}-6s^{2}+4)-162\eta$
$\displaystyle\times(s^{2}-1)(15s^{2}-4)-116(s^{2}-1)(5s^{2}-2)+6\big{[}11\eta^{3}(3s^{4}-6s^{2}$
$\displaystyle+2)+\eta^{2}(9s^{4}-16s^{2}+8)-18\eta(s^{2}-1)(3s^{2}-1)-6(s^{2}-1)$
$\displaystyle\times(5s^{2}-2)\big{]}\cos 2\omega\big{\\}}$ (48)

Finally, we compute $C_{1,2}=\langle{W}_{1,2}\rangle_{M}$ from Eq. (31), to
obtain

$\displaystyle C_{1,2}=$
$\displaystyle\;-a\frac{R_{\Earth}^{4}}{a^{4}}\frac{3}{32\eta}\big{\\{}21\eta^{3}(43s^{4}-24s^{2}+8)+104\eta^{2}s^{2}(s^{2}-2)$
$\displaystyle-15\eta(45s^{4}-24s^{2}+8)-232s^{2}(s^{2}-1)$
$\displaystyle+6(\eta-1)\left[7\eta^{2}(3s^{2}-2)+\eta(27s^{2}-22)-12c^{2}\right]s^{2}\cos
2\omega\big{\\}},$ (49)

from which
$\Phi_{6,0,2}=\langle\tilde{\Phi}_{6,0,2}\rangle_{M}-\frac{3}{2}(n/a)C_{1,2}$,
and hence

$\displaystyle\Phi_{6,0,2}=$
$\displaystyle\;n\frac{R_{\Earth}^{4}}{a^{4}}\frac{3}{64\eta}\big{\\{}459\eta^{3}(43s^{4}-24s^{2}+8)+1456\eta^{2}s^{2}(s^{2}-2)-315\eta$
$\displaystyle\times(45s^{4}-24s^{2}+8)-2784s^{2}(s^{2}-1)+6s^{2}\big{[}153\eta^{3}(3s^{2}-2)$
$\displaystyle+28\eta^{2}(3s^{2}-4)-105\eta(3s^{2}-2)-144(s^{2}-1)\big{]}\cos
2\omega\big{\\}}.$ (50)

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\geq 2$ we readily obtain $C_{j,2}=0$.

For instance, for the semimajor axis we obtain

$a_{0,2}=a\frac{R_{\Earth}^{4}}{a^{4}}\frac{3}{32}\frac{1}{1+\eta}\sum_{i=0}^{2}\sum_{j=1+\delta_{2,i}}^{2i+2}\sum_{k=\delta_{2,i}}^{2+2\delta_{2i,j}}P_{i,j,k}\eta^{k-1}e^{|j-2i|}s^{2i}\cos(2i\omega+jM),$
(51)

where $\delta_{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\geq 1$. On the
contrary, the inverse, osculating to mean correction to the semimajor axis
takes the form

$a_{0,2}^{\prime}=a_{0,2,\mathrm{long}}^{\prime}+a\frac{R_{\Earth}^{4}}{a^{4}}\frac{3}{16}\sum_{i=0}^{2}\sum_{j=1}^{2i+2}\sum_{k=1-\delta_{2,i}}^{1+2\delta_{2i,j}}P^{\prime}_{i,j,k}\eta^{k-1}e^{|j-2i|}s^{2i}\cos(2i\omega+jM),$
(52)

with the new inclination polynomials $P^{\prime}_{i,j,k}$ of Table 2, and the
long-period terms

$a_{0,2,\mathrm{long}}^{\prime}=a\frac{R_{\Earth}^{4}}{a^{4}}\frac{3}{16}\sum_{i=0}^{1}\sum_{k=0}^{3}P^{\prime}_{i,0,k}\eta^{k-1}s^{2i}\cos
2i\omega,$ (53)

with coefficients $P^{\prime}_{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
$\langle{a}_{0,2}^{\prime}\rangle_{M}=a_{0,2,\mathrm{long}}^{\prime}$. The
same features are indeed observed in the direct and inverse corrections to the
other orbital variables.

Table 1: Inclination polynomials $P_{i,j,k}$ in Eq. (51). ${}_{0,1,0}:-208s^{2}(s^{2}-1)$ | ${}_{1,2,2}:8(255s^{2}-142)$ | ${}_{2,3,1}:3$
---|---|---
${}_{0,1,1}:-4(289s^{4}-196s^{2}+48)$ | ${}_{1,2,3}:-8(405s^{2}-298)$ | ${}_{2,3,2}:-3$
${}_{0,1,2}:-8(125s^{4}-72s^{2}+24)$ | ${}_{1,2,4}:-968(3s^{2}-2)$ | ${}_{2,4,1}:-54$
${}_{0,2,0}:224s^{2}(s^{2}-1)$ | ${}_{1,3,0}:40(s^{2}-1)$ | ${}_{2,4,2}:-54$
${}_{0,2,1}:2(247s^{4}-40s^{2}-24)$ | ${}_{1,3,1}:2(281s^{2}-194)$ | ${}_{2,4,3}:70$
${}_{0,2,2}:2(191s^{4}+72s^{2}-24)$ | ${}_{1,3,2}:162(3s^{2}-2)$ | ${}_{2,4,4}:70$
${}_{1,1,0}:168(s^{2}-1)$ | ${}_{1,4,0}:-336(s^{2}-1)$ | ${}_{2,5,1}:-63$
${}_{1,1,1}:6(109s^{2}-82)$ | ${}_{1,4,1}:-84(7s^{2}-6)$ | ${}_{2,5,2}:-77$
${}_{1,1,2}:198(3s^{2}-2)$ | ${}_{1,4,2}:-140(3s^{2}-2)$ | ${}_{2,6,1}:-147$
${}_{1,2,0}:720(s^{2}-1)$ | ${}_{2,2,1}:-3$ | ${}_{2,6,2}:-147$
${}_{1,2,1}:72(43s^{2}-32)$ | ${}_{2,2,2}:-3$ |
Table 2: Non-null coefficients $P^{\prime}_{i,j,k}$ in Eqs. (52) and (53). ${}_{0,0,0}:-232s^{2}(s^{2}-1)$ | ${}_{0,1,1}:96(3s^{4}-3s^{2}+1)$ | ${}_{1,3,1}:-108(3s^{2}-2)$
---|---|---
${}_{0,0,1}:-15(45s^{4}-24s^{2}+8)$ | ${}_{0,2,1}:6(11s^{4}-24s^{2}+8)$ | ${}_{1,4,1}:-84(3s^{2}-2)$
${}_{0,0,2}:104s^{2}(s^{2}-2)$ | ${}_{1,1,0}:24(s^{2}-1)$ | ${}_{2,2,1}:3$
${}_{0,0,3}:21(43s^{4}-24s^{2}+8)$ | ${}_{1,1,1}:-12(3s^{2}-2)$ | ${}_{2,3,1}:-12$
${}_{1,0,0}:-72(s^{2}-1)$ | ${}_{1,2,0}:-24(s^{2}-1)$ | ${}_{2,4,1}:-30$
${}_{1,0,1}:-30(3s^{2}-2)$ | ${}_{1,2,1}:-96(3s^{2}-2)$ | ${}_{2,4,3}:42$
${}_{1,0,2}:12(3s^{2}-4)$ | ${}_{1,2,3}:72(3s^{2}-2)$ | ${}_{2,5,1}:84$
${}_{1,0,3}:42(3s^{2}-2)$ | ${}_{1,3,0}:-56(s^{2}-1)$ | ${}_{2,6,1}:147$

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
$\Phi_{j,k,0}=0$ for $k\geq 2$, from repeated iterations of recursion (11) we
obtain

$\begin{array}[]{l}\Phi_{j,0,3}=\mathcal{L}^{*}_{j,1}\left(\Phi_{j,0,2}\right)+\Phi_{j,1,2}\\\
\Phi_{j,1,2}=\mathcal{L}^{*}_{j,2}\left(\Phi_{j,0,1}\right)+\mathcal{L}^{*}_{j,1}\left(\Phi_{j,1,1}\right)+\Phi_{j,2,1}\\\
\Phi_{j,2,1}=\mathcal{L}^{*}_{j,3}\left(\Phi_{j,0,0}\right)+2\mathcal{L}^{*}_{j,2}\left(\Phi_{j,1,0}\right)+\Phi_{j,3,0},\end{array}$

where, from the same recursion,
$\Phi_{j,1,1}=\Phi_{j,0,2}-\mathcal{L}^{*}_{j,1}\left(\Phi_{j,0,1}\right)$.
Then, the computable terms are
$\tilde{\Phi}_{j,0,3}=\mathcal{L}^{*}_{j,1}\left(\Phi_{j,0,2}+\Phi_{j,1,1}\right)+\mathcal{L}^{*}_{j,2}\left(\Phi_{j,0,1}+2\Phi_{j,1,0}\right)$,
and, like before, we compute the third order components of the mean variations
as $\Phi_{j,0,3}=\langle\tilde{\Phi}_{j,0,3}\rangle_{M}$ for $j<6$, whereas
the computation of
$\Phi_{6,0,3}=\langle\tilde{\Phi}_{6,0,3}-\frac{3}{2}(n/a)W_{1,3}\rangle_{M}$
requires the preliminary computation of $C_{1,3}=\langle{W}_{1,3}\rangle_{M}$
from Eq. (32). For the latter we obtain

$\displaystyle\Phi_{1,0,3}=$
$\displaystyle\;an\left(\\!\frac{R_{\Earth}}{a}\\!\right)^{6}\frac{c^{2}s^{2}e^{2}}{\eta^{2}(1+\eta)}\frac{81}{16}\left[7\eta^{2}(3s^{2}-2)+\eta(27s^{2}-22)+12(s^{2}-1)\right]$
$\displaystyle\times\sin 2\omega.$ (54)

The non-vanishing character of the mean variation of the semimajor axis is a
crucial fact that remains hidden when the perturbation approach is truncated
to the lower orders. Indeed, due to the Lyapunov instability of the Keplerian
motion, the semimajor axis is the more sensitive element in the propagation of
errors. So neglecting these terms of the variation of the mean semimajor axis,
how tiny they may be, soon or later will introduce observable errors in the
mean elements propagation. This feature of Theory 1 will be clearly
illustrated in the semi-analytical propagations presented in Section 6.

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 $\kappa=e\cos\omega$, $\sigma=e\sin\omega$, 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
$\kappa=(e^{\prime}+\epsilon{e}_{0,1}+\frac{1}{2}\epsilon^{2}e_{0,2})\cos(\omega^{\prime}+\epsilon\omega_{0,1}+\frac{1}{2}\epsilon^{2}\omega_{0,2})$,
which, after expansion and truncation to $\mathcal{O}(\epsilon^{3})$ yields

$\displaystyle\kappa=$
$\displaystyle\;e^{\prime}\cos\omega^{\prime}+\epsilon\left(e_{0,1}\cos\omega^{\prime}-\omega_{0,1}e^{\prime}\sin\omega^{\prime}\right)+\frac{1}{2}\epsilon^{2}\left(e_{0,2}\cos\omega^{\prime}-\omega_{0,2}e^{\prime}\sin\omega^{\prime}\right)$
$\displaystyle-\frac{1}{2}\epsilon^{2}\left(2e_{0,1}\omega_{0,1}\sin\omega^{\prime}+\omega_{0,1}^{2}e^{\prime}\cos\omega^{\prime}\right),$
(55)

where the products $e_{0,1}\omega_{0,1}$ and $\omega_{0,1}^{2}$ 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.

##  5 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}=\langle{W}_{j,m}\rangle_{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
$\Phi_{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,

$\displaystyle\Phi_{6,0,2}=$
$\displaystyle\;n\frac{R_{\Earth}^{4}}{a^{4}}\frac{3}{32\eta}\big{\\{}198\eta^{3}(43s^{4}-24s^{2}+8)+572\eta^{2}s^{2}(s^{2}-2)-135\eta(45s^{4}$
$\displaystyle-24s^{2}+8)-1044s^{2}(s^{2}-1)+6s^{2}\big{[}66\eta^{3}(3s^{2}-2)+11\eta^{2}(3s^{2}-4)$
$\displaystyle-45\eta(3s^{2}-2)-54(s^{2}-1)\big{]}\cos 2\omega\big{\\}},$ (56)

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,\mathrm{short}}+a_{0,2,\mathrm{long}},$ (57)

where $a_{0,2,\mathrm{short}}$ is the same as Eq. (51), which only comprises
short-period terms, and the long-period terms $a_{0,2,\mathrm{long}}$ are
exactly one half of those reported in Eq. (53). Analogously, now

$a^{\prime}_{0,2}=a^{\prime}_{0,2,\mathrm{short}}+a^{\prime}_{0,2,\mathrm{long}},$
(58)

where the short-period terms $a^{\prime}_{0,2,\mathrm{short}}$ are those
obtained from Eq. (52) as $a_{0,2}^{\prime}-a_{0,2,\mathrm{long}}^{\prime}$,
whereas $a^{\prime}_{0,2,\mathrm{long}}=a_{0,2,\mathrm{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].

##  6 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^{\circ}$, and semimajor
axis $a=9500$ km. The gravity field parameters used in the propagations are
$\mu=398600.4415\;\mathrm{km^{3}/s^{2}}$, $R_{\Earth}=6378.1363$ km, 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-to-osculating 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 $(m+1)$-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
$(m+1)$-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^{\prime}_{0,2}$, cf. Eqs. (52)–(53) and Eq. (58), as well as in
$\Phi_{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.

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

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 $\sim 10^{-6}$ for the test case, and to values below
one tenth of arc second for the orbital angles.

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

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 $\mathcal{O}(J_{2})$ 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
$\mathcal{O}(10^{-8})$ relative to distance.

Figure 3: Test case. Position errors of the 2nd-order of Theory 1 (left) and 2
(right).

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 zero-average 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.

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

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).

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

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.

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)

Finally, we must mention that, beyond the illustration purposes of this
Section, both perturbation solutions of the toy model have been also obtained
in non-singular 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.

### Acknowledgements

The research has received funding from the Spanish State Research Agency and
the European Regional Development Fund (Projects PID2020- 112576GB-C22 and
PID2021-123219OB-I00, AEI/ERDF, EU)

## References

  * [1] D. A. Vallado.  Fundamentals of Astrodynamics and Applications.  Microcosm, 4940 West 147th Street, Hawthorne, CA 90250-6708 USA, 4th edition, 2013.
  * [2] Y. Kozai.  Second-order solution of artificial satellite theory without air drag.  The Astronomical Journal, 67:446–461, September 1962.
  * [3] D. Brouwer.  Solution of the problem of artificial satellite theory without drag.  The Astronomical Journal, 64:378–397, November 1959.
  * [4] Y. Kozai.  The motion of a close earth satellite.  The Astronomical Journal, 64:367–377, November 1959.
  * [5] S. L. Coffey, H. L. Neal, A. M. Segerman, and J. J. Travisano.  An analytic orbit propagation program for satellite catalog maintenance.  In K. Terry Alfriend, I. Michael Ross, Arun K. Misra, and C. Fred Peters, editors, AAS/AIAA Astrodynamics Conference 1995, volume 90 of Advances in the Astronautical Sciences, pages 1869–1892, P.O. Box 28130, San Diego, California 92198, USA, 1996. American Astronautical Society, Univelt, Inc.
  * [6] M. Lara.  Improving efficiency of analytic orbit propagation (IAC-21,C1,7,2,x65390).  In Proceedings of the 72nd International Astronautical Congress (IAC), Dubai, United Arab Emirates, 25–29 October 2021, pages 1–10. International Astronautical Federation (IAF), International Astronautical Federation (IAF), 2021.
  * [7] C. A. Wagner, B. C. Douglas, and R. G. Williamson.  The ROAD program.  Technical Report NASA-TM-X-70676, Goddard Space Flight Center, Greenbelt, Maryland, January 1974.
  * [8] G. Métris and P. Exertier.  Semi-analytical theory of the mean orbital motion.  Astronomy and Astrophysics, 294:278–286, February 1995.
  * [9] R. S. Bhat, R. B. Frauenholz, and Patrick E. Cannell.  TOPEX/POSEIDON orbit maintenance maneuver design.  In Astrodynamics 1989, pages 645–670, January 1990.
  * [10] J. R. Guinn.  Periodic gravitational perturbations for conversion between osculating and mean orbit elements (AAS 91-430).  In Bernard Kaufman, Kyle T. Alfriend, Ronald L. Roehrich, and Robert R. Dasenbrock, editors, AAS/AIAA Astrodynamics Specialist Conference, Durango, CO, Aug. 19-22, 1991, volume 76 of Advances in the Astronautical Sciences, pages 1–25, P.O. Box 28130, San Diego, California 92198, USA, 1992. American Astronautical Society, Univelt, Inc.
  * [11] Wayne D. McClain.  A Recursively Formulated First-Order Semianalytic Artificial Satellite Theory Based on the Generalized Method of Averaging, Volume 1: The Generalized Method of Averaging Applied to the Artificial Satellite Problem.  NASA CR-156782. NASA, Greenbelt, Maryland, 2nd edition, 1977.  Also, A Semianalytic Artificial Satellite Theory. Volume 1, 1992.
  * [12] D. A. Danielson, B. Neta, and L. W. Early.  Semianalytic Satellite Theory (SST): Mathematical algorithms.  Technical Report NPS-MA-94-001, Naval Postgraduate School, Naval Postgraduate School, Monterey, CA. Dept. of Mathematics., January 1994.
  * [13] M. Lara, J. F. San-Juan, Z. J. Folcik, and P. Cefola.  Deep Resonant GPS-Dynamics Due to the Geopotential.  Journal of the Astronautical Sciences, 58(4):661–676, October 2011\.
  * [14] Martin Lara, Juan F. San-Juan, Luis M. López-Ochoa, and Paul Cefola.  Long-term evolution of Galileo operational orbits by canonical perturbation theory.  Acta Astronautica, 94(2):646–655, February 2014.
  * [15] M. Lara, J. F. San-Juan, and L. M. López-Ochoa.  Proper Averaging Via Parallax Elimination (AAS 13-722).  In Stephen B. Broschart, James D. Turner, Kathleen C. Howell, and Felix R. Hoots, editors, Astrodynamics 2013, volume 150 of Advances in the Astronautical Sciences, pages 315–331, P.O. Box 28130, San Diego, California 92198, USA, January 2014. American Astronautical Society, Univelt, Inc.
  * [16] M. Lara, J. F. San-Juan, and L. M. López-Ochoa.  Delaunay variables approach to the elimination of the perigee in Artificial Satellite Theory.  Celestial Mechanics and Dynamical Astronomy, 120(1):39–56, September 2014.
  * [17] J. A. Morrison.  Generalized Method of Averaging and the Von Zeipel Method.  In Raynor L. Duncombe and Victor G. Szebehely, editors, Methods in Astrodynamics and Celestial Mechanics, volume 17 of Progress in Astronautics and Rocketry, pages 117–138. Elsevier, 1966.
  * [18] S. Ferraz-Mello.  Do Average Hamiltonians Exist?  Celestial Mechanics and Dynamical Astronomy, 73:243–248, January 1999.
  * [19] S. Coffey and K. T. Alfriend.  An analytical orbit prediction program generator.  Journal of Guidance, Control and Dynamics, 7(5):575–581, 1984.
  * [20] Z. Folcik and P. Cefola.  A General Solution to the Second Order $J_{2}$ Contribution in a Mean Equinoctial Element Semianalytic Satellite Theory.  In S. Ryan, editor, Advanced Maui Optical and Space Surveillance Technologies Conference, page 45, September 2012.
  * [21] M. Lara, J. F. San-Juan, and D. Hautesserres.  HEOSAT: a mean elements orbit propagator program for highly elliptical orbits.  CEAS Space Journal, 10:3–23, March 2018.
  * [22] A. H. Nayfeh.  Perturbation Methods.  Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, 2004.
  * [23] A. Deprit.  Canonical transformations depending on a small parameter.  Celestial Mechanics, 1(1):12–30, 1969.
  * [24] A. A. Kamel.  Perturbation Method in the Theory of Nonlinear Oscillations.  Celestial Mechanics, 3:90–106, March 1970.
  * [25] J. Henrard.  On a perturbation theory using Lie transforms.  Celestial Mechanics, 3:107–120, March 1970.
  * [26] D. Brouwer and G. M. Clemence.  Methods of Celestial Mechanics.  Academic Press, New York and London, 1961.
  * [27] G. Hori.  Theory of General Perturbations for Non-Canonical Systems.  Publications of the Astronomical Society of Japan, 23:567–587, 1971\.
  * [28] R. H. Lyddane.  Small eccentricities or inclinations in the Brouwer theory of the artificial satellite.  Astronomical Journal, 68(8):555–558, October 1963.
  * [29] A. Deprit and A. Rom.  The Main Problem of Artificial Satellite Theory for Small and Moderate Eccentricities.  Celestial Mechanics, 2(2):166–206, June 1970.
  * [30] R. H. Lyddane and C. J. Cohen.  Numerical comparison between Brouwer’s theory and solution by Cowell’s method for the orbit of an artificial satellite.  Astronomical Journal, 67:176–177, April 1962.
  * [31] J. V. Breakwell and J. Vagners.  On Error Bounds and Initialization in Satellite Orbit Theories.  Celestial Mechanics, 2:253–264, June 1970.
  * [32] M. Lara.  Brouwer’s satellite solution redux.  Celestial Mechanics and Dynamical Astronomy, 133(47):1–20, September 2021.
  * [33] M. Lara.  Note on the analytical integration of circumterrestrial orbits.  Advances in Space Research, 69:4169–4178, April 2022.


# A solution to the 7-prince’s tour problem

Lorenzo Mella  Dip. di Scienze Fisiche, Informatiche, Matematiche, Università
degli Studi di Modena e Reggio Emilia, Via Campi 213/A, I-41125 Modena, Italy
<EMAIL_ADDRESS>
###### Abstract.

In [1] the authors studied the closed tour problem on the $8\times 8$
chessboard of a chess piece, called $k$-prince, leaving open the existence of
such a tour when $k=7$. In this note we find a solution to this open case.

###### Key words and phrases:

Chessboard, Knight’s Tour, Prince’s Tour

###### 2020 Mathematics Subject Classification:

05C12, 05C45

##  1\. On the prince’s tour problem on a square chessboard

Chessboard related problems have always been a widely studied topic in
Discrete Mathematics, and often they arise from recreational purposes, see for
instance [3]. In particular, one of the most famous problems in this area is
the following:

######  The Knight’s Tour Problem.

Can a knight make a closed tour on an $8\times 8$ chessboard visiting each
cell exactly once?

The affirmative solution to this problem has been known for centuries, and
from there researchers considered knight’s tours over chessboards of different
shapes, or tours of different chess pieces.

In [1] the authors, starting from an irregular domination problem over grid
graphs, proposed the movement of a prince over a chessboard. Given a positive
integer $k$, a $k$-prince (in short $k$P) can move $k$ squares away from its
current position: namely, if the $k$-prince is in the $(i,j)$-th cell of a
chessboard, then it can move in the $(u,v)$-th cell if and only if

$|i-u|+|j-v|=k.$

Given two cells $(i,j)$ and $(u,v)$, we say that $(u,v)$ is attacked by
$(i,j)$ if a $k$P can jump directly from $(i,j)$ to $(u,v)$. Note that if
$(i,j)$ attacks $(u,v)$, then $(u,v)$ attacks $(i,j)$. Then, a $k$P’s closed
tour over an $m\times n$ board is a sequence
$(c_{1},c_{2},\dotsc,c_{mn})\in\\{[1,m]\times[1,n]\\}^{mn}$ such that:

  * (1)

$\bigcup_{i=1}^{mn}c_{i}=\\{[1,m]\times[1,n]\\}$;

  * (2)

$c_{i}$ attacks $c_{i+1}$ for each $i\in[1,mn-1]$;

  * (3)

$c_{mn}$ attacks $c_{1}$.

If Condition (3) is dropped, then the tour is said to be open. In [1], the
authors proposed the following problem:

######  The Prince’s Tour Problem.

Let $k$ be a positive integer. Can a $k$-prince make a closed tour on an
$8\times 8$ chessboard visiting each cell exactly once?

Also, they proved that the set of admissible values for $k$ is
$\\{1,3,5,7\\}$, and remarked that the existence of a $1$P’s closed tour is a
fairly easy exercise (indeed, it exists for every $2n\times 2n$ chessboard,
with $n\geq 1$), while a $3$P’s closed tour can be found by noticing that the
standard chessboard knight is a $3$P that cannot move purely horizontally or
vertically: since a closed tour of a knight exists, then a $3$P’s closed tour
exists. In Figure 10 of [1] the authors directly write a $5$P’s closed tour,
leaving as an open problem the existence of a $7$P’s closed tour. Here, we
show the $7$P’s closed tour that we found by considering rotations of a
smaller one:

$\begin{array}[]{|c|c|c|c|c|c|c|c|}\hline\cr 52&41&48&27&36&1&54&13\\\
\hline\cr 43&46&25&34&7&30&51&56\\\ \hline\cr 64&3&60&15&10&5&40&49\\\
\hline\cr 29&58&55&12&21&18&63&38\\\ \hline\cr 6&31&50&53&44&23&26&61\\\
\hline\cr 17&8&37&42&47&28&35&32\\\ \hline\cr 24&19&62&39&2&57&14&11\\\
\hline\cr 45&22&33&4&59&16&9&20\\\ \hline\cr\end{array}$

This result, combined with the ones of [1], gains the following:

######  Proposition 1.1.

There exists a $k$-prince’s tour on an $8\times 8$ chessboard if and only if
$k$ is odd and not greater than $7$.

We can moreover prove the following proposition:

######  Proposition 1.2.

For $n\in[1,4]$ there exists a closed $k$P’s tour over a $2n\times 2n$
chessboard for every odd $k\leq 2n-1$.

###### Proof.

As already remarked, there exists a closed $1$P’s tour on a $2n\times 2n$
chessboard for every integer $n\geq 1$. If $n=2$, then we need to check the
existence of a closed $3P$’s tour over a $4\times 4$ chessboard (we remark
that a closed knight’s tour on this chessboard does not exist, see [3]); an
example of such a tour is the following:

$\begin{array}[]{|c|c|c|c|}\hline\cr 3&6&9&14\\\ \hline\cr 12&15&4&7\\\
\hline\cr 5&2&13&10\\\ \hline\cr 16&11&8&1\\\ \hline\cr\end{array}$

If $n=3$, then we need to check the statement for $k=3,5$: since there exists
a closed knight’s tour then there exists a $3$-prince’s closed tour, while an
example for $k=5$ is the following one:

$\begin{array}[]{|c|c|c|c|c|c|}\hline\cr 9&12&5&24&29&18\\\ \hline\cr
20&17&34&31&26&21\\\ \hline\cr 15&22&1&10&7&14\\\ \hline\cr
32&25&28&19&4&33\\\ \hline\cr 3&8&13&16&35&2\\\ \hline\cr 36&11&6&23&30&27\\\
\hline\cr\end{array}$

The case $n=4$ follows from Proposition 1.1. ∎

Moreover, in [2] the author analized the closed tour problem of generalized
knights over chessboards: given a positive integer $k$, the generalized
$k$-knight is the piece that is allowed to move along either the horizontal or
vertical direction by $k-1$ steps, and along the other direction by $1$ step.
In particular, it was proven that for every $k\geq 1$ there exists a closed
tour of a generalized $k$-knight on a $2n\times 2n$ chessboard for
sufficiently large values of $n$. Here, we consider Theorem 3 of [2], where we
slightly modify the statement in the context of princes:

######  Theorem 1.3.

Let $n$ be a positive integer, and let $k$ be any odd integer such that
$n\geq(k-1)(3k-2)$. Then, there exists a $k$P’s tour over a $2n\times 2n$
chessboard.

In view of the results contained in [1] and the ones presented in this note,
the following question naturally arises:

######  Open Problem.

Does there exist a closed $k$P’s tour over a $2n\times 2n$ chessboard for
every odd $k\leq 2n-1$?

## Acknowledgements

The author was partially supported by INdAM-GNSAGA.

## References

  * [1] G. Chartrand and P. Zhang, A chessboard problem and irregular domination, Bull. Inst. Combin. Appl. 98, 43–59 (2023).
  * [2] N. Kamčev, Generalised Knight’s Tours, Electron. J. Combin. 21, P1.31 (2014).
  * [3] J.J. Watkins, Across the board: The Mathematics of Chessboard Problems, Princeton university press (2012).