Long-term evolution of mid-altitude quasi-satellite orbits

Quasi-satellite orbits are of great interest for the exploration of planetary moons because of their dynamical features and close proximity with respect to the surface of scientifically relevant objects like Phobos and Deimos. This paper explores the equations of the elliptical Hill problem, offering a new analytical insight into the long-term evolution of mid-altitude quasi-satellite orbits. Our developments are based on the Yamanaka–Ankersen solution of the Tschauner–Hempel equations and capture the effects of the secondary’s gravity and orbital eccentricity on the shape and orientation of near-equatorial retrograde relative trajectories. The analytical solution of the in-plane and out-of-plane components of the secular motion is achieved by averaging over the relative longitude of a spacecraft as seen from the co-rotating frame of the two primaries. Developments are validated against the numerical integration of quasi-periodic trajectories that densely cover the surface of three-dimensional invariant tori. This analysis confirms the stable nature of quasi-satellite orbits and provides new tools for future spacecraft missions such as the Martian Moons eXploration envisaged by JAXA.


Introduction
The study of retrograde relative trajectories around the secondary of two attracting masses has always been of interest to astronomers and orbital mechanicists since their original discovery in 1897 [1]. Already in 1913, Jackson was theorizing that this type of orbits may be more stable than direct trajectories [2], a result that was eventually confirmed by Hénon in his numerical investigations of the restricted three-body problem [3]. Based on the results of Hénon's analysis, retrograde relative orbits, also known in the literature as distant retrograde or quasi-satellite orbits (QSO), have gained the attention of several space agencies and astrodynamics researchers looking for efficient trajectories to explore the solar system. Many spacecraft missions have recently been proposed to fly on retrograde relative orbits around planetary moons, including the Russian Phobos-Grunt, NASA's ARM, and the Euro-pean DePhine [4][5][6]. JAXA is also planning the Martian Moons eXploration mission (MMX) in order to finally settle the debate on the moons' origin and further deepen our understanding of the Martian system. The spacecraft will be placed on quasi-satellite trajectories near the largest of the Martian moons Phobos, studying its dynamical and geophysical features for more then 3 years [7].
In order to better understand motion near small irregular moons and prepare for the proximity operations of MMX, this paper focuses on the design and longterm evolution of quasi-satellite orbits. Recent developments on this matter can be separated in numerical and analytical investigations of a point mass subject to the gravitational attraction of a massive primary and a much lighter secondary. In the case of numerical investigations, shooting techniques are typically applied to calculate periodic or quasi-periodic solutions in the vicinity of planetary satellites [8][9][10]. Once these particular solutions are obtained, stability analysis can be carried out using Floquet theory and/or chaos indicators [11][12][13][14]. As for analytical developments, disturbing functions and variational equations are typically averaged over the time scales of fast variables to obtain qualitative information on the system dynamics. Belonging to the first category are the contributions inspired by Mikkola and Innanen, who first considered quasi-satellite trajectories as slightly perturbed Keplerian orbits in a 1:1 mean motion resonance [15]. Building on this idea, Namouni, Brasser et al., and Mikkola et al. later found out that transitions between QSOs and other 1:1 resonant orbit types are not only possible but also frequent [16][17][18]. Necessary conditions for transitioning between different 1:1 orbital regimes were recently developed by Sidorenko et al. using Hamiltonian theory and double-averaging [19]. A second approach is to consider quasi-satellite orbits as perturbed relative trajectories that deviate from 2:1 ellipses because of the gravitational attraction of the secondary body. Starting from the equations of the circular restricted three-body problem (CRTBP), Kogan used linearization and averaging to prove the stability of QSO orbits when Phobos is in a circular orbit around Mars [20]. Lidov and Vashkov'yak extended this approach to the elliptical three-body problem, but failed to come up with an analytical solution that was valid in the orbital regime of the Martian moon [21]. Owing to small mass ratios and QSO altitudes, Lara more recently considered Lie-Deprit transformations within the framework of the planar circular Hill problem [22]. Differently from the CRTBP, the Hill problem does not depend on any external parameter and is the ideal laboratory for dynamical investigations of many planetary systems. Cabral recognized this advantage and included eccentricity and out-of-plane oscillations to increase the fidelity of his analytical developments [23]. Unfortunately, despite thorough derivations with the variation of parameters, the differential equations of QSO amplitudes, offsets, and phase angles do not seem to adequately capture the long-term evolution of point masses in this orbital regime.
In this paper, we investigate the long-term evolution of mid-altitude quasi-satellite orbits within the framework of the spatial elliptical Hill problem (EHP). Our developments are also based on the analytical solution of the Tschauner-Hempel equations derived by Yamanaka and Ankersen in [24], but yields different equations than the ones found in Cabral's thesis. The analytical derivations are validated by comparison with the numerical integration of initial conditions on the surface of three-dimensional quasi-periodic invariant tori [25]. These invariant manifolds are obtained through the invariant tori of a stroboscopic mapping [26,27] and well represent the equivalent of QSO orbits whenever eccentricity and cross-track oscillations are taken into account. Next, we average over one QSO period and derive a new set of equations of motion for the long-term evolution of satellites in this orbital regime. An analytical solution of these equations is then achieved. A near-identity transformation is implemented to accurately map averaged elements into their osculating counterpart. The solution is finally validated by comparison with the numerical integration of the osculating equations of motion.
The article is organized as follows. First, dynamical models and QSO are introduced in Sect. 2. A novel analytical solution for both the in-plane and out-of-plane components of motion is achieved in Sect. 3 after averaging the dynamical system over one orbital period. The near-identity transformation is presented in Sect. 4, along with the numerical validation of the analytical solution.

Dynamical models
Dynamical models and assumptions necessary for the analytical developments of Sect. 3 are hereby intro-duced. Starting from the equations of the elliptical restricted three-body problem (ERTBP), we linearize the dynamics in the vicinity of the secondary body and discuss validity regions for the equations of the EHP. Examples of quasi-satellite orbits are then illustrated along with their dynamical substitutes when the eccentricity of planetary moon is taken into account (Sect. 2.1). In Sect. 2.2, we momentarily neglect the gravitational attraction of the secondary body to overview the analytical solution of the unperturbed problem as developed by Yamanaka and Ankersen in [24]. We then introduce Gauss variational equations (GVEs) in Sect. 2.3 to discuss the effects of the moon's gravity on the solution of the integrable problem. The results of this investigation are validated in Sect. 2.4 by comparison with the numerical integration of the EHP over 100 orbital revolutions (i.e., approximately one month in the Mars-Phobos system).

Elliptical Hill problem
Consider the motion of a massless satellite subject to the gravitational attraction of a moon and its host planet. As seen from the barycenter of the secondary body, the dynamics of the spacecraft is better expressed in a rotating reference frame S centered on the moon and such thatx is constantly aligned with the separatrix between the two primaries,ẑ is parallel to the orbital angular momentum of the moon, andŷ =ẑ ×x completes the right-handed triad. The equations of motion are given by the equations of ERTBP: where dots denote differentiation with respect to the true anomaly of the planetary satellite, ν; r andṙ are the position and velocity components of the spacecraft in S , respectively; I is the identity matrix; [ẑ] is a skew-symmetric matrix such that [ẑ] r =ẑ × r; and is the effective potential, with G, m, and M as the gravitational constant, mass of the secondary, and mass of the primary, respectively. Also note that is the angular velocity of the rotating reference frame (alongẑ), whereas r 12 = p γx is the relative position vector of the secondary as seen from the barycenter of the primary body. The variables p = a (1 − e 2 ) and γ = (1 + e cos ν) depend on the semimajor axis and eccentricity of the moon, hereby referred to as a and e. Finally, r 12 = r 12 = p γ −1 and r 1 = r 1 , where r 1 = r + r 12 is the relative position vector of the spacecraft as seen from the barycenter of the primary. It is assumed that the mass of the secondary is significantly smaller than the mass of the primary (m/M << 1), as well as that the relative distance between the spacecraft and the secondary is much smaller than the distance between the two attractors (r << r 12 ). Based on these assumptions, the term proportional to r −1 1 in Eq. (2) can be expanded up to the second order in r to obtain Moreover, since ω 2 G M p −3 γ 4 = G M γ r −3 12 , the effective potential may be rewritten as where = (m/M) 1/3 is the cubic root of the mass ratio parameter. Figure 1 illustrates the logarithm of |V − U|/U for several three-body systems and mass ratio parameters when e = 0. Red contours mark the boundary of a region where the relative error between the three-body problem and the so-called Hill approximation of the effective potential is four orders of magnitude smaller than the absolute value of U. As a result, the truncation errors can be considered negligible as long as (m/M) < 10 −4 and r < 0.05 a.
Within the validity region of the Hill approximation, the equation of motion (1) can be further simplified by introducing pulsating coordinates r = p γ −1r . This yields the first-order system of ordinary differential equations in ν known as the EHP [28]: Only the periodic orbits whose period is resonant with 2 π survive when e = 0 Here, x, y, z and u, v, w are the dimensionless position and velocity components ofr andṙ = γ p ṙ +γ γ r , respectively, whereas g = −r γr 3 stands for the gravitational attraction of the planetary moon in pulsating normalized units.
As it can be seen from Eq. (5), the EHP depends explicitly on the independent variable ν through γ = (1 + e cos ν). Because of this dependency, periodic orbits must be resonant with 2 π and cannot be organized in families as in the circular case [29]. This change of paradigm is illustrated in the chart of Fig. 2, showing the orbital period of QSO of the Mars-Phobos system as a function of their positive x-axis crossing when e = 0. Among the infinite family members, only those whose period is resonant with 7.66 hrs-the orbital period of Phobos around Mars-survive when the eccentricity of Phobos, e = 0.0151, is taken into account (i.e., intersections with the horizontal lines of Fig. 2).
When the resonant condition is not met, periodic orbits are replaced by two-dimensional quasi-periodic invariant tori that belong to a family of quasi-periodic retrograde orbits [30]. To compute these manifolds, we adapt the numerical algorithm outlined in Gómez & Mondelo [26] and Olikara & Scheeres [27] and hereby referred to as "GMOS". The main idea of GMOS is that quasi-periodic invariant tori can be calculated via invariant curves of a stroboscopic mapping by solving the boundary value problem described in [31]. The interested reader may find more details on the methodology in the authors' original papers as well as in [31,32]. In particular, Olikara provides an exhaustive explanation of the collocation method used in this work to march along the QSO family branch and generate a variety of quasi-periodic invariant tori with e = 0.0151 [32]. Several examples of planar quasi-periodic QSO orbits are shown in Fig. 3.
An advantage of the GMOS algorithm is that additional frequencies can be included to search for higherdimensional solutions of Eq. (5). This feature was demonstrated by Baresi and Scheeres in [25] and is hereby adopted for calculating quasi-periodic QSOs that extend beyond the equatorial plane of Phobos. We refer to this type of orbits as "3D QSOs" and discuss an example of them in Sect. 2.4.
Despite the GMOS algorithm succeeds in generating three-dimensional QSOs, the numerical procedure remains cumbersome and computationally expensive. Accordingly, the purpose of this paper is to study the long-term evolution of 3D QSOs via GVE and the averaging process so as to -gain analytical insight into the dynamical environment of planetary moons; -provide high-quality initial guesses for trajectory optimization without GMOS; -support onboard spacecraft operations with averaged analytical solutions to the relative motion problem.

General solution of the Keplerian relative motion problem
The EHP admits two equilibrium points in correspondence ofr * = ±(1/3) (1/3) 0 0 T . These points of equilibrium are often referred to as Lagrangian points and delimit a region of the phase space where the dynamics is dominated by the gravity field of the secondary body. Since the magnitude of the gravitational force due to the secondary is proportional to the inverse squared distance between its barycenter and the spacecraft, the dynamics inside the red curves of Fig. 1 remains dominated by the primary body as long as r * = p γ −1r * is significantly smaller than 0.05 a. As a preliminary rule of thumb, let us consider the cases where In what follows, we consider the Mars-Phobos case as an example of a large class of three-body systems that comply with the assumptions of the elliptical Hill problem while satisfying the negligible mass conditions given by Eq. (6). Our goal is to describe the long-term evolution of mid-altitude quasi-satellite orbits crossing the positive x-axis between 0.01 a and 0.025 a. We note that this corresponds to 93.77 km and 234.425 km for the Mars-Phobos system, which is the orbital region under consideration for the proximity operations of the Martian Moons eXploration mission [7].
Assuming that the gravitational attraction of Phobos is negligible, the dynamics of a satellite near the planetary satellite would be governed by the simplified set of equations given by System (7) is equivalent to the Tschauner-Hempel equations describing the relative motion of two neighboring satellites in eccentric Keplerian orbits and was analytically solved by Yamanaka and Ankersen [24,33]. The authors provide an expression for the relative position and velocity vectors in terms of the six integrals of motion K 1 , K 2 , K 3 , K 4 , K 5 , K 6 . The solution reads as where s * = cos ν + e cos 2 ν, c * = −(sin ν + e sin 2 ν), and J is an integral term defined as J = ν ν 0 1/γ (τ ) 2 d τ .
Equation (8) can be easily rewritten in matrix form as T . We note that the deter- , so that the matrix is nonsingular for any e ∈ [0, 1). It follows that the values of K can be inferred from rectangular coordinates through the inverse transformation [34] Following Cabral [23], auxiliary variables can be now introduced to gain geometrical insight into the analytical solution of the Tschauner-Hempel equations. Specifically, let Equation (11) reveals that a mass particle in orbit around a massless Phobos would remain in the vicinity of the Martian moon if A, δ x , and δ y have no secular growth. This is guaranteed when K 4 = 0 as A y and δ y would both remain constant regardless of the integral term appearing in Eqs. (10b) and (10d). If K 4 = 0, the trajectory described by the spacecraft in the equatorial plane becomes a pulsating ellipse with frequency ω xy = 1 and semimajor and semiminor axes equal to (1 + γ ) A and γ A, respectively. The phase angles α and β represent the longitude and latitude of the satellite at periarion, i.e., when ν = 0, as specified in Fig. 4. Meanwhile, the out-of-plane motion is decoupled from the planar one and follows an harmonic oscillation with amplitude B and frequency ω z = 1.

Gauss variational equations
As specified in Eq. (9), values of K can be inferred from the Cartesian coordinates of a spacecraft while knowing the true anomaly of Phobos around Mars. In consequence, the total differential of K = K (r, v, ν) reads aṡ and must be equal to zero when the dynamics is governed by the unperturbed system (7), wherev = a = Let us now consider the perturbed dynamics where u = u x , u y , u z T is a vector of disturbing forces per unit mass in the pulsating synodic reference frame S . Substituting Eq. (13) in Eq. (12) yields the GVE for the K integrals of motion given byK = ∂ K ∂v u [34].
The rates of change for orbit elements A, α, δ x , δ y , B, and β can be deduced from the application of the chain rule and yield We note that the rate of change of the latitude at periarion is inversely proportional to the magnitude of the out-of-plane motion B. To avoid singularities, we consideṙ instead of Eqs. (14e) and (14f) when studying the dynamical evolution of planar QSO orbits such as the one in Fig. 4. The final set of equations of motion for the osculating orbit elements where M is the phase space of the orbit elements oe, and G(ν, oe) is a six-by-three matrix with nonzero components Equation (16) can be now integrated to obtain the evolution of the osculating orbit elements without converting from the Cartesian coordinates adopted in Eq. (13). The system of ordinary differential equations will be referred to as GVE and validated against the numerical integration of the EHP in the next section. Furthermore, analytical expressions for the effects of Phobos' gravity can be now derived and analyzed. The result of these investigations will produce a linear model that will be averaged and solved analytically in Section 6 to grasp the secular evolution of the QSO mean relative orbit elements.

The effects of Phobos' gravity
Let us replace the terms u x , u y , u z appearing in (16) with the gravitational attraction of a spherical object with mass m significantly smaller than the mass of the primary (M >> m) in the pulsating reference frame S : where r (ν) = x(ν) 2 + y(ν) 2 + z(ν) 2 . Based on Eq. (11), r changes as a function of the relative orbit elements oe through a fairly complicated expression that yields the instantaneous relative distance between the satellite and the barycenter of the secondary body, e.g., Phobos. Nevertheless, a satellite on a QSO orbit remains close to the equatorial plane of the Martian moon if δ x and δ y are close to zero. In addition, the eccentricity of the Phobos is relatively small and allows for simplifications in the expressions of the GVE presented in Sect. 3. These considerations justify the first-order expansion of the equations of motion in e, A under the assumption that δ, χ, η 5 , η 6 O(e) << 1. This assumption was investigated by Cabral in [23] and seems to be reliable for mid-altitude QSO orbits based on numerical experiments. In this orbital regime, the expressions of x, y, and z in Eq. (11) imply r (e) A d 2 + 2 e cos ν cos 2 θ + 2 sin 2 θ where e = e, δ, χ T , d = cos 2 θ + 4 sin 2 θ , and Replacing the u x , u y , and u z terms of Eq.
Similarly, the other GVEs becomė oṙ oe = g(ν, oe) (23) in vectorial form. Equation (23) summarizes a linear model (LM) in e, δ, χ , η 5 , and η 6 that accurately predicts the dynamics of mass particles in mid-and high-altitude QSO orbits for several orbital periods around Phobos. This claim is supported by Figs. 5 6 7 and 8, which compare trajectories of the original and linearized systems obtained by numerical integration of the same set of initial conditions: in Cartesian coordinates and osculating orbit elements, respectively, with ν 0 = 324.8780 deg. Table 1 summarizes the models and relevant assumptions introduced throughout these subsections. As depicted in Fig. 8a, the error between the integration of Eqs. (16) and (13) is purely numerical and negligible after 100 orbital revolutions of Phobos around Mars (i.e., almost 32 days). In contrast, the position error of the linearized model exhibits a secular growth, but never exceeds 5 kilometers over the same interval of time. Such an error is deemed reasonable for the purpose of this article, as navigation errors and mismodeled dynamics would likely overtake the inaccuracies of the linearized model over the time span of a few days. Accordingly, the LM is considered as a good candidate for describing the long-term evolution of satellites in mid-altitude QSO orbits and can be used for the analysis of the secular motion as described in the following section. Fig. 6 Variables A, α, δ x , δ y obtained with the LM (blue) and GVEs (red) models

Averaged equations
The long-term evolution of quasi-QSO orbits is now analyzed by averaging Eq. (23) over one orbital period around Phobos, i.e.,oe = 2 π 0oe dθ . We find that the averaged orbit elementsĀ,ᾱ,δ x ,δ y ,K 5 , andK 6 evolve with the true anomaly aṡ where K 2.156516 and E 1.211056 are the complete elliptical integrals of the first and second kind, respectively, evaluated with modulus k = √ 3/2 [35].
We note here that our expressions approach the ones of Kogan [20] whenever e → 0, but differ from the ones found in Cabral's thesis [23]. There, the author derived secular equations using Lagrange planetary equations on the averaged disturbing potential due to the gravity field of Phobos, omitting the − 3 2δ x term appearing in Eq. (25d). Coefficients that model the long-term evolution of the relative orbit elements are also presented differently and in terms of different elliptical integrals, leading to erroneous frequencies of motion.
Equation (25) exhibits cross-coupling between some of the orbit elements, as well as long-term oscillations due to the secular growth of α, i.e., the longitude at periarion. Such a growth causes the geometry between the spacecraft, Mars, and Phobos, to change at every periarion passage, i.e., when the gravitational attraction of Mars is at its strongest. To grasp the secular effects  Fig. 9 Comparison between the numerical value of n QSO and its theoretical prediction 1 +α. The time unit T U is equal to 1.288 h of this phenomenon, observe thatᾱ grows linearly with the true anomaly and with a rate ω α = (K/π )Ā −3 0.686440Ā −3 that depends on the amplitude of the inplane motion: whereᾱ 0 =ᾱ(ν 0 ), and ν 0 is the true anomaly of Phobos at epoch. Observe that ω α is constant due to Eq. (25a). Furthermore, the value of ω α defines the period of the QSO orbit via T = 2 π/n QSO and n QSO =θ = ν +α 1+ω α . This relationship is demonstrated with the chart of Fig. 9, illustrating the comparison between the analytical prediction of n QSO and the numerical value computed with the GMOS algorithm of Sect. 2.1. We find that the proposed linearization and averaging approach are valid as long asĀ > 3.36, corresponding to 80 km for the Mars-Phobos system.
In the validity regime of Eq. (25), the in-plane and out-of-plane components of motion can be decoupled and reduced to the linear systems where and Using the Laplace transform, we find that the secular evolution ofδ x andδ y becomes whereν = (ν − ν 0 ),δ x,0 =δ x (ν 0 ),δ y,0 =δ y (ν 0 ), and The motion of the center of the in-plane ellipse is therefore characterized by two frequencies whose ratio is given by ω α /ω d 1.251283Ā −1.5 . In the limiting case e → 0, such a center traces out ellipses in the (δ x ,δ y ) plane with period 2 π/ω d . The maximum deviations inδ x andδ y occur at (ω dν )δ x = arctan D xδy,0 ω dδx,0 Within this region, the eccentricity of the secondary body plays an important role by introducingδ x and δ y oscillations regardless of the initial conditionsδ x,0 , δ y,0 ,ᾱ 0 . To quantify these effects, we investigate the magnitudes of the terms proportional to the eccentricity in Eq. (31) over a uniform grid ofĀ ∈ [3.36 , 7] values (corresponding to approximately 80 and 160 km for the Mars-Phobos system). The results of this analysis are portrayed in Fig. 12 and reveal thatδ y /Ā oscillations can be as large as 4.45 e even whenδ x,0 =δ y,0 = 0. In the worst-case scenario,ᾱ 0 = ±π/2, the validity of Eq. (31) holds as long as the eccentricity of the secondary is no larger than 0.022. In contrast, whenᾱ 0 = 0 or π , librations ofδ y /Ā do not violate the linearization assumptions of Sect. 5 for eccentricities up to 0.044.
Assuming the validity of Eq. (31) withδ x,0 =δ y,0 = 0, the eccentricity of the secondary body causes quasiperiodic oscillations inδ x andδ y that densely cover the surface of a two-dimensional invariant torus with frequencies ω α < ω d . The general case (e,δ x,0 ,δ y,0 = 0) is a linear combination of the two limiting scenarios in which the center of the in-plane relative trajectories remains bounded with respect to the origin of the coordinate system. An interesting situation occurs for Fig. 11 Maximumδ y deviations as a function of the initial conditionsδ x,0 ,δ y,0 when e = 0 Fig. 12 Relative magnitude of the terms proportional to the eccentricity of the secondary body as a function ofĀ. The top panel refers to terms inδ x , whereas the bottom panel is for terms inδ y those values ofĀ for which ω d /ω α is rational. In this resonant case, the motion of the center becomes fully periodic in theδ x ,δ y plane as illustrated in Fig. 13d.
Concerning the out-of-plane component, we observe that neither of Eqs. (22d), (22e), (25e), and (25f) depend on the eccentricity of Phobos' orbit. The same situation occurs when considering the expressions ofḂ andβ instead of Eqs. (22d) and (22e): where ϕ = (α − β) is the relative phase between the longitude and latitude at periarion. By averaging over one orbital period, the differential equations ofB andβ read aṡ respectively. The true anomaly rate of change ofφ becomesφ and admits as an analytical solution with By plugging Eq. (37) in Eq. (35), we find where . This result illustrates that the amplitude of the out-of-plane motion oscillates with a frequency of 2 ω ϕ 1.4 ω α and is bounded byB max =B 0 1.20859B 0 . Accordingly, in order to satisfy the linearization assumptions of Sect. 5,B 0 must be less than or equal to 0.08274Ā (corresponding to approximately 8 km in the Mars-Phobos system whenĀ = 4.18 100 km). Also notice thatβ varies with a time-varying frequency that is always smaller than ω α . Moreover, analytical expressions for the evolution ofK 5 andK 6 can be finally derived: and K 6 (ν) =B(ν) sinβ(ν) respectively. We note that both of Eqs. (31) and (40) depend on sines and cosines multiplied by time-periodic amplitudes at worst. This observation confirms that 3D quasi-QSO orbits are stable for as long as the assumptions of the LM model are valid. To validate this hypothesis, we apply the analytical solution of System (25) to obtain the dynamical evolution of a spacecraft under the assumption thatoe(ν 0 ) = oe(ν 0 ). Figure 14 shows the behavior of the mean orbit elementsδ x andδ y obtained from the osculating orbit elements of Eq. (24) over 400 orbital revolutions. Although the secular evolutions of δ x andδ y differ quite significantly from the osculating values, we find that the trajectory of the spacecraft does not drift away from the vicinity of Phobos' orbit. This shows that, in the averaged system, trajectories are bounded even whenδ x is not exactly zero. At the same time, we can also conclude that an appropriate transformation between osculating and mean orbit element is mandatory in order to adequately capture the secular evolution of spacecraft in mid-altitude QSOs. The problem of finding accurate initial conditions for the averaged system is addressed in the next section along with the numerical validation and analysis of Eqs. (31) and (40).

Mean-to-osculating orbit element mapping
In order to provide a detailed mapping between mean and osculating orbit elements, different methods of perturbation theories can be considered (see for instance [36], [22] and references therein). We choose to introduce the near-identity transformation T : [0, 2 π ] × M → M so as to reconstruct the unbiased oscillations of motion of the original system from the averaged trajectories. That is, where is a formal small parameter (consider A −2 for practical purposes) emphasizing the slow evolution of the mean orbit elementsoe with respect to the fast variables ν. Under this assumption, the averaged dynamics given by Eq. (25) may be rewritten in vectorial form asoe = ḡ (oe), whereaṡ in agreement with Eq. (23) and the formal definition g(ν, oe) = g (ν, oe). From the true anomaly rate of change of Eq. (41), it now follows that g (ν,oe + T (ν,oe)) = ḡ (oe) which implies up to the first order in . To obtain unbiased oscillations, let us also impose the integral constraint so that the analytical solution of the problem (44)-(45) may be written as [37]: Equation (46) is a Fourier series that depends on the Fourier coefficients g (k) of g (ν,oe) −ḡ (oe). The series is practically truncated at the N -th order in computer implementations, where N is an adequate high number. We highlight that transforming δ x first is mandatory for the problem at hand. On the one hand, the oscillations ofδ x are -small as long as the satellite remains on a mid-altitude 3D QSO orbit like the one in Fig. 5. This argument can be verified by inspection of the osculating and averaged evolution of δ x portrayed in Fig. 6. On the other hand, the differential equation ofδ y contains a linear term inδ x that is not multiplied by -small coefficients, i.e., 3 2δ x . Hence, an -small error would be introduced in the -slow dynamics ofδ y ifδ x is not properly handled (i.e., if oscillations of δ x are biased with respect toδ x at epoch).
To improve the accuracy of the near-identity transformation, we propose the nested transformation where T δ x and T δ y represent the δ x and δ y components of the T map. The advantages of the nested approach are illustrated in the plots of Fig. 15, portraying the reconstruction of δ y fromδ y with and without the preemptive transformation of δ x (Fig. 15a, b, respectively).
As it can be seen, the nested transformation performs much better than the classic one by accurately reconstructing the short-periodic oscillations of δ y up to 25 Phobos revolutions around Mars (8 days circa). In contrast, the biased value ofδ x (ν 0 ) at the epoch of Fig. 15b causes larger oscillations inδ y that translate into inaccurate osculating values. Based on this result, we apply the nested approach to calculate an approximate value ofoe and compare the numerical integration of (25) with Eqs. (31) and (40). The averaged initial conditions obtained from Eq. (24) arē and yield the relative errors of Fig. 16. We find that the relative error is negligible and use this information to validate our analytical developments. We can also investigate the long-term evolution of mid-altitude QSO orbits and compare it with the osculating orbit elements obtained in Section 2.4. Starting from the in-plane component of motion, it is found that satellites in orbit around Phobos trace out prolate  In-plane offset evolution over 100 Phobos orbits ellipses whose center oscillates in the δ y direction. This is verified in the plot of Fig. 17, illustrating the evolution of a typical QSO orbit in the δ x -δ y plane along with the analytical solution presented in Eq. (31). Similar comparisons can be carried out for the values of K 5 and K 6 by noting that the out-of-plane component of motion should follow a pulsating circle with radiusB(ν). This behavior is demonstrated in the chart of Fig. 18 after integrating the initial conditions of Fig. 5 for 400 orbital periods of Phobos (128 days circa). The numerical simulation provides us with the time histories of the α, β, and ϕ angles disclosed in Fig. 19, thus proving that the true anomaly rate of change of the α and ϕ angles is also in good agreement with the analytical values derived in Eqs. (26) and (38), respectively.
Finally, Fig. 20 shows the osculating and averaged values of the in-plane amplitude A, thereby demonstrating that the secular evolution of mid-altitude 3D QSOs can be adequately captured by Eq. (46).

Conclusions
This paper investigated quasi-satellite orbits around small secondary bodies to gain analytical insight into the elliptical restricted three-body problem. Starting from gravitational perturbation of the Tschauner-Hempel equations, we considered small eccentricities, in-plane, and out-of-plane oscillations to simplify the expressions of the QSO relative orbit elements rates of change with respect to the true anomaly of the secondary body. The linearized model obtained through these developments depends on the relative longitude between the spacecraft and the secondary attractor, which is a fast variable comparing to QSO orbit amplitude, offsets, and phase angles at epoch. Based on this observation, the equations of motion were averaged over one orbital period and solved analytically. We found that the eccentricity of the secondary body plays a key role by triggering librations in the radial and along-track directions of the primaries co-rotating frame. It was also verified that the amplitude of the relative motion remains constant as long as the assumptions of the linearized model are valid. These results extend the stability of near-equatorial retrograde orbits to the elliptical Hill problem, thereby validating the selection of QSO orbits for the proximity operations of future spacecraft missions to eccentric planetary moons (e.g., the Phobos sample retrieval MMX mission).
Future work will explore the advantages offered by the analytical solution of the averaged relative motion problem to the mission design and operations around small eccentric bodies. Apart from providing better initial guesses for trajectory design in fullephemeris model, our analytical derivations will be used for onboard spacecraft operations and orbit maintenance analyses of the MMX mission. Indeed, avoiding numerical integrations is crucial because of the limited computational resources of satellites' processors. New guidance algorithms based on mean relative orbit elements will be studied to deal with knowledge and dynamical errors in complex deep-space environments. Unfortunately, questions remain whether similar derivations can be carried out for low-altitude QSOs, which is why future research will be focused on extending the applicability of the averaged approach near the surface of the target bodies.
to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.