# On the role of the fast oscillations in the secular dynamics of the lunar
coplanar perturbation on Galileo satellites

Elisa Maria Alessi  IMATI-CNR, Istituto di Matematica Applicata e Tecnologie
informatiche “E. Magenes”, Consiglio Nazionale delle Ricerche, Via Alfonso
Corti 12, 20133 Milano, Italy  Inmaculada Baldomá  Departament de Matemàtiques
& IMTECH, Universitat Politècnica de Catalunya, Diagonal 647, 08028 Barcelona,
Spain  Centre de Recerca Matemàtica, Campus de Bellaterra, Edifici C, 08193
Barcelona, Spain  Mar Giralt  IMCCE, CNRS, Observatoire de Paris, Université
PSL, Sorbonne Université, 77 Avenue Denfert-Rochereau, 75014 Paris, France
Marcel Guardia  Departament de Matemàtiques i Informàtica, Universitat de
Barcelona, Gran Via, 585, 08007 Barcelona, Spain  Centre de Recerca
Matemàtica, Campus de Bellaterra, Edifici C, 08193 Barcelona, Spain  Alexandre
Pousse  IMATI-CNR, Istituto di Matematica Applicata e Tecnologie informatiche
“E. Magenes”, Consiglio Nazionale delle Ricerche, Via Alfonso Corti 12, 20133
Milano, Italy

###### Abstract

Motivated by the practical interest in the third-body perturbation as a
natural cleaning mechanism for high-altitude Earth orbits, we investigate the
dynamics stemming from the secular Hamiltonian associated with the lunar
perturbation, assuming that the Moon lies on the ecliptic plane. The secular
Hamiltonian defined in that way is characterized by two timescales. We compare
the location and stability of the fixed points associated with the secular
Hamiltonian averaged with respect to the fast variable with the corresponding
periodic orbits of the full system. Focusing on the orbit of the Galileo
satellites, it turns out that the two dynamics cannot be confused, as the
relative difference depends on the ratio between the semi-major axis of
Galileo and the one of the Moon, that is not negligible. The result is
relevant to construct rigorously the Arnold diffusion mechanism that can drive
a natural growth in eccentricity that allows a satellite initially on a
circular orbit in Medium Earth Orbit to reenter into the Earth’s atmosphere.

00footnotetext: _E-mail adresses:_<EMAIL_ADDRESS>(E.M. Alessi),
<EMAIL_ADDRESS>(I. Baldomá<EMAIL_ADDRESS>(M. Giralt),
<EMAIL_ADDRESS>(M. Guardia).

###### Contents

  1. 1 Introduction
  2. 2 From the Secular to the Averaged Hamiltonians
    1. 2.1 A good system of coordinates
    2. 2.2 The hierarchy of models
  3. 3 Main results
    1. 3.1 The critical points of the averaged system
    2. 3.2 The coplanar Hamiltonian
  4. 4 The circular critical points of the $h$-averaged Hamiltonian
    1. 4.1 Rewriting the $h$-averaged system
    2. 4.2 Circular critical points: Proof of Theorem 3.2
    3. 4.3 Modulus of the eigenvalues: Proof of Proposition 3.3
  5. 5 Eccentric critical points of the $h$-averaged system
    1. 5.1 Numerical approach
    2. 5.2 Existence of eccentric critical points: Formulation of the problem
    3. 5.3 A preliminary existence result
    4. 5.4 The fixed points of the form $(0,y)$
    5. 5.5 The fixed points of the form $(\pi,y)$
    6. 5.6 Linearization around the critical points
  6. 6 Periodic orbits of the coplanar Hamiltonian
    1. 6.1 Existence of periodic orbits. Proof of Theorem 3.7
      1. 6.1.1 Proof of Corollary 3.10
    2. 6.2 Perturbative analysis.
      1. 6.2.1 Proof of Proposition 3.12
      2. 6.2.2 Proof of Theorem 3.15
      3. 6.2.3 Proof of Lemma 3.18
      4. 6.2.4 Proof of Theorem 3.19
  7. A Expression of the Hamiltonian

##  1 Introduction

The third-body gravitational perturbation on a bounded orbit at the Earth is
known to yield a long-term variation in eccentricity (e.g., [4]). In the past
years, a specific interest on this phenomenon has arisen, as a possible
mechanism to facilitate the disposal of artificial satellites at the end-of-
life and dilute the probability of collision (e.g., [12]). Special emphasis in
this regard is given to high-altitude orbits, in particular to Medium Earth
Orbits (MEO), where the Global Navigation Satellite Systems (i.e., GPS,
Galileo, Beidou and GLONASS) are located.

Considering an equatorial geocentric reference system, the eccentricity growth
has been investigated mostly associated with resonances that involve the
argument of pericenter and the longitude of ascending node of the orbit of the
satellite and the longitude of the ascending node of the orbit of the third
body, namely, Moon or Sun. Pillar works on the identification of these
resonances, the corresponding phase space and their practical role on high-
altitude orbits are [17, 5, 10, 3, 22, 7]. Several numerical studies (e.g.,
[2, 20, 18]) have highlighted that indeed there exist regions in the phase
space such that a circular orbit can become very eccentric in the limit to
reenter to the Earth’s atmosphere without the need of any propulsion system.
However, from the theoretical point of view, a firm conclusion on how to
explain analytically the numerical results is still missing.

Motivated by this need, in [1] we show how an Arnold diffusion mechanism can
be built to explain the phenomenon. That work lies in the same context of a
series of works that are based on the analysis of the Kaula’s expansion of the
third-body perturbation [14], that is averaged over the mean anomaly of the
orbit of the spacecraft and the one of the orbit of the third body in order to
bring out the secular dynamics. We can say that these studies follow two main
different, but related, concepts: 1.) the eccentricity growth is due to a
chaotic behavior because more than one resonance is important in a given
region [21, 7]; 2.) the growth is due to a single resonance, whose phase space
is modified by other terms of the Kaula’s expansion [9, 6, 15].

In [1] we follow the latter, by considering the full quadrupolar secular
expansion of the lunar perturbation. The aim is to show how to construct an
Arnold diffusion mechanism, starting from the analysis of an autonomous
reduced model. That is, we analyze first what we call the “coplanar
Hamiltonian” (namely, assuming that the Moon lies in the ecliptic plane) and
use it as a first order for the full model. Then, we analyze the secular
Hamiltonian perturbatively taking the inclination of the Moon as a small
parameter. Note that Arnold diffusion can only take place for Hamiltonians
systems of at least 3 degrees of freedom and therefore the coplanar
Hamiltonian has too low dimension. However, it is a good first order to detect
the hyperbolic structures that enhance the drift of eccentricity for the full
secular Hamiltonian.

The main difference between our work and previous ones is that the first order
that we consider does not neglect the fast oscillations of the coplanar
dynamics. Let us explain what we mean by that.

As already mentioned, the mechanism considered in [1] to achieve eccentricity
growth is to drift along a resonance. More concretely we consider the $2:1$
resonance between the argument of the pericenter and the longitude of
ascending node of the orbit of the satellite. Then, the model has two
timescales: the slow one of the resonant angle and the fast one which is given
by the oscillations of the argument of the pericenter or of the longitude of
the ascending node.

One can consider as an effective model the coplanar Hamiltonian averaged over
the fast oscillations along the resonance, which is integrable. This averaged
model has been widely analyzed in literature (see [6, 15, 16]). It already
presents hyperbolic orbits and invariant manifolds which are the “highways” to
achieve eccentricity growth.

The purpose of the present paper is to analyze “how good” this averaged model
is. That is how the dynamics, and in particular, the hyperbolic invariant
objects in phase space of the coplanar Hamiltonian deviate from those of the
averaged coplanar model. In other words, how the fast oscillations in the
coplanar Hamiltonian affect its global dynamics.

More concretely:

  1. 1.

We analyze rigorously the dynamics of the averaged Hamiltonian. It turns out
that one can describe analytically in a rather precise way the location and
character of the critical points for a rather large range of semi-major axis
(which includes that of Galileo). Of major importance for [1] are the critical
points which lie at circular motion.

  2. 2.

We explain the need of overcoming the averaged approximation, that, as far as
we understood, lay the foundations of the works [6, 15, 16]. The main result
that we will present in this direction is to show that, at a fixed energy
level, the distance between hyperbolic periodic circular orbits of the
averaged coplanar system and that of the full coplanar Hamiltonian is
proportional to $a^{3}/a_{\mathrm{M}}^{3}$ with $a$ the semi-major axis of the
orbit and $a_{M}$ the semi-major axis of the Moon. This is done through a
perturbative analysis, that is, assuming $a^{3}/a_{\mathrm{M}}^{3}$ small
enough. Note, however, that for the semimajor axis of Galileo
$a^{3}/a_{\mathrm{M}}^{3}$ is not so small. That is, for this values one
should expect that the dynamics of the averaged Hamiltonian is rather “far”
from that of the full model even at short timescales. In particular, the
inclination of the periodic orbits presents a significant difference.

  3. 3.

We analyze the eigenvalues of the circular periodic orbits of the averaged
Hamiltonian and show that they are hyperbolic in a small interval, parabolic
at the endpoints of this interval and elliptic otherwise. We perform the same
analysis for the coplanar model by carrying out a perturbative analysis in
$a^{3}/a_{\mathrm{M}}^{3}$ of the eigenvalues of the periodic orbits.

The development presented is focused on the orbit of the Galileo
constellation, but the results are general and can be applied to other cases.

The paper is structured as follows. In Section 2 we introduce the models that
we consider and a suitable set of coordinates to analyze it. We follow closely
the exposition of the companion paper [1]. In Section 3, we state the main
results of the paper. In Section 3.1 those concerning the existence and
character of the critical points of the averaged system. Then, in Section 3.2,
we present the results concerning the circular periodic orbits of the coplanar
Hamiltonian emanating from the critical points of the averaged model. Sections
4 and 5 are devoted to prove, respectively, the results concerning to the
circular and non-circular (eccentric) equilibrium points of the averaged
system. Finally, in Section 6 we prove the results related with the periodic
orbits for the coplanar Hamiltonian system.

## Acknowledgements

E.M. Alessi and A. Pousse acknowledge support received by the project entitled
“co-orbital motion and three-body regimes in the solar system”, funded by
Fondazione Cariplo through the program “Promozione dell’attrattività e
competitività dei ricercatori su strumenti dell’European Research Council –
Sottomisura rafforzamento”.

I. Baldomá has been supported by the grant PID-2021- 122954NB-100 funded by
the Spanish State Research Agency through the programs
MCIN/AEI/10.13039/501100011033 and “ERDF A way of making Europe”.

M. Giralt has been supported by the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement No
101034255. M. Giralt has also been supported by the research project PRIN
2020XB3EFL “Hamiltonian and dispersive PDEs”.

M. Guardia has been supported by the European Research Council (ERC) under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No. 757802). M. Guardia is also supported by the Catalan Institution
for Research and Advanced Studies via an ICREA Academia Prize 2019. This work
is part of the grant PID-2021-122954NB-100 funded by
MCIN/AEI/10.13039/501100011033 and “ERDF A way of making Europe”.

This work is also supported by the Spanish State Research Agency, through the
Severo Ochoa and María de Maeztu Program for Centers and Units of Excellence
in R&D (CEX2020-001084-M).

##  2 From the Secular to the Averaged Hamiltonians

Let us consider a spacecraft that moves under the gravitational attraction of
the Earth, the perturbation due to the Earth’s oblateness and the lunar
gravitational perturbation. Assuming a geocentric equatorial reference system,
the orbital elements, namely, semi-major axis $a$, eccentricity $e$,
inclination $i$, longitude of the ascending node $\Omega$ and argument of
pericenter $\omega$, that define the ellipse where the motion occurs, change
in time because of the perturbations. The orbit of the Moon is defined in the
geocentric ecliptic reference system by the corresponding orbital elements
$(a_{\mathrm{M}},e_{\mathrm{M}},i_{\mathrm{M}},\Omega_{\mathrm{M}},\omega_{\mathrm{M}})$,
where

$a_{\mathrm{M}}=384400\,\mbox{km},\qquad e_{\mathrm{M}}=0.0549006,\qquad
i_{\mathrm{M}}=5.15^{\circ},$

while the longitude of the ascending node of the Moon with respect to the
ecliptic plane varies approximately linearly with time in a period
$T_{\Omega_{\mathrm{M}}}$ of 1 Saros (about 6585.321347 days [23, 19]) due to
the solar gravitational perturbation, namely,

$\Omega_{\mathrm{M}}(t)=\Omega_{{\mathrm{M}},0}+n_{\Omega_{\mathrm{M}}}t,\quad\quad\quad
n_{\Omega_{\mathrm{M}}}=2\pi/T_{\Omega_{\mathrm{M}}}.$ (2.1)

In Delaunay action-angle variables $(L,G,H,\ell,g,h)$, see, for instance,
[24], considering as small parameter of the problem

$\alpha=\frac{a}{a_{\mathrm{M}}},$

which characterizes the distance of the Moon with respect to the orbit of the
satellite, the secular dynamics, see [7], is described by the non-autonomous
(recall (2.1)) Hamiltonian

$\mathtt{H}(L,G,H,g,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})={\mathrm{H}}_{\mathrm{K}}(L)+\tilde{{\mathrm{H}}}_{0}(L,G,H)+\alpha^{3}\tilde{{\mathrm{H}}}_{1}(L,G,H,g,h,\Omega_{\mathrm{M}};i_{\mathrm{M}}),$
(2.2)

where each term is defined as follows.

The first term,

${\mathrm{H}}_{\mathrm{K}}(L)=-\frac{1}{2}\frac{\mu^{2}}{L^{2}},$

is the constant term associated with the Earth’s monopole, being

$\mu=398600.44\,\mbox{km}^{3}/\mbox{s}^{2}$ (2.3)

the mass parameter of the Earth.

The second term,

$\tilde{{\mathrm{H}}}_{0}(L,G,H)=\frac{1}{4}\frac{\rho_{0}}{L^{3}}\frac{G^{2}-3H^{2}}{G^{5}},$
(2.4)

is the perturbative term associated with the Earth’s oblateness, averaged over
the orbital period of the spacecraft, being $\rho_{0}=\mu^{4}J_{2}R^{2}$, with
$J_{2}=1.08\times 10^{-3}$ the coefficient of the second zonal harmonic in the
geopotential and $R=6378.14\,\mbox{km}$ the mean equatorial radius of the
Earth.

Finally, the (truncated) secular perturbative term due to the Moon (see [13,
11]) corresponds to

$\begin{split}\tilde{{\mathrm{H}}}_{1}&(L,G,H,g,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})\\\
=&-\frac{\rho_{1}}{L^{2}}\sum_{m=0}^{2}\sum_{p=0}^{2}{\tilde{D}}_{m,p}(L,G,H)\sum_{s=0}^{2}c_{m,s}F_{2,s,1}(i_{\mathrm{M}})\\\
&\times\left[U_{2}^{m,-s}(\epsilon)\cos\left(\tilde{\psi}_{m,p,s}(g,h,\Omega_{\mathrm{M}})\right)+U_{2}^{m,s}(\epsilon)\cos\left(\tilde{\psi}_{m,p,-s}(g,h,\Omega_{\mathrm{M}})\right)\right].\end{split}$
(2.5)

In the expression above,

  * •

the function $U_{2}^{m,\mp s}(\epsilon)$ corresponds to the Giacaglia function
(see Tab. 4 in Appendix A) being $\epsilon=23.44^{\circ}$ the obliquity of the
ecliptic with respect to the equatorial plane.

  * •

One has that

${\tilde{D}}_{m,p}(L,G,H)=\tilde{F}_{m,p}(G,H)\tilde{X}_{p}(L,G),$ (2.6)

with

$\begin{split}\tilde{F}_{m,p}(G,H)=F_{2,m,p}\left(\arccos\frac{H}{G}\right),\quad\tilde{X}_{p}(L,G)=X_{0}^{2,2-2p}\left(\sqrt{1-\frac{G^{2}}{L^{2}}}\right),\end{split}$

where $F_{2,m,p}(i)$ are the Kaula’s inclination functions (see [13]),
$X_{0}^{2,2-2p}(e)$ are the zero-order Hansen coefficients (see [11]) and we
have used the definitions of inclination and eccentricity in terms of Delaunay
coordinates, that is

$H=G\cos i\qquad\text{and}\qquad G=L\sqrt{1-e^{2}}.$ (2.7)

  * •

The angle

$\tilde{\psi}_{m,p,s}(g,h,\Omega_{\mathrm{M}})=2(1-p)g+mh+s\left(\Omega_{\mathrm{M}}-\frac{\pi}{2}\right)-{\mathrm{y}}_{\left|s\right|}\pi$
(2.8)

defines the relative orientation satellite-Moon.

  * •

The other constants are defined as

$\displaystyle\rho_{1}$ $\displaystyle=$
$\displaystyle\frac{\mu\mu_{\mathrm{M}}}{(1-e_{\mathrm{M}}^{2})^{3/2}},$
$\displaystyle\varepsilon_{n}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}1&\mbox{if $n=0$}\\\ 2&\mbox{if
$n\neq 0$}\end{array}\right.,$ (2.9) $\displaystyle c_{m,s}$ $\displaystyle=$
$\displaystyle(-1)^{\lfloor
m/2\rfloor}\frac{\varepsilon_{m}\varepsilon_{s}}{2}\frac{(2-s)!}{(2+m)!},{}$
$\displaystyle{\mathrm{y}}_{s}$ $\displaystyle=$
$\displaystyle\left\\{\begin{array}[]{ll}0&\mbox{if $s$ is even}\\\
1/2&\mbox{if $s$ is odd}\end{array}\right.,$

where $\mu_{\mathrm{M}}=4902.87\,\mbox{km}^{3}/\mbox{s}^{2}$ is the mass
parameter of the Moon.

######  Remark 2.1.

By Kaula’s inclination functions (see [13]), one obtains that

$\displaystyle F_{2,s,1}(0)=\left\\{\begin{array}[]{ll}-\frac{1}{2}&\mbox{if
$s=0$,}\\\\[5.0pt] 0&\mbox{if $s=1,2$.}\end{array}\right.$

That is, assuming that the Moon lies on the ecliptic plane (i.e.,
$i_{\mathrm{M}}=0$), the Hamiltonian $\mathtt{H}$ in (2.2) is autonomous,
because the angle $\tilde{\psi}_{m,p,s}(g,h,\Omega_{\mathrm{M}})$ does not
depend on $\Omega_{\mathrm{M}}$ for $s=0$ (recall (2.1)).

######  Remark 2.2.

In the numerical computations that will be presented throughout the text, all
the variables will be taken in non-dimensional units defined in such a way
that the semi-major axis (equal to 29600 km, because we focus on Galileo) is
the unit of distance and the corresponding orbital period is 2$\pi$.

In what follows, we will omit the Keplerian term of the Hamiltonian
${\mathrm{H}}_{{\mathrm{K}}}$, since it does not contribute to the variation
of the orbital elements. Similarly, we omit the dependence of the Hamiltonian
on the variable $L$, since it is a constant of motion.

###  2.1 A good system of coordinates

The Galileo constellation orbits at an inclination, $i\approx 56^{\circ}$,
such that the dominant term in the lunar perturbation is the $2g+h$ resonance.
This resonance is defined by the unperturbed Hamiltonian
$\tilde{{\mathrm{H}}}_{0}$ in (2.4) (i.e., $\alpha=0$) and the condition
$2\dot{g}+\dot{h}=0$. Namely,
$2\partial_{g}\tilde{{\mathrm{H}}}_{0}+\partial_{h}\tilde{{\mathrm{H}}}_{0}=0$.
This is satisfied provided that

$5H^{2}-G^{2}-HG=0,$

that is, for all $G\neq 0$ (i.e., $e<1$) and $H=G\cos i_{\star}$ that satisfy

$5X^{2}-1-X=0,\quad\quad\quad X=\cos i_{\star}.$

Hence, we can distinguish two situations:

  * •

the prograde case for
$i_{\star}=\arccos\left(\frac{1+\sqrt{21}}{10}\right)\simeq 56,06^{\circ}$,

  * •

the retrograde case for
$i_{\star}=\arccos\left(\frac{1-\sqrt{21}}{10}\right)\simeq 110,99^{\circ}$.

In what follows, $i_{\star}$ will be mentioned as the inclination of the
“exact $2g+h$”-resonance and we will focus on the prograde case.

The unperturbed Hamiltonian highlights that $x=2g+h$ is constant when
$i=i_{\star}$, while it circulates for $i\neq i_{\star}$. Moreover, in a small
enough neighborhood of the exact resonance the angular variables evolve at
different rates: $g$ and $h$ are “fast” angles compared to $x$ which undergoes
a “slow” drift of order ${\mathcal{O}}(i-i_{\star})$.

##### Slow-fast Delaunay coordinates $(y,x)$.

Instead of using the Delaunay action-angle variables, in order to take
advantage of the timescales separation, we introduce the symplectic
transformation ${(G,H,g,h)=\Upsilon^{\mathrm{Del}}(y,\Gamma,x,h)}$ where

$y=\frac{G}{2},\qquad\qquad\Gamma=H-\frac{G}{2},\qquad\qquad x=2g+h.$ (2.10)

Note that, the resonant action $y$ does not depend on the inclination. Hence,
the action-angle variables $(y,x)$ are associated with the variation in
eccentricity.

In slow-fast Delaunay variables, the secular Hamiltonian can be written as

${\mathrm{H}}(y,\Gamma,x,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})={\mathrm{H}}_{0}(y,\Gamma)+\alpha^{3}{\mathrm{H}}_{1}(y,\Gamma,x,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})$

where

$\begin{split}{\mathrm{H}}_{0}(y,\Gamma)&=(\tilde{{\mathrm{H}}}_{0}\circ\Upsilon^{\mathrm{Del}})(y,\Gamma)=\frac{\rho_{0}}{128}\frac{y^{2}-6y\Gamma-3\Gamma^{2}}{L^{3}y^{5}}\\\
{\mathrm{H}}_{1}(y,\Gamma,x,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})&=(\tilde{{\mathrm{H}}}_{1}\circ\Upsilon^{\mathrm{Del}})(y,\Gamma,x,h,\Omega_{\mathrm{M}};i_{\mathrm{M}}).\end{split}$
(2.11)

In the coordinates just defined, the $2g+h$–resonance becomes the
$x$–resonance, which is defined by

$\omega(y,\Gamma)=\frac{\partial{\mathrm{H}}_{0}(y,\Gamma)}{\partial
y}=\frac{3}{128}\frac{\rho_{0}}{L^{3}}\frac{5\Gamma^{2}+8y\Gamma-y^{2}}{y^{6}}=0.$

This resonance takes place at the two lines

$\Gamma=y\frac{-4+\sqrt{21}}{5}\qquad\text{and}\qquad\Gamma=y\frac{-4-\sqrt{21}}{5},$
(2.12)

with $(y,\Gamma)\neq(0,0)$, that are associated with prograde and retrograde
orbits, respectively.

##### Poincaré coordinates $(\eta,\xi)$.

The slow-fast variables (2.10) derived from the Delaunay variables (as happens
with the original Delaunay variables) are singular at $e=0$. Thus, in order to
study the dynamics of the secular Hamiltonian in a neighborhood of circular
orbits, i.e., $0<e\ll 1$, we introduce the Poincaré coordinates
${(y,\Gamma,x,h)=\Upsilon^{\mathrm{Poi}}(\eta,\Gamma,\xi,h)}$ where

$\xi=\sqrt{2L-4y}\cos\left(\frac{x}{2}\right),\qquad\qquad\qquad\eta=\sqrt{2L-4y}\sin\left(\frac{x}{2}\right),$

which are symplectic. Notice that $\xi$ and $\eta$ are respectively equivalent
to $e\cos(x/2)$ and $e\sin(x/2)$ for quasi-circular orbits.

In these coordinates, the secular Hamiltonian becomes

$\mathcal{H}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})=\mathcal{H}_{0}(\eta,\Gamma,\xi)+\alpha^{3}\mathcal{H}_{1}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})$
(2.13)

where

$\displaystyle\mathcal{H}_{0}(\eta,\Gamma,\xi)$
$\displaystyle=({\mathrm{H}}_{0}\circ\Upsilon^{\mathrm{Poi}})(\eta,\Gamma,\xi)$
(2.14)
$\displaystyle=\frac{\rho_{0}}{2}\frac{(2L-\xi^{2}-\eta^{2})^{2}-24(2L-\xi^{2}-\eta^{2})\Gamma-48\Gamma^{2}}{L^{3}(2L-\xi^{2}-\eta^{2})^{5}}$
$\displaystyle\mathcal{H}_{1}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})$
$\displaystyle=({\mathrm{H}}_{1}\circ\Upsilon^{\mathrm{Poi}})(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}}).$
(2.15)

###  2.2 The hierarchy of models

In the present paper, we consider a hierarchy of models which stems from the
secular Hamiltonian (2.13). They are what we call the _coplanar secular
Hamiltonian_ and the _$h$ -averaged (coplanar secular) Hamiltonian_. To
construct these “intermediate models”, we rely on the fact that the model
depends on two parameters: the inclination of the Moon with respect to the
ecliptic plane, $i_{\mathrm{M}}$, and the semi-major axis ratio between the
satellite and the Moon, $\alpha$.

In the companion paper [1], using the same hierarchy of models, we construct
Arnold diffusion orbits for the secular Hamiltonian (2.13), that can lead to a
drastic increase of the eccentricity of the satellite. A crucial point in our
construction is the analysis of certain hyperbolic structures (normally
hyperbolic invariant manifolds, stable and unstable invariant manifolds) of
the coplanar secular Hamiltonian which are persistent in (2.13). Comparing the
coplanar model with the full one allows us to construct the drifting orbits.

On the contrary, in the present paper we focus on the comparison between the
coplanar and the $h$-averaged Hamiltonians. There is now an extensive
literature (see, e.g., [6, 16]) analyzing the hyperbolic critical points of
the $h$-averaged Hamiltonian and its invariant manifolds as a sign of
existence of unstable motions for the full Hamiltonian (2.13). In the present
paper, we compare analytically certain dynamics of the $h$-averaged and
coplanar models focusing both on what dynamics persist and which “deviations”
has one model with respect to the other.

##### First reduction: the Coplanar Model.

Since the inclination of the Moon with respect to the ecliptic plane is
relatively small (that is, $i_{\mathrm{M}}=5.15^{\circ}$), the first reduction
that one can do to have an intermediate model is to take $i_{\mathrm{M}}=0$,
which corresponds to assume that the orbit of the Moon is coplanar to that of
the Earth.

When doing this reduction, the Hamiltonian $\mathcal{H}$ in (2.13) becomes
$\Omega_{\mathrm{M}}$ independent and thus autonomous (see Remark 2.1).
Starting from (2.13), we define the $2$ degrees of freedom Hamiltonian

$\mathcal{H}_{\mathrm{CP}}(\eta,\Gamma,\xi,h)=\mathcal{H}(\eta,\Gamma,\xi,h,\Omega_{M};0),$
(2.16)

where the subscript $\mathrm{CP}$ stands for coplanar. Since this Hamiltonian
is autonomous, $\mathcal{H}_{\mathrm{CP}}$ is a first integral of the system.
Moreover, it can be written as

$\mathcal{H}_{\mathrm{CP}}(\eta,\Gamma,\xi,h)=\mathcal{H}_{0}(\eta,\Gamma,\xi)+\alpha^{3}\mathcal{H}_{\mathrm{CP},1}(\eta,\Gamma,\xi,h),$

where $\mathcal{H}_{0}$ is the Hamiltonian introduced in (2.14) and
$\mathcal{H}_{\mathrm{CP},1}$ is the Hamiltonian $\mathcal{H}_{1}$ in (2.15)
with $i_{\mathrm{M}}=0$, that is,

$\mathcal{H}_{\mathrm{CP},1}(\eta,\Gamma,\xi,h)=\mathcal{H}_{1}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};0).$
(2.17)

See Appendix A, and in particular (A.1), for the explicit expression of
$\mathcal{H}_{\mathrm{CP},1}$.

##### Second reduction: the $h$-averaged problem.

When $\alpha$ is also taken as a small parameter, the autonomous Hamiltonian
$\mathcal{H}_{\mathrm{CP}}$ in (2.16) has a timescale separation between the
slow and fast angles. Indeed,

$\dot{\eta},\dot{\xi}\sim\alpha^{3}\qquad\text{whereas}\qquad\dot{h}\sim 1.$

A classical way to exploit this feature is to simplify the Hamiltonian
$\mathcal{H}_{\mathrm{CP}}$ by averaging out the fast oscillations with
respect to the longitude of the ascending node $h$. That is,

$\mathcal{H}_{\mathrm{AV}}(\eta,\Gamma,\xi)=\frac{1}{2\pi}\int_{0}^{2\pi}\mathcal{H}_{\mathrm{CP}}(\eta,\Gamma,\xi,h){\mathrm{d}}h.$

Note that, since $\mathcal{H}_{0}$ in (2.14) is $h$-independent,
$\mathcal{H}_{\mathrm{AV}}$ can be written as

$\mathcal{H}_{\mathrm{AV}}(\eta,\Gamma,\xi)=\mathcal{H}_{0}(\eta,\Gamma,\xi)+\alpha^{3}\mathcal{H}_{\mathrm{AV},1}(\eta,\Gamma,\xi),$
(2.18)

with

$\mathcal{H}_{\mathrm{AV},1}(\eta,\Gamma,\xi)=\frac{1}{2\pi}\int_{0}^{2\pi}\mathcal{H}_{\mathrm{CP},1}(\eta,\Gamma,\xi,h){\mathrm{d}}h.$

More specifically, the particular form of $\mathcal{H}_{\mathrm{CP},1}$ in
(A.1) implies that

$\begin{split}\mathcal{H}_{\mathrm{AV},1}(\eta,\Gamma,\xi)=\frac{\rho_{1}}{L^{2}}\bigg{[}&\frac{1}{2}U_{2}^{0,0}\mathcal{D}_{0,1}(\eta,\Gamma,\xi)\left(8L^{2}+12L(\xi^{2}+\eta^{2})-3(\xi^{2}+\eta^{2})^{2}\right)\\\
&+\frac{1}{6}U_{2}^{1,0}\mathcal{D}_{1,0}(\eta,\Gamma,\xi)(\xi^{2}-\eta^{2})\bigg{]},\end{split}$
(2.19)

where the constants $U_{2}^{m,0}$ and the functions $\mathcal{D}_{m,p}$ are
given in Tables 4 and 6, respectively, in Appendix A. See also this appendix
for the deduction of this expression.

To analyze the dynamics of the $h$-averaged Hamiltonian, it will be also
convenient to consider it in the slow-fast Delaunay coordinates (2.10). In
these variables, it reads

${\mathrm{H}}_{\mathrm{AV}}(y,\Gamma,x)={\mathrm{H}}_{0}(y,\Gamma)+{\mathrm{H}}_{\mathrm{AV},1}(y,\Gamma,x)$
(2.20)

where ${\mathrm{H}}_{0}$ is the Hamiltonian introduced in (2.11),

${\mathrm{H}}_{\mathrm{AV},1}(y,\Gamma,x)=\frac{\rho_{1}}{L^{2}}\left[\frac{1}{2}U_{2}^{0,0}D_{0,1}(y,\Gamma)+\frac{1}{3}U_{2}^{1,0}D_{1,0}(y,\Gamma)\cos
x\right],$

and the functions $D_{m,p}$ are given in Table 5. We also refer to Appendix A,
for the deduction of the expression of ${\mathrm{H}}_{\mathrm{AV},1}$.

##  3 Main results

In this section, we present the main results of the paper, which deal with
both the $h$-averaged Hamiltonian and the coplanar Hamiltonian. They are
stated respectively in Sections 3.1 and Section 3.2.

######  Remark 3.1.

Along the paper we also present some numerical results. As a rule, the
variables will be taken in non-dimensional units defined in such a way that
the semi-major axis of the satellite (equal to 29600 km) is the unit of
distance and the corresponding orbital period is 2$\pi$.

However, in some places, with the purpose of keeping the exposition clear, we
have kept the standard unit system. It will be pointed out along the text.

The results for the $h$-averaged Hamiltonian are related with the existence
and stability character of the stationary points. To be more precise,

  * •

We study the stability character of the origin of the $h$-averaged Hamiltonian
in Theorem 3.2. Then, in Proposition 3.3, we provide some properties of the
eigenvalues of the linearization of $\mathcal{H}_{\mathrm{AV}}$ at the origin
(in Poincaré variables).

  * •

Theorem 3.5 deals with the existence and the linear stability behavior of the
eccentric critical points of the $h$-averaged Hamiltonian
${\mathrm{H}}_{\mathrm{AV}}$. We state these results in Delaunay coordinates.

These results are proven for a specific range of the semi-major axis $a$,
namely, we prove the results for values of $a\in[a_{\min},a_{\max}]$ with

$a_{\min}=6378.14\text{ km}\qquad\text{and}\qquad a_{\max}=30000\text{ km}.$
(3.1)

The minimum value corresponds to a collision orbit with the Earth assuming
$e=0$. Note also that this interval contains the semi-major axis of Galileo
$a=29600\text{ km}$.

The coplanar Hamiltonian $\mathcal{H}_{\mathrm{CP}}$ can be seen as a
perturbation of the $h$-averaged Hamiltonian $\mathcal{H}_{\mathrm{AV}}$. From
this point of view and recalling that the origin $\xi=\eta=0$ is a fixed point
of $\mathcal{H}_{\mathrm{AV}}$, in Section 3.2, we state the results for the
coplanar Hamiltonian $\mathcal{H}_{\mathrm{CP}}$ related to the periodic
orbits emanating from the origin. Some of these results are perturbative
results with respect to the parameter $\delta$ defined as

$\delta=\rho\alpha^{3}L^{4},\qquad\mathrm{with}\qquad\alpha=\frac{a}{a_{{\mathrm{M}}}},\qquad\rho=\frac{\rho_{1}}{\rho_{0}}$
(3.2)

with $\alpha,\rho_{0},\rho_{1}$ introduced in Section 2.

Next, we study the existence and main properties of the periodic orbits of the
coplanar Hamiltonian:

  * •

In Theorem 3.7, we characterize the energy levels of the coplanar Hamiltonian
$\mathcal{H}_{\mathrm{CP}}$ containing a periodic orbit lying in
$\\{\xi=\eta=0\\}$.

  * •

For these values of energy, $\mathbf{E}$, Corollary 3.10 relates the $\Gamma$
value of the point of the periodic orbit at the energy level
$\\{\mathcal{H}_{\mathrm{CP}}=\mathbf{E}\\}$ at the transverse section
$\\{h=0\\}$ with the corresponding parameter $\Gamma$ such that the
$h$-averaged Hamiltonian has the same energy level, namely
$\\{\mathcal{H}_{\mathrm{AV}}=\mathbf{E}\\}$. We see that, for small enough
$\delta$, the two values are order $\delta$ separated. Since for Galileo
$\delta\sim 0.11$ (computed in the unities specified in Remark 3.1), this
shows that $h$-averaged model differs significantly from the coplanar one. At
circular motions, deviations in $\Gamma$ are equivalent to deviations in (the
cosine of the) inclination.

  * •

Proposition 3.12 computes a first order approximation with respect to $\delta$
of the periodic orbits and their corresponding periods.

Even when there is no physical meaning for values of $a<a_{\min}$, these
results hold true for values of $a\in(0,a_{\max}]$.

The analysis of the stability character of the periodic orbits is the only one
that needs the parameter $\delta$ to be small enough.

  * •

Theorem 3.15, and its refinement Theorem 3.19, give an approximation, at least
of order $\mathcal{O}(\delta^{3/2})$, of the eigenvalues of the monodromy
matrix associated to the periodic orbits by means of the eigenvalues of the
linearized part at the origin for the $h$-averaged Hamiltonian. As a
consequence we are able to elucidate, for small values of $\delta$, the
character of the periodic orbits.

###  3.1 The critical points of the averaged system

We now present the results about the critical points of the $h$-averaged
Hamiltonian. We notice that, to simplify the exposition, we state some results
referred to the $h$-averaged Hamiltonian in Poincaré variables (see (2.18)),
and others referred to the Hamiltonian expressed in slow-fast Delaunay
variables (see (2.20)).

The specific range of the semi-major axis $[a_{\min},a_{\max}]$ is given in
(3.1). We also introduce the notation

$L_{\min}=\sqrt{\mu a_{\min}},\qquad L_{\max}=\sqrt{\mu
a_{\max}},\qquad\alpha_{\max}=\frac{a_{\max}}{a_{\mathrm{M}}},\qquad\delta_{\max}=\rho\alpha_{\max}^{3}L_{\max}^{4},$
(3.3)

where $\mu$ has been introduced in (2.3). Finally, we define

$m=\frac{-4+\sqrt{21}}{5}\sim 0.116515138991168,$ (3.4)

the slope of the resonance prograde line is $\Gamma=my$ (see (2.12)).

In addition, we notice that, since we want to study the prograde resonance we
keep our analysis to $\Gamma>0$ and, from the expressions of
$\mathcal{D}_{m,p}$ in Table 6 , we deduce that

$0\leq\Gamma\leq\frac{L}{2}.$

The first result deals with the character of the origin in Hamiltonian (2.18).
Its proof is carried out in Section 4.2.

######  Theorem 3.2.

The origin is a critical point of the $h$-averaged Hamiltonian
$\mathcal{H}_{\mathrm{AV}}$, see (2.18). Moreover, there exist two functions
${\Gamma}_{-},{\Gamma}_{+}:[L_{\min},L_{\max}]\to\left(0,\frac{L}{2}\right)$,
such that the character of the critical point is described as follows for
$L\in[L_{\min},L_{\max}]$:

  * •

if either $\Gamma\in\big{(}0,\Gamma_{-}(L)\big{)}$ or
$\Gamma\in\left(\Gamma_{+}(L),\frac{L}{2}\right)$, then $(0,0)$ is a center.

  * •

If $\Gamma\in\big{(}\Gamma_{-}(L),\Gamma_{+}(L)\big{)}$, then $(0,0)$ is a
saddle.

  * •

If $\Gamma=\Gamma_{-}(L)$ or $\Gamma=\Gamma_{+}(L)$, the origin is a
degenerated equilibrium point with nilpotent linearization.

In addition, $L^{-1}\Gamma_{\pm}(L)$ can be expressed as

$\frac{\Gamma_{\pm}(L)}{L}=\hat{\Gamma}_{\pm}(\delta),\qquad\delta=\rho\alpha^{3}L^{4},$

where $\hat{\Gamma}_{\pm}:[0,\delta_{\max}]\to\mathbb{R}$, see (3.2) and
(3.3), and there exists a constant $C_{*}$ such that for all
$L\in[L_{\min},L_{\max}]$,

$\left|\hat{\Gamma}_{\pm}(\delta)-\frac{m}{2}\mp\frac{5U_{2}^{1,0}(m+3)\sqrt{3-2m-4m^{2}}}{{3(16+20m)}}\delta\right|\leq
C_{*}\delta^{2}.$ (3.5)

We denote by $\pm\lambda_{\mathrm{AV}}(\Gamma;L)$ the eigenvalues of the
linearization around the origin of the Hamiltonian system associated to
(2.18). We assume that, in the saddle case,
$\lambda_{\mathrm{AV}}(\Gamma;L)>0$ and in the center case,
$-i\lambda_{\mathrm{AV}}(\Gamma;L)>0$. Notice that
$\lambda_{\mathrm{AV}}(\Gamma_{-}(L);L)=\lambda_{\mathrm{AV}}(\Gamma_{+}(L);L)=0$
where $\Gamma_{-}(L)$ and $\Gamma_{+}(L)$ are the functions defined in Theorem
3.2.

######  Proposition 3.3.

For any $L\in[L_{\min},L_{\max}]$ we have the following.

  * •

There exists a unique
$\Gamma_{*}(L)\in\big{(}\Gamma_{-}(L),\Gamma_{+}(L)\big{)}$ such that for all
$\Gamma\in\big{(}\Gamma_{-}(L),\Gamma_{+}(L)\big{)}$,
$\lambda_{\mathrm{AV}}(\Gamma;L)\leq\lambda_{\mathrm{AV}}(\Gamma_{*}(L);L)$
and

$\frac{5\rho_{0}}{16L^{7}}U_{2}^{1,0}\delta(m+3)\sqrt{3-2m-m^{2}}\leq\lambda_{\mathrm{AV}}(\Gamma_{*}(L);L)\leq\frac{15\sqrt{3}\rho_{0}}{16L^{7}}U_{2}^{1,0}\delta.$

  * •

For $\Gamma\in\big{(}0,\Gamma_{-}(L))$ the function
$|\lambda_{\mathrm{AV}}(\Gamma;L)|$ is smooth and decreasing with respect to
$\Gamma$ (and $L$ fixed).

  * •

For $\Gamma\in\left(\Gamma_{+}(L),\frac{L}{2}\right)$, the function
$|\lambda_{\mathrm{AV}}(\Gamma;L)|$ is smooth and increasing with respect to
$\Gamma$.

In addition, there exists a function
$\hat{\lambda}_{\mathrm{AV}}:\left[0,\frac{1}{2}\right]\times[0,\delta_{\max}]\to\mathbb{C}$
such that

$\lambda_{\mathrm{AV}}(\Gamma;L)=\frac{3\rho_{0}}{4L^{7}}\hat{\lambda}_{\mathrm{AV}}\left(\frac{\Gamma}{L};\delta\right),\qquad\delta=\rho\alpha^{3}L^{4}.$

######  Remark 3.4.

The eigenvalues $\lambda_{\mathrm{AV}}(\Gamma;L)$ can be explicitly computed,
see Remark 4.1.

This result is proven in Section 4.3. We also provide some numerical
computations.

To finish with the results related to the equilibrium points of the
$h$-averaged system, we state the following theorem concerning “non-circular”
critical points, namely, those critical points of the $h$-averaged system
(2.20) (in Delaunay coordinates) satisfying $e\neq 0$. We consider slow-fast
Delaunay variables and the Hamiltonian in (2.20).

Note that we are interested in a region of the phase space satisfying the
following conditions:

  * •

We want the satellite to be in a prograde “elliptic regime”. Equivalently, we
want the Delaunay coordinates to be well-defined and, therefore, we consider

$(x,y)\in\mathbb{T}\times(0,L/2).$ (3.6)

  * •

We want the actions $y$ and $\Gamma$ to be close to the resonance (see, for
instance, the first resonance (2.12)), namely the prograde resonance, so that
the $h$-averaged Hamiltonian is a good approximation of the coplanar model.
For this reason, we consider

$\frac{\Gamma}{y}\in\left(0,\frac{1}{2}\right).$ (3.7)

We recall that the resonance prograde line is $\Gamma=my$ with $m$ defined in
(3.4) and we notice that, from Table 5 the values of $\Gamma,y$ are restricted
to $\Gamma<y$.

  * •

To avoid collision, the semi-major axis $a$ and the eccentricity $e$ have to
satisfy $a(1-e)\geq R=a_{\min}$. Then $e\leq 1-\frac{R}{a}$. Moreover, since
$\frac{y}{L}=\frac{1}{2}\sqrt{1-e^{2}}$, we are only interested in values

$\frac{y}{L}\geq\frac{1}{2}\sqrt{2\frac{R}{a}-\frac{R^{2}}{a^{2}}}\geq\frac{1}{2}\sqrt{2\frac{R}{a_{\max}}-\frac{R^{2}}{a_{\max}^{2}}}=:\hat{y}_{\mathrm{col}}=0.308224183446436.$
(3.8)

In conclusion, we consider values of
$\frac{y}{L}\in\left[\hat{y}_{\mathrm{col}},\frac{1}{2}\right]$.

######  Theorem 3.5.

There exist two constants111The value $L_{*}$ (in international unities)
corresponds to a semi-major axis $a_{*}\sim 23893.56218133389\,\text{km}$.

$L_{*}\sim 97590.90325766560,\qquad\hat{\Gamma}_{0}\sim 0.013584391815073,$

two decreasing functions $\hat{\Gamma}_{1}(L)$, $\hat{\Gamma}_{2}(L)$ and one
increasing function $\hat{\Gamma}_{3}(L)$ defined for
$L\in[L_{\min},L_{\max}]$ such that the Hamiltonian (2.20) has the following
critical points satisfying (3.6), (3.7), (3.8):

  * •

For $L\in[L_{\min},L_{*}]$ and $\Gamma\in(0,L/2]$,

    1. 1.

if $\Gamma\notin[L\hat{\Gamma}_{1}(L),L\hat{\Gamma}_{3}(L)]$, has no critical
points.

    2. 2.

If $\Gamma\in[L\hat{\Gamma}_{1}(L),L\hat{\Gamma}_{0})$, has a unique critical
point that is a center and of the form $(\pi,y)$.

    3. 3.

If $\Gamma\in[L\hat{\Gamma}_{0},L\hat{\Gamma}_{2}(L)]$, has only two critical
points $(0,y_{1}),(\pi,y_{2})$, which are of saddle and center type
respectively, except when $L=L_{*}$ and $\Gamma=L_{*}\hat{\Gamma}_{2}(L_{*})$
that $(\pi,y_{2})$ is parabolic.

    4. 4.

If $\Gamma\in(L\hat{\Gamma}_{2}(L),L\hat{\Gamma}_{3}(L)]$, has a unique
critical point that is a saddle of the form $(0,y)$.

  * •

If $L\in(L_{*},L_{\max}]$ and $\Gamma\in(0,L/2]$, there exists a function
$\hat{\Gamma}_{*}(L)\in[\hat{\Gamma}_{2}(L),\hat{\Gamma}_{3}(L)]$ such that,

    1. 1.

if $\Gamma\notin[L\hat{\Gamma}_{1}(L),L\hat{\Gamma}_{3}(L)]$, has no critical
points.

    2. 2.

If $\Gamma\in[L\hat{\Gamma}_{1}(L),L\hat{\Gamma}_{0})$, has a unique critical
point that is a center and of the form $(\pi,y)$.

    3. 3.

If $\Gamma\in[L\hat{\Gamma}_{0},L\hat{\Gamma}_{2}(L))$, has only two critical
points $(0,y_{1}),(\pi,y_{2})$, of saddle and center type respectively.

    4. 4.

If $\Gamma\in[L\hat{\Gamma}_{2}(L),L\hat{\Gamma}_{*}(L))$, has only three
critical points $(0,y),(\pi,y_{1,2})$ where $y_{2}>y_{1}>0$ with $(0,y)$,
$(\pi,y_{2})$ saddles and $(\pi,y_{1})$ a center.

    5. 5.

If $\Gamma=L\hat{\Gamma}_{*}(L)$, has only two critical points $(0,y)$,
$(\pi,y_{*})$ with $(0,y)$ a saddle and $(\pi,y_{*})$ parabolic.

    6. 6.

If $\Gamma\in(L\hat{\Gamma}_{*}(L),L\hat{\Gamma}_{3}(L)]$, has a unique
critical point that is a saddle and of the form $(0,y)$.

The complete proof of this result relies on elementary analytic tools.
However, due to the complexity of the expression of the functions involved in
the definition of the Hamiltonian (2.18), it is not straightforward at all and
it is postponed to Section 5.

######  Remark 3.6.

Along the proof of this result we check some extra properties of the functions
involved. In particular we deduce that,

$\hat{\Gamma}_{1}(L)\in\left[\frac{2}{5}\hat{\Gamma}_{0},\hat{\Gamma}_{0}\right],\qquad\hat{\Gamma}_{2}(L)\in\left[\frac{1}{5}m,\frac{1}{2}m\right],\qquad\hat{\Gamma}_{3}(L)\in\left[\frac{1}{2}m,\frac{3}{4}m\right].$

The constant $\hat{\Gamma}_{0}$ is defined as
$\hat{\Gamma}_{0}=m\hat{y}_{\min}$ with $\hat{y}_{\min}$ defined below in
(5.6).

The exact value of $L_{*}$ is not used along the proof of Theorem 3.5; it has
been computed numerically just for completeness.

###  3.2 The coplanar Hamiltonian

The coplanar Hamiltonian
${\mathcal{H}}_{\mathrm{CP}}=\mathcal{H}_{0}+\alpha^{3}{\mathcal{H}}_{\mathrm{CP},1}$
(see (2.16)) can be seen, when $\alpha^{3}$ is small enough, as a perturbation
of the averaged one. Enlarging the dimension, one can ask if, in some way, the
equilibrium points of the $h$-averaged system, survive as periodic orbit in
the coplanar system and, in this case, if they conserve the character (saddle
or elliptic type) of the equilibrium point.

From expressions of $\mathcal{H}_{0}$ in (2.14) and
${\mathcal{H}}_{\mathrm{CP},1}$ in (A.1) in Appendix A, $\Pi=\\{\eta=\xi=0\\}$
is invariant by the flow associated to the coplanar Hamiltonian
${\mathcal{H}}_{\mathrm{CP}}$. From now on, we focus our analysis on these
types of periodic orbits, which are the ones used in [1] to obtain drift in
eccentricity.

The results presented in this section attempt to give a theoretical framework
of the previous numerical study in [1] by means of perturbative arguments with
respect to the (small) parameter $\delta=\rho\alpha^{3}L^{4}$ defined in
(3.2).

The following result proves the existence of periodic orbits in suitable
energy level of the Hamiltonian $\mathcal{H}_{\mathrm{CP}}$. It is proven in
Section 6.1.

######  Theorem 3.7.

For $L\in(0,L_{\max}]$, the functions $\mathbf{E}_{\min,\max}$ defined by

$\mathbf{E}_{\max}(L)={\mathcal{H}}_{\mathrm{CP}}(0,0,0,0),\qquad\mathbf{E}_{\min}(L)={\mathcal{H}}_{\mathrm{CP}}(0,0.49L,0,\pi)$

can be written as

$\mathbf{E}_{\max}(L)=\frac{\rho_{0}}{16L^{6}}\widehat{\mathbf{E}}_{\max}(\delta),\qquad\mathbf{E}_{\min}(L)=\frac{\rho_{0}}{16L^{6}}\widehat{\mathbf{E}}_{\min}(\delta),$

with $\delta=\rho\alpha^{3}L^{4}$ introduced in (3.2), and
$\widehat{\mathbf{E}}_{\min,\max}:[0,\delta_{\max}]\to\mathbb{R}$.

Then, for any energy level $\mathbf{E}$ such that
$\mathbf{E}_{\min}(L)\leq\mathbf{E}\leq\mathbf{E}_{\max}(L)$, there exists a
periodic orbit $(0,{\Gamma}(t;L),0,h(t;L))$ satisfying that $h(0;L)=0$,
$\Gamma(t;L)\in[0,0.49L]$ for all $t\geq 0$,

$\mathcal{H}_{\mathrm{CP}}(0,{\Gamma}(t;L),0,h(t;L))=\mathbf{E}$

and the differential equations

$\displaystyle\dot{h}=\frac{dh}{dt}=$
$\displaystyle-\frac{3\rho_{0}(L+2\Gamma)}{4L^{8}}-\alpha^{3}\frac{3\rho_{1}}{8L^{4}}\Big{(}U_{2}^{0,0}(L+2\Gamma)+\frac{4}{3}U_{2}^{1,0}\frac{L^{2}-4L\Gamma-4\Gamma^{2}}{\sqrt{(L-2\Gamma)(3L+2\Gamma)}}\cos
h$ $\displaystyle-\frac{1}{3}U_{2}^{2,0}(L+2\Gamma)\cos(2h)\Big{)}$
$\displaystyle\dot{\Gamma}=\frac{d\Gamma}{dt}=$
$\displaystyle-\alpha^{3}\frac{\rho_{1}}{16L^{4}}\Big{(}2U_{2}^{1,0}\sqrt{(L-2\Gamma)(3L+2\Gamma)}(2\Gamma+L)\sin
h$ $\displaystyle+U_{2}^{2,0}(L-2\Gamma)(3L+2\Gamma)\sin(2h)\Big{)}.$

######  Remark 3.8.

The values $\mathbf{E}_{\min}(L),\mathbf{E}_{\max}(L)$ can be easily computed.
For the case $a=29600\,\mathrm{km}$ these are

$\mathbf{E}_{\min}(1)=-2.558100888960067\cdot
10^{-5},\qquad\mathbf{E}_{\max}(1)=2.477266122798186\cdot 10^{-6}.$

Recall that the units are the ones in Remark 2.2, namely the semi-major axis
of Galileo is $a=1$ and the period of Galileo is $2\pi$ (this implies $L=1$).
Compare with the numerical results in [1] for hyperbolic periodic orbits in
Figure 1.

Figure 1: Examples of periodic orbits for $\mathcal{H}_{\mathrm{CP}}$ in the
case $a=29600\,\mathrm{km}$ (Galileo semi-major axis) in non-dimensional
units. The colorbar reports the value of ${\mathcal{H}}_{\mathrm{CP}}$.

We emphasize that, by Theorem 3.7, the initial conditions
$\Gamma(0;L)=\Gamma_{0}$ for a given value $L\in(0,L_{\max}]$ for a periodic
orbit are the ones satisfying

$\mathbf{E}_{\min}(L)\leq\mathcal{H}_{\mathrm{CP}}(0,\Gamma_{0},0,0)\leq\mathbf{E}_{\max}(L),$
(3.9)

with $\mathbf{E}_{\min,\max}(L)$ defined in Theorem 3.7.

For $\sigma~{}\in\left[0,\frac{1}{2}\right)$, we define

$\displaystyle\hat{a}_{0}(\sigma;\delta)$
$\displaystyle=-1-2\sigma-\delta\frac{1}{2}U_{2}^{0,0}(1+2\sigma),$ (3.10)
$\displaystyle\hat{c}_{0}(\sigma;\delta)$
$\displaystyle=-\frac{1}{24\hat{a}_{0}(\sigma;\delta)}\left[4U_{2}^{1,0}\sqrt{(1-2\sigma)(3+2\sigma)}(2\sigma+1)+U_{2}^{2,0}(1-2\sigma)(3+2\sigma)\right].$

######  Remark 3.9.

After some tedious computations (see Section 6.1.1) one can check that

$\displaystyle\hat{a}_{0}\left(\frac{\Gamma}{L};\rho\alpha^{3}L^{4}\right)$
$\displaystyle=\frac{4L^{7}}{3\rho_{0}}\partial_{\Gamma}\mathcal{H}_{\mathrm{AV}}(0,\Gamma,0),$
$\displaystyle\hat{c}_{0}\left(\frac{\Gamma}{L};\rho\alpha^{3}L^{4}\right)$
$\displaystyle=\frac{\rho}{L^{5}}\frac{\mathcal{H}_{\mathrm{CP},1}(0,\Gamma,0,0)-\mathcal{H}_{\mathrm{AV},1}(0,\Gamma,0)}{\partial_{\Gamma}\mathcal{H}_{\mathrm{AV}}(0,\Gamma,0)}.$

The next result is a corollary of Theorem 3.7.

######  Corollary 3.10.

There exists a constant $C_{*}>0$ such that for a given $L\in(0,L_{\max}]$ and
an energy level $\mathbf{E}\in[\mathbf{E}_{\min}(L),\mathbf{E}_{\max}(L)]$, if
$\Gamma_{\mathrm{AV}}^{\mathbf{E}},\Gamma_{\mathrm{CP}}^{\mathbf{E}}$ satisfy

$\mathbf{E}=\mathcal{H}_{\mathrm{AV}}(0,\Gamma_{\mathrm{AV}}^{\mathbf{E}},0)=\mathcal{H}_{\mathrm{CP}}(0,\Gamma_{\mathrm{CP}}^{\mathbf{E}},0,0)$

then there exists a function
$\gamma_{2}:\left(0,\frac{1}{2}\right)\times[0,\delta_{\max}]\to\mathbb{R}$
such that,

$\Gamma_{\mathrm{AV}}^{\mathbf{E}}=\Gamma_{\mathrm{CP}}^{\mathbf{E}}+L\delta\hat{c}_{0}\big{(}\widehat{\Gamma}_{\mathrm{CP}}^{\mathbf{E}};\delta\big{)}+L\delta^{2}\gamma_{2}(\widehat{\Gamma}_{\mathrm{CP}}^{\mathbf{E}};\delta),\quad\widehat{\Gamma}_{\mathrm{CP}}^{\mathbf{E}}=L^{-1}\Gamma_{\mathrm{CP}}^{\mathbf{E}},\quad\delta=\rho\alpha^{3}L^{4},$

with $\hat{c}_{0}$ in (3.10) and
$|\gamma_{2}(\hat{\Gamma}_{\mathrm{CP}}^{\mathbf{E}};\delta)|\leq C_{*}$.

The proof of this result is postponed to Section 6.1.1.

######  Remark 3.11.

Corollary 3.10 gives us a first order in $\delta$ correction of the constant
value of $\Gamma$ to consider in the $h$-averaged Hamiltonian with respect to
the initial condition of the periodic orbit in order to conserve the same
energy level. Such a correction will have a relevant role when the stability
of the periodic orbit is analyzed.

As already mentioned, for Galileo $\delta\sim 0.11$. Then, Corollary 3.10
shows that the two periodic orbits have a significant separation. Moreover,
since at circular motions $G=L$ (see (2.7)), taking (2.10) into account, one
has that

$\Gamma=H-\frac{L}{2}=L\cos i-\frac{L}{2}$

Then, since $L$ is a constant for the secular Hamiltonian, deviations in
$\Gamma$ correspond to deviations of $\cos i$. Namely, Corollary 3.10 implies
that the inclination between the two periodic orbits is significantly
different.

Our next goal is to provide a first order approximation, with respect to the
parameter $\delta$, defined in (3.2), of the periodic orbits
$(0,\Gamma(t;L),0,h(t;L))$ of the coplanar Hamiltonian system
$\mathcal{H}_{\mathrm{CP}}$ emanating from the origin of the $h$-averaged
Hamiltonian (2.18).

To keep our analysis as concrete as possible, before stating the result, we
introduce the functions defined for $\sigma\in\left[0,\frac{1}{2}\right)$

$\displaystyle\hat{c}_{1}(\sigma;\delta)$
$\displaystyle=\frac{1}{6\hat{a}_{0}(\sigma;\delta)}U_{2}^{1,0}\sqrt{(1-2\sigma)(3+2\sigma)}(2\sigma+1),$
(3.11) $\displaystyle\hat{c}_{2}(\sigma;\delta)$
$\displaystyle=\frac{1}{24\hat{a}_{0}(\sigma;\delta)}U_{2}^{2,0}(1-2\sigma)(3+2\sigma),$
$\displaystyle\hat{d}_{1}(\sigma;\delta)$
$\displaystyle=\frac{1}{3\hat{a}_{0}(\sigma;\delta)}U_{2}^{1,0}\frac{5-12\sigma-12\sigma^{2}}{\sqrt{(1-2\sigma)(3+2\sigma)}},$
$\displaystyle\hat{d}_{2}(\sigma;\delta)$
$\displaystyle=\frac{1}{24\hat{a}_{0}(\sigma;\delta)}U_{2}^{2,0}\frac{1-12\sigma-12\sigma^{2}}{1+2\sigma}.$

######  Proposition 3.12.

There exists a constant $C_{*}>0$ such that if $L\in(0,L_{\max}]$ and
$\Gamma_{0}$ satisfy (3.9), then the periodic orbit $(0,\Gamma(t;L),0,h(t;L))$
of $\mathcal{H}_{\mathrm{CP}}$ with initial condition
$\Gamma(0;L)=\Gamma_{0}:=L\hat{\Gamma}_{0}$ is of the form

$\Gamma(t;L)=L\hat{\Gamma}\left({\frac{3\rho_{0}}{4L^{7}}}t;\delta\right),\qquad
h(t;L)=\hat{h}\left({\frac{3\rho_{0}}{4L^{7}}}t;\delta\right)$

where $\delta=\rho\alpha^{3}L^{4}$ as in (3.2) and the functions
$\hat{\Gamma},\hat{h}:[0,\delta_{\max}]\to\mathbb{R}$ are of the form

$\hat{\Gamma}(s;\delta)=\hat{\Gamma}_{0}+\delta\hat{\Gamma}_{1}(s;\delta)+\delta^{2}\hat{\Gamma}_{2}(s;\delta),\qquad\hat{h}(s)=\hat{h}_{0}(s;\delta)+\delta\hat{h}_{1}(s;\delta)+\delta^{2}\hat{h}_{2}(s;\delta)$

with $|\hat{\Gamma}_{2}(s;\delta)|,|\hat{h}_{2}(s;\delta)|\leq C_{*}$ and

$\displaystyle\hat{\Gamma}_{1}(s;\delta)$
$\displaystyle=\hat{c}_{0}(\Gamma_{0};\delta)+\hat{c}_{1}(\hat{\Gamma}_{0};\delta)\cos\big{(}\hat{a}_{0}(\hat{\Gamma}_{0};\delta)s\big{)}+\hat{c}_{2}(\Gamma_{0};\delta)\cos\big{(}2\hat{a}_{0}(\hat{\Gamma}_{0};\delta)s\big{)},$
$\displaystyle\hat{h}_{0}(s;\delta)$
$\displaystyle=\hat{a}_{0}(\hat{\Gamma}_{0};\delta)s,$
$\displaystyle\hat{h}_{1}(s;\delta)$
$\displaystyle=\partial_{\hat{\Gamma}}\hat{a}_{0}(\hat{\Gamma}_{0};\delta)\hat{c}_{0}(\hat{\Gamma}_{0};\delta)s+\hat{d}_{1}(\hat{\Gamma}_{0};\delta)\sin\big{(}\hat{a}_{0}(\hat{\Gamma}_{0};\delta)s\big{)}+\hat{d}_{2}(\hat{\Gamma}_{0};\delta)\sin\big{(}2\hat{a}_{0}(\hat{\Gamma}_{0};\delta)s\big{)}.$

The constants $\hat{a}_{0},\hat{c}_{0,1,2}$ and $\hat{d}_{1,2}$ are defined in
(3.10) and (3.11).

The period $\mathcal{T}(\Gamma_{0};L)$ of the periodic orbit satisfies that

$\mathcal{T}(\Gamma_{0};L)=\frac{4L^{7}}{3\rho_{0}}\widehat{\mathcal{T}}(\hat{\Gamma}_{0};\delta),\qquad\widehat{\mathcal{T}}(\hat{\Gamma}_{0};\delta)=\frac{2\pi+\delta^{2}\widehat{\mathcal{T}_{2}}(\hat{\Gamma}_{0};\delta)}{\big{|}\hat{a}_{0}(\hat{\Gamma}_{0};\delta)+\delta\partial_{\hat{\Gamma}}\hat{a}_{0}(\hat{\Gamma}_{0};\delta)\hat{c}_{0}(\hat{\Gamma}_{0};\delta)\big{|}}$
(3.12)

with
$\widehat{\mathcal{T}}:\left(0,\frac{1}{2}\right)\times[0,\delta_{\max}]\to\mathbb{R}$
and $|\widehat{\mathcal{T}}_{2}(\hat{\Gamma}_{0};\delta)|\leq M$.

This proposition is proven in Section 6.2.1.

######  Remark 3.13.

We notice that the dependence on $L$ of the periodic orbit comes through the
dependence on suitable functions with respect to the “small” parameter
$\delta$. In [1], the period of the hyperbolic periodic orbits are numerically
computed by some values of $\Gamma/L$ for Galileo ($L=1$). Indeed, in Fig. 2
there is the comparison between the numerically computed value of the period
and the approximated value, $\mathcal{T}^{\mathrm{approx}},$ given by formula
(3.12) with $\widehat{\mathcal{T}}_{2}\equiv 0$, namely

$\mathcal{T}^{\mathrm{approx}}(\Gamma_{0};L)=\frac{4L^{7}}{3\rho_{0}}\frac{2\pi}{\big{|}\hat{a}_{0}(\hat{\Gamma}_{0};\delta)+\delta\partial_{\hat{\Gamma}}\hat{a}_{0}(\hat{\Gamma}_{0};\delta)\hat{c}_{0}(\hat{\Gamma}_{0};\delta)\big{|}}$

Figure 2: On the right, the numerically computed period (left) $\mathcal{T}$
of the hyperbolic periodic orbits as a function of
${\mathcal{H}}_{\mathrm{CP}}$. On the left the approximated period,
$\mathcal{T}^{\mathrm{aprox}}$ computed discarding the
$\mathcal{O}(\delta^{2})$ terms in (3.12).

######  Remark 3.14.

From definition (3.10) of $\hat{a}_{0}$, we obtain that, if $\delta$ is small
enough

$\widehat{\mathcal{T}}(\hat{\Gamma}_{0};\delta)=\frac{2\pi}{2\hat{\Gamma}_{0}+1}+\mathcal{O}(\delta)\in\big{(}\pi+\mathcal{O}(\delta),2\pi+\mathcal{O}(\delta)\big{)}.$

Our last result is the comparison between the stability of a periodic orbit of
the coplanar Hamiltonian $\mathcal{H}_{\mathrm{CP}}$ and the stability of the
origin in the $h$-averaged Hamiltonian belonging to the same energy level (see
Corollary 3.10).

In order to analyze the character of the periodic orbits, we denote by
$X_{\mathrm{CP}}$ the vector field of the Hamiltonian
$\mathcal{H}_{\mathrm{CP}}$. For $(0,\Gamma(t;L),0,h(t;L))$, a periodic orbit,
the variational equation around the periodic orbit is given by

$\dot{z}=DX_{\mathrm{CP}}(0,\Gamma(t;L),0,h(t;L))z.$

Let $\Phi(t;\Gamma_{0},L)$ with $\Gamma_{0}=\Gamma(0;L)$, be the fundamental
matrix satisfying the initial condition $\Phi(0;\Gamma_{0},L)=\mathrm{Id}$.
The monodromy matrix is defined just by

$\Phi(\mathcal{T}(\Gamma_{0};L);\Gamma_{0},L)$ (3.13)

with $\Gamma_{0}=\Gamma(0;L)$ and $\mathcal{T}(\Gamma_{0};L)$, the period of
the considered periodic orbit. Since $\mathcal{H}_{\mathrm{CP}}$ is a
2-degrees-of-freedom Hamiltonian, the eigenvalues of the monodromy matrix are
of the form

$\left(\mu_{\mathrm{CP}}(\Gamma_{0};L),\big{(}\mu_{\mathrm{CP}}(\Gamma_{0};L)\big{)}^{-1},1,1\right).$
(3.14)

We are interested now in comparing the eigenvalues of the monodromy matrix
with the ones coming from the constant linear part around the origin of
$\mathcal{H}_{\mathrm{AV}}$, the $h$-averaged Hamiltonian (2.18). That is, we
want to compare $\mu_{\mathrm{CP}}$ with
$e^{\lambda_{\mathrm{AV}}\mathcal{T}}$ where $\lambda_{\mathrm{AV}}$ has been
analyzed in Proposition 3.3.

In the following result, we provide a first order (up to an error of order
$\mathcal{O}(\delta^{3/2})$) approximation of the eigenvalues of the monodromy
matrix.

######  Theorem 3.15.

For $L\in(0,L_{\max}]$ and $\Gamma\in\left(0,\frac{L}{2}\right)$, we define

$\Gamma_{\mathrm{CP}}^{(0)}(\Gamma;L)=L\hat{\Gamma}_{\mathrm{CP}}^{(0)}(\hat{\Gamma};\delta):=L\big{(}\hat{\Gamma}+\delta\hat{c}_{0}(\hat{\Gamma};\delta)),$
(3.15)

with $\hat{c}_{0}$ defined in (3.10) and $\delta=\rho\alpha^{3}L^{4}$. We also
introduce

$\mu_{\mathrm{CP}}^{(0)}(\Gamma;L)=\mathrm{exp}\big{(}\mathcal{T}^{(0)}(\Gamma;L)\lambda_{\mathrm{AV}}(\Gamma;L)\big{)},\qquad\mathcal{T}^{(0)}(\Gamma;L)=\frac{3\rho_{0}}{4L^{7}}\frac{2\pi}{\big{|}\hat{a}_{0}(\hat{\Gamma};\delta)\big{|}},$
(3.16)

where $\hat{a}_{0}$ is defined in (3.10) and $\lambda_{\mathrm{AV}}(\Gamma;L)$
(see Proposition 3.3 ) is an eigenvalue of the linearization around the origin
of the $h$-averaged Hamiltonian $\mathcal{H}_{\mathrm{AV}}$ for constant
values $\Gamma,L$ .

Then, there exist $\delta_{0}>0$ small enough and a constant $C_{*}>0$ such
that, if $L\in(0,L_{\max}]$, $\Gamma_{0}=L\hat{\Gamma}_{0}$ satisfy (3.9) and
$\delta=\rho\alpha^{3}L^{4}\in[0,\delta_{0}]$, we have that

$|\mu_{\mathrm{CP}}(\Gamma_{0};L)-\mu_{\mathrm{CP}}^{(0)}(\Gamma_{\mathrm{CP}}^{(0)}(\Gamma_{0};L);L)|\leq
C_{*}\delta^{3/2},$

where $\mu_{\mathrm{CP}}(\Gamma_{0};L)$ is an eigenvalue of the monodromy
matrix in (3.13) (see also (3.14)) corresponding to the linearized system
around $(0,\Gamma(t;L),0,h(t;L))$, the periodic orbit of
$\mathcal{H}_{\mathrm{CP}}$ with initial condition $\Gamma(0;L)=\Gamma_{0}$,
given in Theorem 3.7.

The proof of this result is postponed to Section 6.2.2 and relies on a first
order perturbation analysis.

######  Corollary 3.16.

For any $\hat{\Gamma}_{0}\in\left(0,\frac{1}{2}\right)$ there exists
$\delta_{0}>0$ such that if $L,\Gamma_{0}=\hat{\Gamma}_{0}L$ satisfy (3.9) and
$\delta=\rho\alpha^{3}L^{4}\in[0,\delta_{0}]$, then the periodic orbit with
initial condition $\Gamma_{0}=\hat{\Gamma}_{0}L$ is of saddle type if
$\hat{\Gamma}_{0}\in(\hat{\Gamma}_{-}(\delta),\hat{\Gamma}_{+}(\delta))$, with
$\hat{\Gamma}_{\pm}(\delta)$ defined in Theorem 3.2, and of elliptic type
otherwise.

In addition, there exists $\delta_{0}>0$, uniform for
$\hat{\Gamma}_{0}\in\left\\{\frac{m}{2}\right\\}\cup\left(0,\frac{m}{4}\right]\cup\left[m,\frac{1}{2}\right)$,
such that the periodic orbit with initial condition
$\Gamma_{0}=L\hat{\Gamma}_{0}$ is of saddle type if
$\hat{\Gamma}_{0}=\frac{m}{2}$ and of elliptic type if
$\hat{\Gamma}_{0}\in\left(0,\frac{m}{4}\right]\cup\left[m,\frac{1}{2}\right)$.

This result is a direct consequence of Theorems 3.15 and 3.2.

Figure 3: The value of $\mu_{\mathrm{CP}}$ (the eigenvalue greater than 1) as
a function of ${\mathcal{H}}_{\mathrm{CP}}$ for the hyperbolic periodic
orbits. On the right, the numerical computations performed in [1], and on the
left, the corresponding approximated value $\mu_{\mathrm{CP}}^{(0)}$ in
Theorem 3.15.

The theoretical approximation result in Theorem 3.15, agrees with the ones
computed numerically in [1]. Indeed, in Figure 3, we present the comparison
between the first order approximated eigenvalue $\mu_{\mathrm{CP}}^{(0)}$ and
the numerical computation of $\mu_{\mathrm{CP}}$ performed in [1]. The
analysis is restricted to values of $\hat{\Gamma}_{0}$ close to
$\frac{m}{2}\sim 0.058257569495584$ so that the periodic orbit could be (and
from the numerical point of view is) hyperbolic.

######  Remark 3.17.

By Corollary 3.10, we obtain that

$\mathcal{H}_{\mathrm{CP}}(0,\Gamma_{0},0,0)=\mathcal{H}_{\mathrm{AV}}(0,\Gamma_{\mathrm{CP}}^{(0)}(\Gamma_{0};L),0)+\mathcal{O}(\delta^{2})$

with $\Gamma_{\mathrm{CP}}^{(0)}$ defined in (3.15). Therefore, as a
straightforward consequence of Theorem 3.15 we conclude that the eigenvalues
of the monodromy matrix of a periodic orbit of period
$\mathcal{T}_{\mathbf{E}}$ lying in an energy level $\mathbf{E}$ are well
approximated by $e^{\mathcal{T}_{\mathbf{E}}\lambda_{\mathbf{E}}}$ with
$\lambda_{\mathbf{E}}$ an eigenvalue associated to the origin as critical
point of the $h$-averaged hamiltonian in the energy level
$\\{\mathcal{H}_{\mathrm{AV}}=\mathbf{E}\\}$ .

Figure 4: Approximated value of the eigenvalue greater than $1$ associated to
the monodromy matrix as a function of $L^{-1}{\Gamma}$. On the left for the
Hamiltonian $\mathcal{H}_{\mathrm{CP}}$ and on the right the corresponding one
for the $h$-averaged Hamiltonian $\mathcal{H}_{\mathrm{AV}}$.

Note, however, that, in terms of the initial condition of the periodic orbit,
namely $\Gamma(0;L)=\Gamma_{0}$, the information about the eigenvalues of the
monodromy matrix associated to a periodic orbit with initial condition
$\Gamma_{0}$ of the coplanar Hamiltonian comes from the $h$-averaged
Hamiltonian corresponding to the corrected value of the parameter $\Gamma$
given by $\Gamma_{\mathrm{CP}}^{(0)}(\Gamma_{0};L)$ in (3.15). To show this
discrepancy, for $L=1$, in Figure 4 they are represented the eigenvalues of
the monodromy matrix associated to the periodic orbit with respect to the
initial condition $\Gamma/L$ as well as
$e^{\mathcal{T}(\Gamma;L)\lambda_{\mathrm{AV}}(\Gamma;L)}$ the eigenvalues of
the linearization of $\mathcal{H}_{\mathrm{AV}}$ corresponding to the same
parameter $\Gamma/L$.

For Galileo, which corresponds to $\mathcal{H}_{\mathrm{CP}}=0$ and semi-major
axis $a=1$ (in the normalized units). The value
$\hat{\Gamma}_{0}=0.058788134221194$ corresponding to the initial condition
$\mathcal{H}_{\mathrm{CP}}(0,L\hat{\Gamma}_{0},0,0)=0$, can be numerically
computed. Then the corrected value for $\Gamma$ is
$\Gamma_{\mathrm{CP}}^{(0)}(\Gamma_{0};L)=0.079342370619096.$ Therefore, we
have that the approximated eigenvalue is

$\mu_{\mathrm{CP}}^{(0)}(\Gamma_{\mathrm{CP}}^{(0)}(\Gamma_{0};L);L)=1.610955038638576.$

However, if the correction in $\Gamma$ is not considered in the $h$-averaged
Hamiltonian and the period is only approximated by the averaged Hamiltonian,
we obtain

$e^{\lambda_{\mathrm{AV}}(\Gamma_{0};L)\mathcal{T}(\Gamma_{0};L)}=2.038395173606618.$

Therefore, the difference between the eigenvalues of the monodromy and the
corresponding orbit in the $h$-averaged system, namely with the same initial
condition $\Gamma_{0}$, can not be _a priori_ neglected.

The last result of this work can be seen as a refinement of Theorem 3.15. It
gives a more accurate bound for the eigenvalues of the monodromy matrix for
the most part of the values of $L,\Gamma_{0}$. In order to state such a
result, we introduce the values
$\widehat{\Gamma}^{(k)}_{0}\in\left(0,\frac{1}{2}\right)$ satisfying that

$\frac{2\pi}{\hat{a}_{0}\big{(}\widehat{\Gamma}^{(k)}_{0};0\big{)}}|1-16\widehat{\Gamma}^{(k)}_{0}-20(\widehat{\Gamma}^{(k)}_{0})^{2}|=k\pi,\qquad
k=1,2.$ (3.17)

That is
$|1-16\widehat{\Gamma}^{(k)}_{0}-20(\widehat{\Gamma}^{(k)}_{0})^{2}|=2(1+2\widehat{\Gamma}^{(k)}_{0})k$.
These values, $\widehat{\Gamma}^{(k)}_{0}$, can be analytically computed (see
Table 1 for the numerical value) because they are the zeros of suitable
polynomial of degree $2$.

$k$ | $\hat{\Gamma}^{(k)}_{0}$
---|---
$1$ | $0.189897948556636$
$2$ | $0.338516480713450$
Table 1:  The values of $\widehat{\Gamma}_{0}^{(k)}$. For $k=3$,
$\widehat{\Gamma}_{0}^{(3)}=1/2$ and for $k\geq 4$, the values of
$\widehat{\Gamma}_{0}^{(k)}$ are out of the interval
$\left(0,\frac{1}{2}\right)$.

We also introduce the constant $C_{0}$ defined through (3.5) in Theorem 3.2,
as

$C_{0}=\frac{5U_{2}^{1,0}(m+3)\sqrt{3-2m-4m^{2}}}{3(16+20m)}$ (3.18)

in such a way that $\Gamma_{\pm}(L)L^{-1}=\frac{m}{2}\mp
C_{0}\delta+\mathcal{O}(\delta^{2})$. We recall that
$\lambda_{\mathrm{AV}}(\Gamma_{\pm}(L);L)=0$.

######  Lemma 3.18.

There exists $\delta_{0}>0$ such that if $L$ satisfies that
$\delta=\rho\alpha^{3}L^{4}\in[0,\delta_{0}]$ then, there exist
$\Gamma^{(k)}(L)$, $k=1,2$ such that

$\mu_{\mathrm{CP}}^{(0)}(\Gamma^{(k)}(L);L)=(-1)^{k},$

where $\Gamma^{(k)}(L)=L\widehat{\Gamma}^{(k)}(\delta)$ (with
$\delta=\rho\alpha^{3}L^{4}$) satisfy

$\big{|}\widehat{\Gamma}^{(k)}(\delta)-\widehat{\Gamma}^{(k)}_{0}\big{|}\leq
M\delta^{2},$

with the value $\hat{\Gamma}^{(k)}_{0}$ defined by (3.17) (see also Table 1)
and the constant $M$ is independent of $L$.

In addition, if either $\Gamma\neq\Gamma^{(k)}(L)$ or
$\Gamma\neq\Gamma_{\pm}(L)$, then $\mu_{\mathrm{CP}}^{(0)}(\Gamma;L)\neq\pm
1$.

This lemma is proven in Section 6.2.3.

######  Theorem 3.19.

Assume that we are under the conditions of Theorem 3.15 and denote
$\Gamma_{\mathrm{CP}}^{(0)}(\Gamma_{0};L)=L\hat{\Gamma}_{\mathrm{CP}}^{(0)}(\hat{\Gamma}_{0};\delta)$,
as in (3.15).

Fix $C\in(0,C_{0})$ (see (3.18)), $\nu\in[1,2]$ and consider the sets, for
$k=1,2$, defined for $\delta>0$ as

$I_{\pm,\nu}^{(0)}(\delta)=\left\\{\sigma\in\mathbb{R},\,\left|\sigma-\frac{m}{2}\mp
C_{0}\delta\right|\leq
C\delta^{\nu}\right\\},\;I^{(k)}_{\nu}(\delta)=\left\\{\sigma\in\mathbb{R},\,\left|\sigma-\widehat{\Gamma}_{0}^{(k)}\right|\leq
C\delta^{\frac{1+\nu}{2}}\right\\}.$

There exists a constant $C_{*}>0$ and $\delta_{0}>0$ small enough, such that,
for $\delta\in[0,\delta_{0}]$, if
$\hat{\Gamma}_{\mathrm{CP}}^{(0)}(\hat{\Gamma}_{0};\delta)\notin
I_{\pm,\nu}^{(0)}(\delta)\cup I_{\nu}^{(1)}(\delta)\cup I_{\nu}^{(2)}(\delta)$
then

$\big{|}\mu_{\mathrm{CP}}(\Gamma_{0};L)-\mu_{\mathrm{CP}}^{(0)}(\Gamma^{(0)}_{\mathrm{CP}}(\Gamma_{0};L);L)\big{|}\leq
C_{*}\delta^{\frac{5-\nu}{2}}.$

As a consequence when
$\hat{\Gamma}_{\mathrm{CP}}^{(0)}(\hat{\Gamma}_{0};\delta)\notin
I_{\pm,1}^{(0)}(\delta)\cup I_{1}^{(1)}(\delta)\cup I_{1}^{(2)}(\delta)$, then

$\big{|}\mu_{\mathrm{CP}}(\Gamma_{0};L)-\mu_{\mathrm{CP}}^{(0)}(\Gamma^{(0)}_{\mathrm{CP}}(\Gamma_{0};L);L)\big{|}\leq
C_{*}\delta^{2}.$

This result is proven in Section 6.2.4.

The remaining part of this work is devoted to prove all the results stated in
this section.

##  4 The circular critical points of the $h$-averaged Hamiltonian

In Section 4.1, we rewrite the averaged system in a more suitable way to
analyze the critical points. This new expression will be used both to study
the circular critical points and the eccentric ones. After that, in Sections
4.2 and 4.3, we prove Theorem 3.2 and Proposition 3.3 respectively.

###  4.1 Rewriting the $h$-averaged system

We start with the $h$-averaged Hamiltonian in slow-fast Delaunay variables
$(y,x)$. That is, we consider (see (2.20))

$\displaystyle{\mathrm{H}}_{\mathrm{AV}}(y,\Gamma,x;L)$
$\displaystyle={\mathrm{H}}_{0}(y,\Gamma;L)+\alpha^{3}{\mathrm{H}}_{\mathrm{AV},1}(y,\Gamma,x;L)$
$\displaystyle=\frac{\rho_{0}}{128}\frac{y^{2}-6y\Gamma-3\Gamma^{2}}{L^{3}y^{5}}+\alpha^{3}\frac{\rho_{1}}{L^{2}}\left[\frac{1}{2}U_{2}^{0,0}D_{0,1}(y,\Gamma)+\frac{1}{3}U_{2}^{1,0}D_{1,0}(y,\Gamma)\cos
x\right].$

with $D_{0,1}$ and $D_{1,0}$ as given in Table 5. In order to rewrite it in a
more concise form, we introduce the following polynomials in $(y,\Gamma)$,
which depend implicitly on the parameter $L$,

$\begin{array}[]{lllll}P_{0}(y,\Gamma)&=y^{2}-6y\Gamma-3\Gamma^{2},&&Q_{0}(y;L)&=(5L^{2}-12y^{2})y^{3},\\\
P_{1}(y,\Gamma)&=(y-\Gamma)(3y+\Gamma),&&Q_{1}(y,\Gamma;L)&=(3y+\Gamma)(L-2y)(L+2y)y^{3},\end{array}$

and their derivatives with respect to $y$,

$\begin{array}[]{lllll}P_{0}^{\prime}(y,\Gamma)&=2(y-3\Gamma),&&Q_{0}^{\prime}(y;L)&=15(L-2y)(L+2y)y^{2},\\\
P_{1}^{\prime}(y,\Gamma)&=2(3y-\Gamma),&&Q_{1}^{\prime}(y,\Gamma;L)&=y^{2}\bigg{[}3(L-2y)(L+2y)(4y+\Gamma)-8y^{2}(3y+\Gamma)\bigg{]}.\end{array}$

Hence,the $h$-averaged Hamiltonian as well as its derived equations of motion
can be written as follows

$\begin{split}{\mathrm{H}}_{\mathrm{AV}}(y,\Gamma,x;L)=\frac{1}{128}\frac{\rho_{0}}{L^{4}y^{5}}\left({\mathcal{A}}_{0,\alpha}(y,\Gamma)+\alpha^{3}\sqrt{P_{1}(y,\Gamma)}{\mathcal{A}}_{1}(y,\Gamma;L)\cos
x\right)\end{split}$

and

$\begin{split}\dot{x}&=-\frac{1}{256}\frac{\rho_{0}}{L^{4}y^{6}}\frac{1}{\sqrt{P_{1}(y,\Gamma)}}\left[\sqrt{P_{1}(y,\Gamma)}{\mathcal{B}}_{0,\alpha}(y,\Gamma;L)+\alpha^{3}{\mathcal{B}}_{1}(y,\Gamma;L)\cos
x\right]\\\
\dot{y}&=\phantom{-}\frac{1}{128}\frac{\rho_{0}\alpha^{3}}{L^{4}y^{5}}\sqrt{P_{1}(y,\Gamma)}{\mathcal{A}}_{1}(y,\Gamma;L)\sin
x,\end{split}$ (4.1)

where the functions

$\begin{split}{\mathcal{A}}_{0,\alpha}&=P_{0}\left(L+\alpha^{3}d_{0}Q_{0}\right),\\\
{\mathcal{A}}_{1}&=d_{1}Q_{1},\\\
{\mathcal{B}}_{0,\alpha}&=10P_{0}\left(L+\alpha^{3}d_{0}Q_{0}\right)-2y(L+\alpha^{3}d_{0}Q_{0})P_{0}^{\prime}-2\alpha^{3}yd_{0}P_{0}Q_{0}^{\prime},\\\
{\mathcal{B}}_{1}&=d_{1}(10P_{1}Q_{1}-yQ_{1}P_{1}^{\prime}-2yP_{1}Q_{1}^{\prime})\end{split}$

are polynomial in $(y,\Gamma)$, with

$d_{0}=2U_{2}^{0,0}\frac{\rho_{1}}{\rho_{0}},\quad
d_{1}=20U_{2}^{1,0}\frac{\rho_{1}}{\rho_{0}}.$ (4.2)

To analyze the circular critical points (i.e., at $e=0$), we work with
Poincaré variables. In this case, the $h$-averaged Hamiltonian (2.18) as well
as its derived equations of motion can be written as follows

$\begin{split}\mathcal{H}_{\mathrm{AV}}(\xi,\eta)&=\frac{1}{128}\frac{\rho_{0}}{L^{4}y^{5}}\left({\mathcal{A}}_{0,\alpha}+\alpha^{3}d_{1}\sqrt{P_{1}}\tilde{Q}_{1}(\xi^{2}-\eta^{2})\right)\end{split}$

and

$\begin{split}\dot{\xi}&=-\frac{1}{512}\frac{\rho_{0}}{L^{4}y^{6}}\frac{\eta}{\sqrt{P_{1}}}\left[\sqrt{P_{1}}{\mathcal{B}}_{0,\alpha}-8y\alpha^{3}d_{1}P_{1}\tilde{Q}_{1}+\alpha^{3}{\tilde{{\mathcal{B}}}}_{1}\times(\xi^{2}-\eta^{2})\right]\\\
\dot{\eta}&=\frac{1}{512}\frac{\rho_{0}}{L^{4}y^{6}}\frac{\xi}{\sqrt{P_{1}}}\left[\sqrt{P_{1}}{\mathcal{B}}_{0,\alpha}+8y\alpha^{3}d_{1}P_{1}\tilde{Q}_{1}+\alpha^{3}{\tilde{{\mathcal{B}}}}_{1}\times(\xi^{2}-\eta^{2})\right]\\\
\end{split}$ (4.3)

with

$y=\frac{L}{2}-\frac{\xi^{2}+\eta^{2}}{4},$

and

$\begin{split}\tilde{Q}_{1}&=\frac{Q_{1}}{2L-4y}=\frac{1}{2}(3y+\Gamma)(L+2y)y^{3},\\\
\tilde{Q}_{1}^{\prime}&=\frac{1}{2}y^{2}\bigg{[}3L(\Gamma+4y)+2y(15y+4\Gamma)\bigg{]},\\\
{\tilde{{\mathcal{B}}}}_{1}&=d_{1}(10P_{1}\tilde{Q}_{1}-y\tilde{Q}_{1}P_{1}^{\prime}-2yP_{1}\tilde{Q}_{1}^{\prime}).\end{split}$

For future purposes, we also decompose

${\mathcal{B}}_{0,\alpha}=B_{0}+\alpha^{3}B_{1}$ (4.4)

with

$B_{0}=L(10P_{0}-2yP_{0}^{\prime}),\qquad
B_{1}=d_{0}(10P_{0}Q_{0}-2yQ_{0}P_{0}^{\prime}-2yP_{0}Q_{0}^{\prime})$

and we recall that

$U_{2}^{0,0}=0.762646,\qquad
U_{2}^{1,0}=0.547442,\qquad\rho:=\frac{\rho_{1}}{\rho_{0}}=1.762157978551987\cdot
10^{-18}.$ (4.5)

We also recall that $a_{\max}=30000$ km is the maximum semi-major axis we are
going to consider and that $a_{\min}=6378.14$ km (the radius of the Earth) is
the minimum one.

Along the proof of the results this notation will be used extensively without
an explicit mention.

###  4.2 Circular critical points: Proof of Theorem 3.2

The origin, $(\xi,\eta)=(0,0)$, is clearly an equilibrium point of system
(4.3). To prove the statements in Theorem 3.2, we only need to study the
linear part of the averaged system at the origin and elucidate the values of
the parameters $L,\Gamma$ for which the origin is either a saddle, an elliptic
point or a degenerated (parabolic) fixed point.

We observe that, since

$y=y(\xi,\eta)=\frac{L}{2}-\frac{1}{4}(\xi^{2}+\eta^{2})$

satisfies $y(0,0)=L/2$ and $\partial_{\xi}y(0,0)=\partial_{\eta}y(0,0)=0$, the
variational equation of the averaged system around $(\xi,\eta)=(0,0)$ is
$\dot{z}=M(\Gamma;L)z$ where

$M(\Gamma;L)=g(L)\left(\begin{array}[]{cc}0&-X^{\eta}(\Gamma;L)\\\
X^{\xi}(\Gamma;L)&0\end{array}\right)$ (4.6)

and

$\displaystyle g(L)$ $\displaystyle=\frac{\rho_{0}}{8L^{10}}$ $\displaystyle
X^{\xi}(\Gamma;L)$
$\displaystyle=B_{0}\left(\frac{L}{2},\Gamma;L\right)+\alpha^{3}B_{1}\left(\frac{L}{2},\Gamma;L\right)+4\alpha^{3}d_{1}L\sqrt{P_{1}\left(\frac{L}{2},\Gamma\right)}\tilde{Q}_{1}\left(\frac{L}{2},\Gamma;L\right)$
$\displaystyle X^{\eta}(\Gamma;L)$
$\displaystyle=B_{0}\left(\frac{L}{2},\Gamma;L\right)+\alpha^{3}B_{1}\left(\frac{L}{2},\Gamma;L\right)-4\alpha^{3}d_{1}L\sqrt{P_{1}\left(\frac{L}{2},\Gamma\right)}\tilde{Q}_{1}\left(\frac{L}{2},\Gamma;L\right).$

Therefore the eigenvalues $\lambda$ of $M(\Gamma;L)$ satisfy

$\lambda^{2}=-g^{2}(L)X^{\xi}(\Gamma;L)\cdot X^{\eta}(\Gamma;L)$ (4.7)

so that, to determine the character of $(\xi,\eta)=(0,0)$, we need to study
the sign of the product $X^{\xi}(\Gamma;L)\cdot X^{\eta}(\Gamma;L)$ with
respect to the parameters $\Gamma,L$. To this end, we perform the scaling

$\Gamma=L\sigma,\qquad\text{with}\qquad 0<\sigma<\frac{1}{2}$

and compute $X^{\xi}(L\sigma;L)\cdot X^{\eta}(L\sigma;L)$. One can check
(after some tedious computations) that

$\displaystyle
X^{\xi}(L\sigma;L)=\frac{3}{2}L^{3}\left[\beta(L)\mathbf{a}(\sigma)+\frac{1}{12}d_{1}\alpha^{3}L^{4}\mathbf{b}(\sigma)\right]$
(4.8) $\displaystyle
X^{\eta}(L\sigma;L)=\frac{3}{2}L^{3}\left[\beta(L)\mathbf{a}(\sigma)-\frac{1}{12}d_{1}\alpha^{3}L^{4}\mathbf{b}(\sigma)\right]$

with

$\mathbf{a}(\sigma)=1-16\sigma-20\sigma^{2},\qquad\mathbf{b}(\sigma)=(2\sigma+3)\sqrt{3-4\sigma-4\sigma^{2}}$
(4.9)

(note that $p(\sigma):=3-4\sigma-4\sigma^{2}>0$ if $0<\sigma<\frac{1}{2}$ and
$p\left(\frac{1}{2}\right)=0$) and

$\beta(L)=1+\frac{d_{0}\alpha^{3}L^{4}}{4}.$

######  Remark 4.1.

Notice that the eigenvalues $\lambda=\lambda(\sigma;L)$ can be explicitly
computed. Indeed, using (4.7) and definitions (4.2) and (4.5) for
$d_{0},d_{1}$ and $\rho$ respectively

$\lambda^{2}=-\frac{9\rho_{0}^{2}}{256L^{14}}\left[\mathbf{a}(\sigma)+\rho\alpha^{3}L^{4}\mathbf{c}_{-}(\sigma)\right]\left[\mathbf{a}(\sigma)+\rho\alpha^{3}L^{4}\mathbf{c}_{+}(\sigma)\right]$

with

$\mathbf{c}_{\pm}(\sigma)=\frac{U_{2}^{0,0}}{2}\mathbf{a}(\sigma)\pm\frac{5U_{2}^{1,0}}{3}\mathbf{b}(\sigma).$

######  Lemma 4.2.

Consider

$\mathcal{X}_{\pm}(\sigma;L):=\beta(L)\mathbf{a}(\sigma)\pm\frac{1}{12}d_{1}\alpha^{3}L^{4}\mathbf{b}\left(\sigma\right).$

For $L\in[0,L_{\max}]$, the function $\mathcal{X}_{+}(\cdot;L)$ is strictly
decreasing for $\sigma\in\left(0,\frac{1}{2}\right)$ and
$\mathcal{X}_{-}(\cdot;L)$ is strictly decreasing for
$\sigma\in\left(0,\frac{m}{2}\right)$.

Moreover, there exist two functions
$\sigma_{\pm}:[0,L_{\max}]\to(0,\frac{1}{2})$ such that for any
$L\in[0,L_{\max}]$ $\mathcal{X}_{\pm}(\sigma;L)=0$ if and only if
$\sigma=\sigma_{\pm}(L)$. In addition, $\sigma_{\pm}(0)=\frac{m}{2}$ and for
$L>0$

$\sigma_{-}(L)\in\left(0,\frac{m}{2}\right),\qquad\sigma_{+}(L)\in\left(\frac{m}{2},\frac{1}{4}\right),$

where $m$, defined in (3.4), is the slope of the prograde resonance (see also
(2.12)).

###### Proof.

First we notice that $\mathbf{a},\mathbf{b}$ are decreasing functions for
$\sigma\in(0,1/2)$. Indeed, we only need to compute

$\partial_{\sigma}\mathbf{b}(\sigma)=-8\frac{\sigma(2\sigma+3)}{\sqrt{3-4\sigma-4\sigma^{2}}}<0$

(for $\mathbf{a}$ is obvious). Therefore, $\mathcal{X}_{+}(\sigma;L)$ is an
strictly decreasing function (with respect to $\sigma$) and has at most a
unique zero. Using that $\mathbf{a}\left(\frac{m}{2}\right)=0$ we obtain that

$\mathcal{X}_{+}\left(\frac{m}{2};L\right)=\frac{1}{12}d_{1}\alpha^{3}L^{4}\mathbf{b}\left(\frac{m}{2}\right)>0$

and

$\mathcal{X}_{+}\left(\frac{1}{4};L\right)\leq\mathbf{a}\left(\frac{1}{4}\right)+\frac{1}{12}d_{1}\alpha_{\max}^{3}L_{\max}^{4}\mathbf{b}\left(\frac{1}{4}\right)=-3.744000070111437<0,$

with $\alpha_{\max}=a_{\max}/a_{{\mathrm{M}}}$, and by Bolzano theorem we
obtain the conclusion for $\mathcal{X}_{+}$.

Concerning $\mathcal{X}_{-}(\sigma;L)$, we have that

$\mathcal{X}_{-}\left(\frac{m}{2};L\right)=-\frac{1}{12}d_{1}\alpha^{3}L^{4}\mathbf{b}\left(\frac{m}{2}\right)<0$

and

$\mathcal{X}_{-}(0;L)=\beta(L)-\frac{1}{12}d_{1}\alpha^{3}L^{4}3\sqrt{3}=1+\frac{\alpha^{3}L^{4}}{4}\left(d_{0}-d_{1}\sqrt{3}\right).$

Therefore, since $d_{0}-d_{1}\sqrt{3}\sim-17.4387\rho<0$,

$\mathcal{X}_{-}(0;L)\geq
1+\frac{\alpha_{\max}^{3}L_{\max}^{4}}{4}\left(d_{0}-d_{1}\sqrt{3}\right)\geq
0.477808830620940>0.$

Again there exists $\sigma_{-}(L)\in\left(0,\frac{m}{2}\right)$ such that
$\mathcal{X}_{-}(\sigma_{-}(L);L)=0$. In addition, one can easily check that
$\partial_{\sigma}\mathcal{X}_{-}(\sigma;L)<0$ if
$\sigma\in\left(0,\frac{m}{2}\right)$ so that $\sigma=\sigma_{-}(L)$ is the
only solution of $\mathcal{X}_{-}(\sigma;L)=0$ belonging to
$\left(0,\frac{m}{2}\right)$.

To finish with this analysis, we point out that, for
$\sigma\in\left(\frac{m}{2},\frac{1}{2}\right)$,

$\mathcal{X}_{-}(\sigma;L)=\beta(L)\mathbf{a}(\sigma)-\frac{1}{12}d_{1}\alpha^{3}L^{4}\mathbf{b}(\sigma)<0$

provided $\beta(L)>0$, $\mathbf{a}(\sigma)<0$ and $\mathbf{b}(\sigma)>0$ if
$\sigma\in\left(\frac{m}{2},\frac{1}{2}\right)$. ∎

From formula (4.7) of the eigenvalues and the previous lemma, recalling that
$\Gamma=L\sigma$, we have that

  * •

If either $\Gamma\in(0,L\sigma_{-}(L))$ or
$\Gamma\in\left(L\sigma_{+}(L),\frac{L}{2}\right)$, then $(\xi,\eta)=(0,0)$ is
a center equilibrium point.

  * •

If $\Gamma\in(L\sigma_{-}(L),L\sigma_{+}(L))$, then $(\xi,\eta)=(0,0)$ is a
saddle equilibrium point.

  * •

Using that $\mathcal{X}_{+}(\cdot;L)$ and $\mathcal{X}_{-}(\cdot;L)$ are not
zero simultaneously, we deduce that, when $\Gamma=L\sigma_{-}(L)$ or
$\Gamma=L\sigma_{+}(L)$, then $(\xi,\eta)=(0,0)$ is a degenerated equilibrium
point with nilpotent linear part.

Therefore, in order to prove Theorem 3.2, we only need to prove the expansions
of $\sigma_{\pm}(L)$ in (3.5). Indeed, $\sigma_{\pm}(L)$ satisfy

$\mathcal{X}_{\pm}(\sigma_{\pm}(L);L)=\beta(L)\mathbf{a}(\sigma_{\pm}(L))\pm\frac{1}{12}d_{1}\alpha^{3}L^{4}\mathbf{b}(\sigma_{\pm}(L))=0.$
(4.10)

We rewrite condition (4.10) in a more suitable way. To do so, denoting
$\delta=\rho\alpha^{3}L^{4}$, we introduce the functions

$\mathbf{A}_{\pm}(\sigma;\delta):=\mathbf{a}(\sigma)\left(1+\frac{U_{2}^{0,0}}{2}\delta\right)\pm\frac{5U_{2}^{1,0}}{3}\delta\mathbf{b}(\sigma),$

which are smooth for
$(\sigma,\delta)\in\left(0,\frac{1}{2}\right)\times(0,\infty)$. Then, by the
definition of $d_{0},d_{1}$ in (4.2) and also (4.5), we have that,
$\mathcal{X}_{\pm}(\sigma;L)=\mathbf{A}_{\pm}(\sigma;\delta)$ so that
expression (4.10) is equivalent to

$\mathbf{A}_{\pm}(\sigma_{\pm}(L);\delta)=0.$

Since $\mathbf{A}_{\pm}\left(\frac{m}{2};0\right)=0$ and

$\partial_{\sigma}\mathbf{A}_{\pm}\left(\frac{m}{2};0\right)=\partial_{\sigma}\mathbf{a}\left(\frac{m}{2}\right)=-16-20m\neq
0,$

by the Implicit Function Theorem, $\sigma_{\pm}(L)$ is, in fact, a smooth
function of $\delta=\rho\alpha^{3}L^{4}$, namely
$\sigma_{\pm}(L)=\tilde{\sigma}_{\pm}(\delta)$ with
$\tilde{\sigma}(0)=\frac{m}{2}$. Therefore,

$\sigma_{\pm}(L)=\tilde{\sigma}_{\pm}(\delta)=\frac{m}{2}+c_{\pm}\delta+\mathcal{O}(\delta^{2})$

with

$\displaystyle c_{\pm}$
$\displaystyle=\partial_{\delta}\tilde{\sigma}_{\pm}(0)=-\frac{\partial_{\delta}\mathbf{A}_{\pm}\left(\frac{m}{2};0\right)}{\partial_{\sigma}\mathbf{A}_{\pm}\left(\frac{m}{2};0\right)}=\frac{1}{16+20m}\left(\mathbf{a}\left(\frac{m}{2}\right)\frac{U_{2}^{0,0}}{2}\pm\frac{5U_{2}^{1,0}}{3}\mathbf{b}\left(\frac{m}{2}\right)\right)$
$\displaystyle=\pm\frac{1}{16+20m}\frac{5U_{2}^{1,0}}{3}\mathbf{b}\left(\frac{m}{2}\right).$

Using that, by Lemma 4.2, $\mathcal{X}_{+}$ is a strictly decreasing function
(with respect to $\sigma$) and that $\mathcal{X}_{-}$ is also a strictly
decreasing function if $\sigma\in\left(0,\frac{m}{2}\right)$, we conclude that
$\tilde{\sigma}_{\pm}$ are well defined for
$\delta\in\big{[}0,\rho\alpha_{\max}L_{\max}^{4}\big{]}$ and the proof of
Theorem 3.2 is complete.

###  4.3 Modulus of the eigenvalues: Proof of Proposition 3.3

We use, along this section, the notation introduced in Section 4.2. We first
define, from (4.7) and (4.8),

$\mathcal{M}(\sigma;L):=-g^{2}(L)\frac{9L^{6}}{4}\left[\beta^{2}(L)\mathbf{a}^{2}(\sigma)-\widetilde{\beta}^{2}(L)\mathbf{b}^{2}(\sigma)\right],$
(4.11)

with $g(L)=\frac{\rho_{0}}{8L^{10}}$,
$\widetilde{\beta}(L)=\frac{d_{1}\alpha^{3}L^{4}}{12}$ and $\mathcal{X}_{\pm}$
were introduced in Lemma 4.2. It is not difficult to check that

$\displaystyle\partial_{\sigma}\mathcal{M}(\sigma;L)$
$\displaystyle=\frac{-9\rho_{0}^{2}}{256L^{14}}\left[2\beta^{2}(L)\mathbf{a}(\sigma)\mathbf{a}^{\prime}(\sigma)-2\widetilde{\beta}^{2}(L)\mathbf{b}(\sigma)\mathbf{b}^{\prime}(\sigma)\right]$
$\displaystyle=\frac{-9\rho_{0}^{2}}{256L^{14}}\left[\beta^{2}(L)(1600\sigma^{3}+1920\sigma^{2}+432\sigma-32)+\widetilde{\beta}^{2}(L)(16\sigma(2\sigma+3)^{2})\right].$

Therefore $\partial_{\sigma}\mathcal{M}$, for any fixed $L$, is a degree three
polynomial such that its derivative has no positive zero (all its coefficients
are positive). That implies that $\partial_{\sigma}\mathcal{M}(\cdot;L)$ can
at most have one zero for any value of $L$. Since, by definition,
$\mathcal{M}(\sigma_{+}(L);L)=\mathcal{M}(\sigma_{-}(L);L)=0$, by Rolle’s
theorem, there exists
$\sigma_{*}(L)\in\big{(}\sigma_{-}(L),\sigma_{+}(L)\big{)}$ such that
$\partial_{\sigma}\mathcal{M}(\sigma_{*}(L);L)=0$. Moreover, since
$\mathbf{a}\left(\frac{m}{2}\right)=0$,

$\partial_{\sigma}\mathcal{M}(0;L)>0,\qquad\partial_{\sigma}\mathcal{M}\left(\frac{m}{2};L\right)<0,$

we conclude that, for any $L_{\min}\leq L\leq L_{\max}$, the function
$\mathcal{M}(\sigma;L)$ has a maximum at $\sigma_{*}(L)$. Therefore, for a
fixed value of $L$,

  * •

on the interval $[\sigma_{-}(L),\sigma_{+}(L)]$, the function
$\lambda_{+}(\sigma;L)=\sqrt{\mathcal{M}(\sigma;L)}$ has a maximum at
$\sigma=\sigma_{*}(L)$.

  * •

When $\sigma\in(0,\sigma_{-}(L))$ the function
$\sqrt{\big{|}\mathcal{M}(\sigma;L)\big{|}}$ is decreasing with respect to
$\sigma$

  * •

and when $\sigma\in\left(\sigma_{+}(L),\frac{1}{2}\right)$,
$\sqrt{\big{|}\mathcal{M}(\sigma;L)\big{|}}$ is an increasing function.

To finish the proof of Proposition 3.3, we notice that, from (4.11) and using
that $\mathbf{b}$ is a decreasing positive function, we deduce that, if
$\sigma\in(\sigma_{-}(L),\sigma_{+}(L))$, then the corresponding eigenvalues
$\pm\lambda(\sigma;L)$ satisfy

$0<\lambda(\sigma;L)\leq
g(L)\frac{3L^{3}}{24}d_{1}\alpha^{3}L^{4}\mathbf{b}(0)\leq\frac{\rho_{0}}{8L^{10}}\frac{3L^{3}}{24}d_{1}\alpha^{3}L^{4}3\sqrt{3}=\frac{15\sqrt{3}\rho_{0}}{16L^{7}}U_{2}^{1,0}\rho\alpha^{3}L^{4},$

with $\rho=\rho_{1}/\rho_{0}$. In addition, since
$\lambda(\sigma_{*}(L);L)>\lambda\left(\sigma\left(\frac{m}{2}\right);L\right)$,

$\lambda(\sigma_{*}(L);L)\geq
g(L)\frac{3L^{3}}{24}d_{1}\alpha^{3}L^{4}\mathbf{b}\left(\frac{m}{2}\right)=\frac{5\rho_{0}}{16L^{7}}U_{2}^{1,0}\rho\alpha^{3}L^{4}(m+3)\sqrt{3-2m-m^{2}}.$

######  Remark 4.3.

As a consequence of our study, if we want to analyze the value of
$\mathcal{M}(\sigma;L)$ at some interval
$\sigma\in[\sigma_{\min},\sigma_{\max}]\subset[\sigma_{-}(L),\sigma_{+}(L)]$,
its minimum value is located either at $\sigma=\sigma_{\min}$ or
$\sigma=\sigma_{\max}$.

######  Remark 4.4.

Using $L=\sqrt{\mu a}$ and $\alpha=\frac{a}{a_{M}}$, we write from (4.7),
$\widetilde{\mathcal{M}}(\sigma;a):=\mathcal{M}(\sigma;\sqrt{\mu a})$ as

$\displaystyle\widetilde{\mathcal{M}}(\sigma;a)$
$\displaystyle:=-\frac{9\rho_{0}^{2}}{256L^{14}}\left[\left(1+\frac{d_{0}\alpha^{3}L^{4}}{4}\right)^{2}\mathbf{a}^{2}(\sigma)-\frac{1}{144}d_{1}^{2}\alpha^{6}L^{8}\mathbf{b}^{2}(\sigma)\right]$
$\displaystyle=-\frac{9\rho_{0}^{2}}{256\mu^{7}a^{7}}\left[\mathbf{a}^{2}(\sigma)\left(1+\frac{d_{0}a^{5}\mu^{2}}{2a_{M}^{3}}+\frac{d_{0}^{2}a^{10}\mu^{4}}{16a_{M}^{6}}\right)-\frac{1}{144}d_{1}^{2}\mathbf{b}^{2}(\sigma)\frac{a^{10}\mu^{4}}{a_{M}^{6}}\right]$
$\displaystyle=-\frac{9\rho_{0}^{2}}{256\mu^{7}}\left[\mathbf{a}^{2}(\sigma)\left(\frac{1}{a^{7}}+\frac{d_{0}\mu^{2}}{2a^{2}a_{M}^{3}}+\frac{d_{0}^{2}a^{3}\mu^{4}}{16a_{M}^{6}}\right)-\frac{1}{144}d_{1}^{2}\mathbf{b}^{2}(\sigma)\frac{a^{3}\mu^{4}}{a_{M}^{6}}\right].$

We have that $\partial_{a}\widetilde{\mathcal{M}}(\sigma;a)>0$ provided that,
for $a\leq a_{\max}$,

$-\frac{7}{a^{8}}-\frac{d_{0}\mu^{2}}{a^{3}a_{M}^{3}}+\frac{3d_{0}^{2}a^{2}\mu^{4}}{16a_{M}^{6}}=-\frac{1}{a^{8}}\left(7+\frac{d_{0}\mu^{2}}{a_{M}^{3}}a^{5}-\frac{3d_{0}^{2}\mu^{4}}{16a_{M}^{6}}(a^{5})^{2}\right)<0$

(the associated degree two polynomial has a unique positive zero which is
greater than $a_{\max}$). Therefore, for $a_{\min}\leq a\leq a_{\max}=30000$
we have that

$0<\widetilde{\mathcal{M}}(\sigma;a_{\min})\leq\widetilde{\mathcal{M}}(\sigma;a)\leq\widetilde{\mathcal{M}}(\sigma;a_{\max}).$

We present now some numerical results. In Figure 5 it is depicted
$\lambda_{+}(\sigma;L)>0$ as a function of $\sigma$ for some values of $a$.

We can also compute numerically $\sigma_{\pm}$ and $\sigma_{*}$. We compute
$\sigma_{\pm}$ using that $\mathcal{X}_{\pm}(\sigma_{\pm};L)=0$ and
$\sigma_{*}$ as the unique positive zero of

$\beta^{2}(L)(1600\sigma^{3}+1920\sigma^{2}+432\sigma-32)+\widetilde{\beta}^{2}(L)(16\sigma(2\sigma+3)^{2}).$

We show its values in Table 2.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 5: These figures show $\lambda_{+}(\sigma;L)>0$ as functions of $\sigma$ for $a=\\{10000,15000,20000,24000,29600,30000\\}\,\text{km}$, when the origin is a saddle point (the units are the ones made explicit in Remark 3.1). $a$ | $\sigma_{-}(L)$ | $\sigma_{+}(L)$ | $\sigma_{*}(L)$
---|---|---|---
$10000$ | 0.058130693719535 | 0.058384404707503 | 0.058257569492108
$15000$ | 0.057294284878483 | 0.059218520762404 | 0.058257569455989
$20000$ | 0.054201084378365 | 0.062272995357407 | 0.058257569273126
$24000$ | 0.048177348389369 | 0.068087985206979 | 0.058257568831379
29600 | 0.029613649805289 | 0.084971418151141 | 0.058257567158286
$30000$ | 0.027639200647529 | 0.086680627913961 | 0.058257566962293
Table 2: Values of $\sigma_{\pm}$ and $\sigma_{*}$ for some choices of semi-
major axis.

######  Remark 4.5.

We emphasize that the values of $\sigma_{*}(L)$ are close to
$\frac{m}{2}=0.058257569495584$ which is the only positive zero of the
polynomial
$\mathbf{a}(\sigma)\mathbf{a}^{\prime}(\sigma)=1600\sigma^{3}+1920\sigma^{2}+432\sigma-32$.
This is because $\beta(L)\sim 6$ for all the values of semi-major axis $a$
considered, meanwhile the range of values for $\widetilde{\beta}(L)$ goes from
$\mathcal{O}(10^{-10})$ to $\mathcal{O}(10^{-7})$. Therefore, the zeros of
$\partial_{\sigma}\mathcal{M}$ are close to the ones of
$\mathbf{a}(\sigma)\mathbf{a}^{\prime}(\sigma)$.

To finish, in Table 3, we present the maximum value of the positive eigenvalue
in the saddle case, namely $\lambda_{+}(\sigma_{*}(L);L)$.

$a$ | $\lambda_{+}(\sigma_{*}(L);L)$
---|---
$10000$ | $0.97568857909377\cdot 10^{-6}$
$15000$ | $1.79245437499162\cdot 10^{-6}$
$20000$ | $2.75966404251617\cdot 10^{-6}$
$24000$ | $3.62767259360406\cdot 10^{-6}$
$29600$ | $4.96876878785117\cdot 10^{-6}$
$30000$ |  $5.06982657623624\cdot 10^{-6}$.
Table 3: Maximum value of the positive eigenvalue in the saddle case for some
values of the semi-major axis $a$.

##  5 Eccentric critical points of the $h$-averaged system

We devote this section to prove Theorem 3.5 and analyze numerically the
results obtained in this theorem. Along this section, we will use the notation
introduced in Section 4.1 without any explicit mention.

The first step to compute the location of the equilibrium points of the
$h$-averaged problem (4.1) is to analyze the equation $\dot{y}=0$. It implies
that a fixed point exists if it satisfies one of the following conditions:

$x=k\pi,\quad k=0,1\quad\mbox{ or }\quad\sqrt{P_{1}}Q_{1}=0.$ (5.1)

However, since we are assuming that $y\in(0,L/2)$ the latter condition cannot
happen. That is, all critical points of (4.1) must be of the form

$(k\pi,y)\quad\text{with}\quad y\in(0,L/2).$

###  5.1 Numerical approach

For given $L$ and $\Gamma$, equation (5.1) implies that an eccentric orbit can
be a critical point if it satisfies the condition $x=0,\pi$. The location of
the critical point (equivalently its associated eccentricity for a given
$\Gamma$) is obtained by solving the equations

$\begin{split}\sqrt{P_{1}(y,\Gamma)}{\mathcal{B}}_{0,\alpha}(y,\Gamma;L)+\alpha^{3}{\mathcal{B}}_{1}(y,\Gamma;L)&=0,\quad\mbox{
for $x=0$, }\quad\\\
\sqrt{P_{1}(y,\Gamma)}{\mathcal{B}}_{0,\alpha}(y,\Gamma;L)-\alpha^{3}{\mathcal{B}}_{1}(y,\Gamma;L)&=0,\quad\mbox{
for $x=\pi$, }\quad\end{split}$ (5.2)

where $P_{1}$, ${\mathcal{B}}_{0,\alpha}$ and ${\mathcal{B}}_{1}$ are
polynomials. In other words, they are roots of the following polynomial in
$(y,\Gamma)$

$\begin{split}P_{1}(y,\Gamma)({\mathcal{B}}_{0,\alpha}(y,\Gamma;L))^{2}-\alpha^{6}({\mathcal{B}}_{1}(y,\Gamma;L))^{2}=0.\end{split}$

A way to obtain the solutions is to fix $y$ such that $0<y/L<1/2$, compute the
roots of a polynomial of degree 6 in $\Gamma$ and select the ones for which
$0<\Gamma<L/2$ and (5.2) is satisfied. The results, computed for
$a=19000\,\mbox{km}$, $a=24000\,\mbox{km}$, $a=25450\,\mbox{km}$ and
$a=29600\,\mbox{km}$, are depicted in the left panel of the Figure 6 in the
action space $(\Gamma/L,y/L)$-plane.

The stability of the two families of equilibria can be computed from the
Hessian matrix of the Hamiltonian. More precisely, for fixed $L$ and $\Gamma$
and a given fixed point $(x_{0},y_{0})$, the variational equations read

$\begin{pmatrix}\dot{x}\\\
\dot{y}\end{pmatrix}=M^{\Gamma,L}(x_{0},y_{0})\begin{pmatrix}x\\\
y\end{pmatrix}\quad\mbox{ with }\quad
M^{\Gamma,L}=\begin{pmatrix}\frac{\partial^{2}\mathcal{H}^{\Gamma,L}}{\partial
y\partial x}&\frac{\partial^{2}\mathcal{H}^{\Gamma,L}}{\partial y^{2}}\\\
-\frac{\partial^{2}\mathcal{H}^{\Gamma,L}}{\partial
x^{2}}&-\frac{\partial^{2}\mathcal{H}^{\Gamma,L}}{\partial x\partial
y}\end{pmatrix}.$

Since $\mathrm{tr}{M^{\Gamma,L}}=0$, $M^{\Gamma,L}$ possesses two eigenvalues
$(\lambda,-\lambda)$ with either $\lambda\in i\mathbb{R}$ or
$\lambda\in\mathbb{R}$. Hence, the sign of $\det{M^{\Gamma,L}}=-\lambda^{2}$
characterizes the linear stability of the fixed point:

  * •

hyperbolic fixed point (saddle) for $\det{M^{\Gamma,L}}<0$;

  * •

elliptic fixed point (center) for $\det{M^{\Gamma,L}}>0$.

The results are depicted in the right panel of Figure 6.

Figure 6: Evolution of the fixed point location (left) and stability (right)
in the $(\Gamma/L,y/L)$-plane (Action space). In the plots in the first
column, the blue and red curves correspond respectively to the fixed point
families for $x=0$ and $x=\pi$. In the plots in the second column, the violet
and orange curves correspond respectively to the elliptic and hyperbolic part
of the $x=0$ and $x=\pi$ fixed point families.

###  5.2 Existence of eccentric critical points: Formulation of the problem

To prove Theorem 3.5, we have to analyze the existence and character of the
fixed points of the form $(0,y)$ and $(\pi,y)$ with $0<y<\frac{L}{2}$.

First, let us now rewrite condition (5.2) on the fixed points. Using the
decomposition in (4.4), the condition (5.2) can be rephrased as

$\mathcal{F}_{\pm}(y,\Gamma;L):=B_{0}(y,\Gamma;L)+\alpha^{3}B_{1}(y,\Gamma;L)\pm\alpha^{3}\frac{{\mathcal{B}}_{1}(y,\Gamma;L)}{\sqrt{P_{1}(y,\Gamma)}}=0$
(5.3)

where the sign $+$ corresponds to $x=0$ and the $-$ sign corresponds to
$x=\pi$.

We introduce the scaling

$y=L\hat{y},\qquad\Gamma=ys=L\hat{y}s.$

The values of the variables and parameters we are interested in are
$0<\hat{y}<\frac{1}{2}$, $a_{\min}\leq a\leq a_{\max}$ and
$s\in(0,\frac{1}{2})$

We rewrite $\mathcal{F}$ in these new variables $(s,\hat{y})$, namely
$\widehat{\mathcal{F}}_{\pm}(s,\hat{y};L):=\mathcal{F}_{\pm}(L\hat{y},L\hat{y}s;L)$.
Tedious but easy computations lead to

$\displaystyle\widehat{\mathcal{F}}_{\pm}(s,\hat{y};L)=$
$\displaystyle-6L^{3}\hat{y}^{2}(5s^{2}+8s-1)+12d_{0}\alpha^{3}L^{7}\hat{y}^{5}(4\hat{y}^{2}-5s-12s\hat{y}^{2}-5s^{2})$
$\displaystyle\pm
4d_{1}\alpha^{3}L^{7}\hat{y}^{5}\frac{(3+s)(12\hat{y}^{2}-8s\hat{y}^{2}-s^{2})}{\sqrt{3-2s-s^{2}}}$

where $d_{0}$ and $d_{1}$ were introduced in (4.2). We define

$\displaystyle\mathcal{G}_{\pm}(s,\hat{y};L):=$ $\displaystyle
5s^{2}+8s-1+2d_{0}\alpha^{3}L^{4}\hat{y}^{3}(5s^{2}+12s\hat{y}^{2}+5s-4\hat{y}^{2})$
$\displaystyle\pm
2d_{1}\alpha^{3}L^{4}\hat{y}^{3}\frac{(3+s)(s^{2}+8s\hat{y}^{2}-12\hat{y}^{2})}{3\sqrt{3-2s-s^{2}}}.$

Then, recalling that $y=L\hat{y}$ and $\Gamma=ys=L\hat{y}s$,

$\mathcal{F}_{\pm}(y,\Gamma;L)=0\Longleftrightarrow\widehat{\mathcal{F}}_{\pm}(s,\hat{y};L)=0\Longleftrightarrow\mathcal{G}_{\pm}(s,\hat{y};L)=0.$
(5.4)

The function $\mathcal{G}_{\pm}$ can be expressed as

$\mathcal{G}_{\pm}(s,\hat{y};L)=\mathbf{a}(s)+2\alpha^{3}L^{4}\rho[\mathbf{b}_{\pm}(s)\hat{y}^{3}+\mathbf{c}_{\pm}(s)\hat{y}^{5}]$
(5.5)

with

$\displaystyle\mathbf{a}(s)$ $\displaystyle=5{s}^{2}+8{s}-1$
$\displaystyle\mathbf{b}_{\pm}(s)$
$\displaystyle=10U_{2}^{0,0}{s}(1+{s})\pm\frac{20U_{2}^{1,0}(3+{s})}{3\sqrt{3-2{s}-{s}^{2}}}{s}^{2}$
$\displaystyle\mathbf{c}_{\pm}(s)$
$\displaystyle=8U_{2}^{0,0}(3{s}-1)\pm\frac{20U_{2}^{1,0}(3+{s})}{3\sqrt{3-2{s}-{s}^{2}}}(8{s}-12).$

######  Lemma 5.1.

The functions $\mathbf{a}$, $\mathbf{b}_{\pm}$, $\mathbf{c}_{\pm}$ are
strictly increasing functions at $\left[0,\frac{1}{2}\right]$. As a
consequence, for $\hat{y},L>0$, we have that
$\partial_{s}\mathcal{G}_{\pm}(s,\hat{y};L)>0$.

In addition, $\mathbf{a}(m)=0$, $\mathbf{a}(s)<0$ if $s\in[0,m)$ and
$\mathbf{a}(s)>0$ if $s>m$.

###### Proof.

The statements related to $\mathbf{a}$ are immediate. We compute now
$\partial_{s}\mathbf{b}_{\pm},\partial_{s}\mathbf{c}_{\pm}$,

$\displaystyle\partial_{s}\mathbf{b}_{\pm}(s)$
$\displaystyle=10U_{2}^{0,0}(1+2{s})\pm
20U_{2}^{1,0}\frac{2s(9-4s^{2}-s^{3})}{3(3-2{s}-{s}^{2})^{3/2}}$
$\displaystyle\partial_{s}\mathbf{c}_{\pm}(s)$ $\displaystyle=24U_{2}^{0,0}\mp
20U_{2}^{1,0}\frac{8s^{2}}{3(1-s)\sqrt{3-2{s}-{s}^{2}}}.$

The function $s(9-4s^{2}-s^{3})$ is an increasing function for
$s\in\left(0,0.7\right)$ provided

$9-12\cdot 0.7^{2}-4\cdot 0.7^{3}=1.748>0.$

Therefore, denoting $s_{*}=0.35$, if $s\in(0,s_{*}]$,

$\displaystyle\left|\pm
20U_{2}^{1,0}\frac{2s(9-4s^{2}-s^{3})}{3(3-2{s}-{s}^{2})^{3/2}}\right|$
$\displaystyle\leq
20U_{2}^{1,0}\frac{2s_{*}(9-4s_{*}^{2}-s_{*}^{3})}{3(3-2{s_{*}}-{s_{*}}^{2})^{3/2}}$
$\displaystyle=6.731992451918273.$

Since $10U_{2}^{0,0}=7.62646$, we conclude that
$\partial_{s}\mathbf{b}_{\pm}(s)>0$ if $s\in(0,s_{*}]$. When
$s\in\left(s_{*},\frac{1}{2}\right]$,

$\left|\pm
20U_{2}^{1,0}\frac{2s(9-4s^{2}-s^{3})}{3(3-2{s}-{s}^{2})^{3/2}}\right|\leq
12.414817621987043$

and $10U_{2}^{0,0}(1+2s)\geq 10U_{2}^{0,0}(1+2s_{*})=12.96498200000000$.
Therefore, we conclude that $\partial_{s}\mathbf{b}_{\pm}(s)>0$ if
$s\in\left[0,\frac{1}{2}\right]$. On the other hand, if $0<s\leq\frac{1}{2}$,
then

$\left|\mp 20U_{2}^{1,0}\frac{8s^{2}}{3(1-s)\sqrt{3-2{s}-{s}^{2}}}\right|\leq
20U_{2}^{1,0}\frac{2}{\frac{3}{2}\sqrt{\frac{7}{4}}}=11.035393441766260$

and $24U_{2}^{0,0}=18.303504$. Therefore $\partial_{s}\mathbf{c}_{\pm}(s)>0$.

When $s\geq 0$, $\partial_{s}\mathbf{a}(s)=8+10s>0$ and thus
$\partial_{s}\mathcal{G}_{\pm}(s,\hat{y};L)>0$ provided $\hat{y},L>0$. ∎

###  5.3 A preliminary existence result

The next step to compute the eccentric critical points is to solve the
equations in (5.4). We look for zeros of the functions $\mathcal{G}_{\pm}$ in
$s$.

######  Lemma 5.2.

Let $L$ be satisfying $0<L<L_{\max}$, $L=\sqrt{\mu a}$, and

$\hat{y}_{\min}=\sqrt{-\frac{\mathbf{b}_{+}(m)}{\mathbf{c}_{+}(m)}}=0.116589071022807.$
(5.6)

We have that $\mathcal{G}_{+}(m,\hat{y}_{\min};L)=0$ and

  1. 1.

For any $0<\hat{y}\leq\frac{1}{2}$, there exists a unique analytic function
$s_{+}(\hat{y};L)$, which belongs to $\left(0,\frac{3}{2}m\right)$, satisfying
$\mathcal{G}_{+}(s_{+}(\hat{y};L),\hat{y};L)=0$. In addition, for
$\hat{y}_{\min}<\hat{y}\leq\frac{1}{2}$ the function $s_{+}(\hat{y};L)$
belongs to $\left(m,\frac{3}{2}m\right)$ and when $0<\hat{y}<\hat{y}_{\min}$,
$s_{+}(\hat{y};L)\in(0,m)$.

  2. 2.

For any $0<\hat{y}\leq\frac{1}{2}$, there exists a unique analytic function
$s_{-}(\hat{y};L)$, which belongs to $(\frac{2}{5}m,m)$, satisfying
$\mathcal{G}_{-}(s_{-}(\hat{y};L),\hat{y};L)=0$.

###### Proof.

To prove the existence of solutions of the equation
$\mathcal{G}_{\pm}(s,\hat{y};L)=0$, we fix $\hat{y},L$ (or $a$) in the ranges
$0<\hat{y}\leq\frac{1}{2}$ and $0\leq a\leq a_{\max}$ and we use a Bolzano
argument (with respect to the $s$ variable).

We deal first with the $+$ case. Since $\mathbf{b}_{+}(0)=0$,

$\mathcal{G}_{+}(0,\hat{y};L)=-1+2\alpha^{3}L^{4}\hat{y}^{5}\mathbf{c}_{+}(0).$

It is clear from the definition of $\mathbf{c}_{+}$ that $\mathbf{c}_{+}(0)<0$
so that $\mathcal{G}_{+}(0,\hat{y};L)<0$. For fixed $L$, now we prove that the
function $\mathcal{G}_{+}\left(\frac{3}{2}m,\hat{y};L\right)$ is positive. We
introduce

$C_{+}(\hat{y}):=\mathcal{G}_{+}\left(\frac{3}{2}m,\hat{y};L\right)=\hat{\mathbf{a}}+2\alpha^{3}L^{4}\rho\left[\hat{\mathbf{b}}_{+}\hat{y}^{3}+\hat{\mathbf{c}}_{+}\hat{y}^{5}\right]$

with

$\widehat{\mathbf{a}}=\mathbf{a}\left(\frac{3}{2}m\right)=0.550909166052992,\quad\widehat{\mathbf{b}}_{+}=\mathbf{b}_{+}\left(\frac{3}{2}m\right)=1.784508217583372,$
(5.7)
$\widehat{\mathbf{c}}_{+}=\mathbf{c}_{+}\left(\frac{3}{2}m\right)=-78.794344798569341.$
(5.8)

For a fixed $L$, the function $C_{+}(\hat{y})$ has the minimum value either at
the points
$\hat{y}=0,\hat{y}=\sqrt{\frac{-3\widehat{\mathbf{b}}_{+}}{5\widehat{\mathbf{c}}_{+}}}=0.116570155823739$
or $\hat{y}=\frac{1}{2}$. Evaluating at these points, we obtain that
$C_{+}(0)=\widehat{\mathbf{a}}>0$, from the fact that
$\widehat{\mathbf{c}}_{+}<0$:

$C_{+}\left(\sqrt{\frac{-3\widehat{\mathbf{b}}_{+}}{5\widehat{\mathbf{c}}_{+}}}\right)=\widehat{\mathbf{a}}-2\alpha^{3}L^{4}\rho\frac{6\widehat{\mathbf{b}}_{+}^{2}}{25\widehat{\mathbf{c}}_{+}}\sqrt{\frac{-3\widehat{\mathbf{b}}_{+}}{5\widehat{\mathbf{c}}_{+}}}>0$

and finally, using that
$4\widehat{\mathbf{b}}_{+}+\widehat{\mathbf{c}}_{+}<0$,

$\displaystyle C_{+}\left(\frac{1}{2}\right)$
$\displaystyle=\widehat{\mathbf{a}}+\alpha^{3}L^{4}\rho\frac{1}{16}(4\widehat{\mathbf{b}}_{+}+\widehat{\mathbf{c}}_{+})\geq\widehat{\mathbf{a}}+\rho\frac{a_{\max}^{5}\mu^{2}}{16a_{M}^{3}}(4\widehat{\mathbf{b}}_{+}+\widehat{\mathbf{c}}_{+})$
(5.9) $\displaystyle=0.014481715549738>0.$

The previous analysis proves that

$\mathcal{G}_{+}(0,\hat{y};L)\cdot\mathcal{G}_{+}\left(\frac{3}{2}m,\hat{y};L\right)<0$

and therefore, there exists
$s_{+}=s_{+}(\hat{y};L)\in\left(0,\frac{3}{2}m\right)$ such that
$\mathcal{G}_{+}(s_{+},\hat{y};L)=0$.

By Lemma 5.1, $\partial_{s}\mathcal{G}_{+}(s,\hat{y};L)>0$ so that, for fixed
$\hat{y},L$, $s_{+}=s_{+}(\hat{y};L)$ is the unique solution of
$\mathcal{G}_{+}(s,\hat{y};L)=0$. By the implicit function theorem, the
dependence on $\hat{y},L$ of $s_{+}(\hat{y};L)$ is analytic.

To finish the result related to the $+$ case, we consider now

$\mathcal{G}_{+}(m,\hat{y};L)=2\alpha^{3}L^{4}\rho\left[\mathbf{b}_{+}(m)\hat{y}^{3}+\mathbf{c}_{+}(m)\hat{y}^{5}\right].$

Notice that its sign does not depend on $L$. We have that

$\mathbf{b}_{+}(m)=1.085189654741836,\qquad\mathbf{c}_{+}(m)=-79.834380790596938.$

Therefore, for $0<\hat{y}<\hat{y}_{\min}$, $\mathcal{G}_{+}(m,\hat{y};L)>0$
and when $\hat{y}>\hat{y}_{\min}$, then $\mathcal{G}_{+}(m,\hat{y};L)<0$. We
conclude that for any $0\leq a\leq a_{\max}$ and
$\hat{y}_{\min}<\hat{y}\leq\frac{1}{2}$ the function,
$s_{+}(\hat{y};L)\in(m,\frac{3}{2}m)$ and when $0<\hat{y}<\hat{y}_{\min}$, we
have that $s_{+}(\hat{y};L)\in(0,m)$.

For the $-$ case, we proceed analogously. We study the function
$\mathcal{G}_{-}(m,\hat{y};L)$. Using that $\mathbf{a}(m)=0$,

$\mathcal{G}_{-}(m,\hat{y};L)=2\alpha^{3}L^{4}\rho\left[{\mathbf{b}}_{-}(m)\hat{y}^{3}+{\mathbf{c}}_{-}(m)\hat{y}^{5}\right]$

with

$\mathbf{b}_{-}(m)=0.899076688965464,\qquad\mathbf{c}_{-}(m)=71.897315415767750,$

and hence $\mathcal{G}_{-}(m,\hat{y};L)>0$ provided $\hat{y},L>0$. Moreover,
we have that

$\mathbf{a}\left(\frac{2}{5}m\right)=-0.616290933136957,\qquad\mathbf{b}_{-}\left(\frac{2}{5}m\right)=0.357833740529904$
(5.10)

and

$\mathbf{c}_{-}\left(\frac{2}{5}m\right)=70.608271414680971.$ (5.11)

From these values we have that

$\mathcal{G}_{-}\left(\frac{2}{5}m,\hat{y};L\right)\leq\mathcal{G}_{-}\left(\frac{2}{5}m,\frac{1}{2};L_{\max}\right)=-0.076994095997652<0.$

As a consequence, for any $0<\hat{y}\leq\frac{1}{2}$ and $0<a\leq a_{\max}$,
there exists $s_{-}=s_{-}(\hat{y};L)\in(\frac{2}{5}m,m)$ such that
$\mathcal{G}_{-}(s_{-},\hat{y};L)=0$. Then, for fixed $\hat{y},L$, since by
Lemma 5.1, $\partial_{s}\mathcal{G}_{-}>0$, $s_{-}(\hat{y};L)$ is the unique
zero of $\mathcal{G}_{-}(s,\hat{y};L)$ and by the Implicit Function Theorem
$s_{-}$ is an analytic function. ∎

######  Remark 5.3.

We stress that if the semi-major axis $a\to\infty$ (that is $L\to\infty$),
$C_{+}\left(\frac{1}{2}\right)$ in (5.9) is negative provided
$4\widehat{\mathbf{b}}_{+}+\widehat{\mathbf{c}}_{+}<0$. In order to be a
negative quantity, the maximum semi-major axis has to satisfy

$\alpha^{3}L^{4}<\frac{16\widehat{\mathbf{a}}}{|4\widehat{\mathbf{b}}_{+}+\widehat{\mathbf{c}}_{+}|\rho}\Leftrightarrow
a<\left(\frac{16\hat{\mathbf{a}}\
a_{M}^{3}}{|4\hat{\mathbf{b}}_{+}+\hat{\mathbf{c}}_{+}|\rho\mu^{2}}\right)^{1/5}=30160.25822948035\,\mathrm{km},$

that corresponds to $a<1.018927642887850$ in the units specified in Remark
3.1. This fact leads us to think that $a_{\max}$ is (almost) the optimal value
for our arguments to be true.

######  Remark 5.4.

We have that $\mathbf{b}_{-}(\frac{1}{2}m),\mathbf{c}_{-}(\frac{1}{2}m)>0$.
Then for semi-major axis satisfying

$a\leq\left(-\frac{16\mathbf{a}(\frac{1}{2}m)a_{M}^{3}}{\mu^{2}\rho(4\mathbf{b}_{-}(\frac{1}{2}m)+\mathbf{c}_{-}(\frac{1}{2}m))}\right)^{\frac{1}{5}}=29700.20301662558\,\mathrm{km}$

(corresponding to $a\leq 1.003385237048162$ in the units in Remark 3.1) we
have that $\mathcal{G}_{-}(\frac{m}{2},\hat{y};L)<0$ and therefore
$s_{-}(\hat{y};L)\in\left(\frac{m}{2},m\right)$.

A straightforward corollary of this preliminary analysis is the following
result, see (5.4).

######  Corollary 5.5.

Let $a$ be a semi-major axis satisfying $0<a\leq a_{\max}$. Then, given
$L\hat{y}_{\min}<y\leq\frac{L}{2}$ with $L=\sqrt{\mu a}\in(0,L_{\max}]$, there
exist two analytic functions $\Gamma_{\pm}(y;L)$ such that

$\frac{\Gamma_{+}(y;L)}{my}\in\left(1,\frac{3}{2}\right),\;\frac{\Gamma_{-}(y;L)}{my}\in\left(\frac{2}{5},1\right),\qquad\mathcal{F}_{\pm}(y,\Gamma_{\pm}(y;L);L)=0.$

For $y,L$ fixed, $\Gamma_{\pm}:=\Gamma_{\pm}(y;L)$ are the unique solutions of
$\mathcal{F}_{\pm}(y,\Gamma_{\pm};L)=0$ belonging to
$\left(0,\frac{y}{2}\right)$. We also have that
$\partial_{\Gamma}\mathcal{F}_{\pm}(y,\Gamma;L)<0$ for
$\Gamma\in\left(0,\frac{y}{2}\right)$.

###### Proof.

For the existence and regularity of $\Gamma_{\pm}$, we only need to define
$\Gamma_{\pm}(y;L)=ys_{\pm}\left(\frac{y}{L};L\right)$ and to apply Lemma 5.1.
In addition, using (5.4), we have that

$\mathcal{F}_{\pm}(y,\Gamma;L)=\hat{\mathcal{F}}_{\pm}\left(\frac{\Gamma}{y},\frac{y}{L};L\right)=-6Ly^{3}\mathcal{G}_{\pm}\left(\frac{\Gamma}{y},\frac{y}{L};L\right)$

and then

$\partial_{\Gamma}\mathcal{F}_{\pm}(y,\Gamma;L)=-6Ly^{2}\partial_{s}\mathcal{G}_{\pm}\left(\frac{\Gamma}{y},\frac{y}{L};L\right)<0,$

where we have used that, by Lemma 5.1,
$\partial_{s}\mathcal{G}_{\pm}(s,\hat{y};L)>0$. ∎

Corollary 5.5 proves the existence of critical points. However, we can not
elucidate neither the range of the parameters $\Gamma,L$ for which they exist
nor how many critical points the system possesses for given values of the
parameters. To do so, let us study the functions

$\hat{\Gamma}_{\pm}(\hat{y};L)=\frac{\Gamma_{\pm}(L\hat{y};L)}{L}=\hat{y}s_{\pm}(\hat{y};L),$

where $s_{\pm}(\hat{y};L)$ are defined implicitly by the equations
$\mathcal{G}_{\pm}(s,\hat{y};L)=0$ (see Lemma 5.2). That is
$\mathcal{G}_{\pm}(s_{\pm}(\hat{y};L),\hat{y};L)=0$. Therefore, writing
$s_{\pm}=s_{\pm}(\hat{y};L)$ and taking the derivative with respect to
$\hat{y}$ in the equation $\mathcal{G}_{\pm}(s,\hat{y};L)=0$, we have that

$\partial_{\hat{y}}\mathcal{G}_{\pm}(s_{\pm},\hat{y};L)+\partial_{s}\mathcal{G}_{\pm}(s_{\pm},\hat{y};L)\partial_{\hat{y}}s_{\pm}(\hat{y};L)=0.$

We recall that by Lemma 5.1, $\partial_{s}\mathcal{G}_{\pm}>0$. This fact
allows computing $\partial_{\hat{y}}s_{\pm}(\hat{y};L)$ in terms of $s_{\pm}$
and, consequently,

$\displaystyle\partial_{\hat{y}}\hat{\Gamma}_{\pm}(\hat{y};L)$
$\displaystyle=s_{\pm}(\hat{y};L)+\hat{y}\partial_{\hat{y}}s_{\pm}(\hat{y};L)$
$\displaystyle=\frac{1}{\partial_{s}\mathcal{G}_{\pm}(s_{+},\hat{y};L)}\left[s_{\pm}(\hat{y};L)\partial_{s}\mathcal{G}_{\pm}(s_{\pm},\hat{y};L)-\hat{y}\partial_{\hat{y}}\mathcal{G}_{\pm}(s_{\pm},\hat{y};L)\right].$
(5.12)

We want to study the possible changes in the sign of
$\partial_{\hat{y}}\hat{\Gamma}_{\pm}$. We notice that, following the same
kind of computations,

$\partial_{L}\hat{\Gamma}_{\pm}(\hat{y};L)=-\frac{\partial_{L}\mathcal{G}_{\pm}(s_{\pm},\hat{y};L)}{\partial_{s}\mathcal{G}_{\pm}(s_{+},\hat{y};L)}\hat{y}.$
(5.13)

######  Remark 5.6.

Notice that the value $\hat{y}_{\min}$ defined in Lemma 5.2 satisfies that
$\hat{y}_{\min}<\hat{y}_{\mathrm{col}}$, the collision value
$\hat{y}_{\mathrm{col}}$ in (3.8). For that reason, from now on, we will
restrict our analysis to values of
$\hat{y}\in\left[\hat{y}_{\min},\frac{1}{2}\right)$.

As we will see below, it turns out that the analysis of the functions
$\hat{\Gamma}_{+}$ and $\hat{\Gamma}_{-}$ are quite different, being the
corresponding to $\hat{\Gamma}_{-}$ more involved.

We start with the result for the fixed points of the form $(0,y)$, which
corresponds with the $+$ case, in Section 5.4, and we postpone the study of
the fixed points of the form $(\pi,y)$ to Section 5.5.

###  5.4 The fixed points of the form $(0,y)$

For any $L_{\min}\leq L\leq L_{\max}$ we define

$\displaystyle\hat{\Gamma}_{\min}^{+}$
$\displaystyle=\hat{y}_{\min}m=0.013584391815073,$
$\displaystyle\hat{\Gamma}_{\max}^{+}(L)$
$\displaystyle=\hat{\Gamma}_{+}\left(\frac{1}{2};L\right)=\frac{1}{2}s_{+}\left(\frac{1}{2};L\right).$

We emphasize that, by Lemma 5.2,
$s_{+}\left(\frac{1}{2};L\right)\in\left(m,\frac{3}{2}m\right)$ and therefore

$\hat{\Gamma}_{\max}^{+}(L)\in\left[\frac{1}{2}m,\frac{3}{4}m\right]=[0.058257569495584,0.087386354243376].$
(5.14)

######  Proposition 5.7.

The function $\hat{\Gamma}_{\max}^{+}(L)$ is strictly increasing on
$[L_{\min},L_{\max}]$ and therefore

$\hat{\Gamma}_{\max}^{+}(L)\in[\hat{\Gamma}_{\max}^{+}(L_{\min}),\hat{\Gamma}_{\max}^{+}(L_{\max})]=[0.058270961487710,0.086680627913961]$

where these values have been computed numerically222 Notice that there is not
a significant difference with the estimates in (5.14).. Moreover,

$\hat{\Gamma}^{+}(\hat{y};L)\in[\hat{\Gamma}_{\min}^{+},\hat{\Gamma}_{\max}^{+}(L)],\qquad\text{for
all }\;\hat{y}_{\min}\leq\hat{y}\leq\frac{1}{2},$

and $\hat{\Gamma}^{+}(\cdot;L)$ is an injective increasing function.

With respect to the equilibrium points, for any $L\in[L_{\min},L_{\max}]$:

  1. 1.

If $\hat{\Gamma}\notin[\hat{\Gamma}_{\min}^{+},\hat{\Gamma}_{\max}^{+}(L)]$,
system (4.1) has no fixed points of the form $(0,L\hat{y})$ satisfying
$\hat{y}_{\min}<\hat{y}\leq\frac{1}{2}$ and
$\hat{\Gamma}\in\left(0,\frac{1}{2}\hat{y}\right)$.

  2. 2.

If $\hat{\Gamma}\in[\hat{\Gamma}_{\min}^{+},\hat{\Gamma}_{\max}^{+}(L)]$,
there exists a unique $\hat{y}_{+}=\hat{y}_{+}(\hat{\Gamma};L)$ such that
$(0,L\hat{y}_{+})$ is a fixed point of the system (4.1). In addition,
$\hat{y}_{+}(\cdot;L)$ is a strictly increasing function.

###### Proof.

We start with the statement related to $\hat{\Gamma}_{\max}^{+}$. Using
(5.13), we have that

$\partial_{L}\hat{\Gamma}_{\max}^{+}(L)=\partial_{L}\hat{\Gamma}^{+}\left(\frac{1}{2};L\right)=-\frac{\partial_{L}\mathcal{G}_{+}\left(s_{+},\frac{1}{2};L\right)}{2\partial_{s}\mathcal{G}_{+}\left(s_{+},\frac{1}{2};L\right)}.$

The sign of $\partial_{L}\mathcal{G}_{+}\left(s_{+},\frac{1}{2};L\right)$ is
the sign of $\mathbf{b}_{+}(s_{+})+\frac{1}{4}\mathbf{c}_{+}(s_{+})$ (see
definition (5.5)). Then, since by Lemma 5.1, $\mathbf{b}_{+},\mathbf{c}_{+}$
are increasing functions and by Lemma 5.2,
$s_{+}\in\left(m,\frac{3}{2}m\right)$

$\mathbf{b}_{+}(s_{+})+\frac{1}{4}\mathbf{c}_{+}(s_{+})\leq\mathbf{b}_{+}\left(\frac{3}{2}m\right)+\frac{1}{4}\mathbf{c}_{+}\left(\frac{3}{2}m\right)<0,$

where we have used (5.7) and (5.8) to estimate
$\mathbf{b}_{+}\left(\frac{3}{2}m\right)+\frac{1}{4}\mathbf{c}_{+}\left(\frac{3}{2}m\right)$.
Therefore, $\partial_{L}\mathcal{G}_{+}<0$ and, using that
$\partial_{s}\mathcal{G}_{+}>0$ by Lemma 5.1,
$\partial_{L}\Gamma_{\max}^{+}>0$.

Now we prove the rest of the properties. We write

$\partial_{\hat{y}}\mathcal{G}_{+}(s,\hat{y};L)=2\alpha^{3}L^{4}\rho\left[3\mathbf{b}_{+}(s)\hat{y}^{2}+5\mathbf{c}_{+}(s)\hat{y}^{4}\right].$

By Lemma 5.1, $\mathbf{b}_{+},\mathbf{c}_{+}$ are increasing functions and, by
Lemma 5.2, $s_{+}\in\left(m,\frac{3}{2}m\right)$. Therefore, for
$s\in\left(m,\frac{3}{2}m\right)$,

$\partial_{\hat{y}}\mathcal{G}_{+}(s,\hat{y};L)\leq
2\alpha^{3}L^{4}\rho\left[3\mathbf{b}_{+}\left(\frac{3}{2}m\right)\hat{y}^{2}+5\mathbf{c}_{+}\left(\frac{3}{2}m\right)\hat{y}^{4}\right].$

Using the values of
$\mathbf{b}_{+}\left(\frac{3}{2}m\right),\mathbf{c}_{+}\left(\frac{3}{2}m\right)$
in (5.7) and (5.8), if

$\hat{y}\geq\sqrt{-\frac{3\mathbf{b}_{+}\left(\frac{3}{2}m\right)}{5\mathbf{c}_{+}\left(\frac{3}{2}m\right)}}=0.116570155823739,$

we have then $\partial_{\hat{y}}\mathcal{G}_{+}(s,\hat{y};L)\leq 0$. In
particular, the same happens if $\hat{y}\geq\hat{y}_{\min}$ (see Lemma 5.2 for
the exact value of $\hat{y}_{\min}$). As a consequence, from (5.3),
$\partial_{\hat{y}}\hat{\Gamma}_{+}(\hat{y};L)>0$ if
$\hat{y}_{\min}\leq\hat{y}\leq\frac{1}{2}$ and $L_{\min}\leq L\leq L_{\max}$
so that $\hat{\Gamma}_{+}(\cdot;L)$ is an injective (strictly increasing)
function for any fixed $L\in[L_{\min},L_{\max}]$. The range of values of
$\hat{\Gamma}_{+}$ for a given $L$ is then

$\hat{\Gamma}^{+}(\hat{y};L)\in\left[\hat{\Gamma}_{+}(\hat{y}_{\min};L),\hat{\Gamma}_{+}\left(\frac{1}{2};L\right)\right].$

We notice that, since $\mathcal{G}_{+}(m,\hat{y}_{\min};L)=0$ (see Lemma 5.1),

$\hat{\Gamma}_{+}(\hat{y}_{\min};L)=\hat{y}_{\min}s_{+}(\hat{y}_{\min};L)=\hat{y}_{\min}m=\hat{\Gamma}_{\min}^{+}.$

For a given $L$, let $\hat{y}_{+}(\hat{\Gamma};L)$ be such that

$\hat{\Gamma}_{+}(\hat{y}_{+}(\hat{\Gamma};L);L)=\hat{\Gamma}.$

It is clear that $\hat{y}_{+}(\cdot;L)$ is an increasing function and it is
defined for $\hat{\Gamma}\in[\hat{\Gamma}_{\min},\hat{\Gamma}_{\max}^{+}(L)]$.
∎

######  Remark 5.8.

If we want to control the range of $\hat{\Gamma}$ for
$\hat{y}_{\mathrm{col}}\leq\hat{y}\leq\frac{1}{2}$, we notice that, since
$\hat{\Gamma}_{+}(\cdot;L)$ is an increasing function, the “new”
$\hat{\Gamma}_{\min}^{+}=\hat{\Gamma}_{\min}^{+}(L)=\hat{\Gamma}_{+}(y_{\mathrm{col}};L)$.
Then, since
$\hat{\Gamma}_{+}(\hat{y}_{\mathrm{col}};L)=\hat{y}_{\mathrm{col}}s_{+}(\hat{y}_{\mathrm{col}};L)$,

$\hat{\Gamma}_{\min}^{+}(L)\in\left[\hat{y}_{\mathrm{col}}m,\hat{y}_{\mathrm{col}}\frac{3}{2}m\right]=[0.035912783574701,0.053869175362051].$

We emphasize that $\hat{\Gamma}_{+}(\hat{y}_{\mathrm{col}};\cdot)$ is also an
increasing function. Indeed, writing $s_{+}=s_{+}(\hat{y}_{\mathrm{col}};L)$,
the sign of $-\partial_{L}\hat{\Gamma}_{+}(\hat{y}_{\mathrm{col}};\cdot)$ is
the sign of

$\mathbf{b}_{+}(s_{+})+\hat{y}_{\mathrm{col}}^{2}\mathbf{c}_{+}(s_{+})\leq\mathbf{b}_{+}\left(\frac{3}{2}m\right)+\hat{y}_{\mathrm{col}}^{2}\mathbf{c}_{+}\left(\frac{3}{2}m\right)=-5.701123730321832<0$

where, again, we have used (5.7) and (5.8) for
$\mathbf{b}_{+}\left(\frac{3}{2}m\right)$ and
$\mathbf{c}_{+}\left(\frac{3}{2}m\right)$. Therefore, we numerically obtain
that

$\hat{\Gamma}_{\min}^{+}(L)\in[\hat{\Gamma}_{+}(y_{\mathrm{col}};L_{\min}),\hat{\Gamma}_{+}(y_{\mathrm{col}};L_{\max})]=[0.035913449550478,0.037427016836718].$

We notice that in this case, the numerical computation induces a more accurate
range of values of $\hat{\Gamma}_{\min}^{+}(L)$.

###  5.5 The fixed points of the form $(\pi,y)$

Now we pay attention to the $-$ case. In this case, $\hat{\Gamma}_{-}$ is not
an injective function and for that reason, for studying the behavior of
$\hat{\Gamma}_{-}$, we must control the existence of critical points, namely
the values of $\hat{y},L$ such that
$\partial_{\hat{y}}\hat{\Gamma}_{-}(\hat{y};L)=0$. To this end, we introduce

$\mathcal{H}(s,\hat{y};L):=s\partial_{s}\mathcal{G}_{-}(s,\hat{y};L)-\hat{y}\partial_{\hat{y}}\mathcal{G}_{-}(s,\hat{y};L),\qquad\partial_{\hat{y}}\hat{\Gamma}_{-}(\hat{y};L)=\frac{\mathcal{H}(s_{-},\hat{y};L)}{\partial_{s}\mathcal{G}_{-}(s_{-},\hat{y};L)}$
(5.15)

with $s_{-}=s_{-}(\hat{y};L)$. By the uniqueness statement in Lemma 5.1, we
have that $\mathcal{G}_{-}(s,\hat{y};L)=0$ if and only if
$s=s_{-}(\hat{y};L)$. Moreover
$s_{-}(\hat{y};L)\in\left(\frac{2}{5}m,m\right)$. Therefore, the equation

$\partial_{\hat{y}}\hat{\Gamma}_{-}(\hat{y};L)=0$

is equivalent to the existence of solutions of

$\mathcal{H}(s,\hat{y};L)=0,\qquad\mathcal{G}_{-}(s,\hat{y};L)=0$

under the restrictions $s\in\left(\frac{2}{5}m,m\right)$,
$\hat{y}\in\left[\hat{y}_{\min},\frac{1}{2}\right]$ (see Remark 5.6),
$L\in[L_{\min},L_{\max}]$.

From expression (5.5) of $\mathcal{G}_{-}$, we write the system in a more
suitable way

$\displaystyle\mathcal{H}(s,\hat{y};L)=\overline{\mathbf{a}}(s)+2\alpha^{3}L^{4}\rho\big{[}\overline{\mathbf{b}}(s)\hat{y}^{3}+\overline{\mathbf{c}}(s)\hat{y}^{5}]=0$
(5.16)
$\displaystyle\mathcal{G}_{-}(s,\hat{y};L)=\mathbf{a}(s)+2\alpha^{3}L^{4}\rho\big{[}\mathbf{b}_{-}(s)\hat{y}^{3}+\mathbf{c}_{-}(s)\hat{y}^{5}\big{]}=0$

with

$\displaystyle\overline{\mathbf{a}}(s)$
$\displaystyle=s\partial_{s}\mathbf{a}(s)=2s(5s+4)$ (5.17)
$\displaystyle\overline{\mathbf{b}}(s)$
$\displaystyle=s\partial_{s}\mathbf{b}_{-}(s)-3\mathbf{b}_{-}(s)=-10U_{2}^{0,0}s(s+2)-\frac{20U_{2}^{1,0}s^{2}(s^{2}+4s-3)}{3(1-s)\sqrt{3-2s-s^{2}}}$
$\displaystyle\overline{\mathbf{c}}(s)$
$\displaystyle=s\partial_{s}\mathbf{c}_{-}(s)-5\mathbf{c}_{-}(s)=8U_{2}^{0,0}(5-12s)-\frac{80U_{2}^{1,0}(8s^{3}+5s^{2}-60s+45)}{3(1-s)\sqrt{3-2s-s^{2}}}.$

Hence, in order to solve system (5.16), we write from the second equation

$2\alpha^{3}L^{4}\rho=-\frac{\mathbf{a}(s)}{\mathbf{b}_{-}(s)\hat{y}^{3}+\mathbf{c}_{-}(s)\hat{y}^{5}}$

and substituting this value into the first equation, we obtain that

$\mathcal{H}(s,\hat{y};L)=\frac{\overline{\mathbf{a}}(s)\mathbf{b}_{-}(s)-\mathbf{a}(s)\overline{\mathbf{b}}(s)+\hat{y}^{2}\big{(}\overline{\mathbf{a}}(s)\mathbf{c}_{-}(s)-\mathbf{a}(s)\overline{\mathbf{c}}(s)\big{)}}{\mathbf{b}_{-}(s)+\mathbf{c}_{-}(s)\hat{y}^{2}}=0$

so that

$\hat{y}^{2}=\frac{\mathbf{a}(s)\overline{\mathbf{b}}(s)-\overline{\mathbf{a}}(s)\mathbf{b}_{-}(s)}{\overline{\mathbf{a}}(s)\mathbf{c}_{-}(s)-\mathbf{a}(s)\overline{\mathbf{c}}(s)}.$

Therefore system (5.16) has a unique solution for those values of
$s\in\left(\frac{2}{5}m,m\right)$ such that the following restrictions hold

$\displaystyle{\hat{y}_{\min}}\leq\hat{y}$
$\displaystyle=\hat{y}_{*}(s):=\left(\frac{\mathbf{a}(s)\overline{\mathbf{b}}(s)-\overline{\mathbf{a}}(s)\mathbf{b}_{-}(s)}{\overline{\mathbf{a}}(s)\mathbf{c}_{-}(s)-\mathbf{a}(s)\overline{\mathbf{c}}(s)}\right)^{\frac{1}{2}}\leq\frac{1}{2},$
(5.18) $\displaystyle a_{\min}\leq a$
$\displaystyle=a_{*}(s):=\left(-\frac{\mathbf{a}(s)a_{M}^{3}}{2\mu^{2}\rho(\mathbf{b}_{-}(s)\hat{y}_{*}^{3}(s)+\mathbf{c}_{-}(s)\hat{y}_{*}^{5}(s))}\right)^{\frac{1}{5}}\leq
a_{\max}.$

######  Lemma 5.9.

For those values of $s$ such that
$\hat{y}_{*}(s)\in\left(\hat{y}_{\min},\frac{1}{2}\right]$, $\hat{y}_{*}(s)$
is an increasing function and $a_{*}(s)$ is decreasing. Defining implicitly
$s_{\min}$ and $s_{\max}$ such that

$a_{*}(s_{\min})=a_{\max},\qquad\hat{y}_{*}(s_{\max})=\frac{1}{2},$

we have that, for $s\in[s_{\min},s_{\max}]$, $\hat{\Gamma}_{-}(\hat{y};L)$ has
a unique critical point at $(\hat{y},L)=(\hat{y}_{*}(s),L_{*}(s))$ with
$L_{*}(s)=\sqrt{\mu a_{*}(s)}$. The values of $s_{\min}$ and $s_{\max}$ and
$L_{*}(s_{\max}),\hat{y}_{*}(s_{\min})$ can be computed numerically:

${s_{\min}=0.096577225237580,}\qquad s_{\max}=0.096796816334740,$

and then $a_{*}(s_{\max})=23893.56218133389\,\mathrm{km}$,

$L_{*}(s_{\max})=97590.90325766560\,\mathrm{km^{2}/s},\qquad\hat{y}_{*}(s_{\min})=0.397273020602216.$

In addition, if $L\in[L_{\min},L_{*}(s_{\max})]$ there are no critical points
of $\hat{\Gamma}_{-}(\cdot;L)$ belonging to
$\left[\hat{y}_{\min},\frac{1}{2}\right)$.

###### Proof.

One can perform an analytic thorough study of $\hat{y}_{*}(s),a_{*}(s)$ as a
function of $s$. However, we have preferred just to draw the explicit
functions $\hat{y}_{*}(s),a_{*}(s)$ in order to convince the reader about this
result. For those values of $s$ such that $\hat{y}_{*}(s)\in\mathbb{R}$ with
$0<\hat{y}_{*}(s)\leq\frac{1}{2}$, we draw the corresponding semi-major axis
$a$, see Figure 7.

|  .
---|---
Figure 7: On the left, $\hat{y}_{*}(s)$ and on the right $a_{*}(s)$. Only the
values of $s$ such that $\hat{y}_{*}(s)\in\left[0,\frac{1}{2}\right]$ are
considered. In this figure, the semi-major axis $a$ is measured in km.

These figures illustrate that indeed $\hat{y}_{*}(s)$ is an increasing
function of $s$ meanwhile $a_{*}(s)$ is decreasing.

If we only consider the values of $s$ such that
$\hat{y}_{*}(s)\in\left[\hat{y}_{\min},\frac{1}{2}\right]$ and
$a_{*}(s)\in[a_{\min},a_{\max}]$ (which are the ones we are interested in) we
observe the behaviour in Figure 8.

|  .
---|---
Figure 8: As in Figure 7, $\hat{y}_{*}(s)$ (left) and $a_{*}(s)$ (right). The
values $s$ considered as the ones such that
$\hat{y}_{\min}\leq\hat{y}_{*}(s)\leq\frac{1}{2}$ and $a_{\min}\leq a\leq
a_{\max}$ in km.

Summarizing, for the values of $s_{\min},s_{\max}$ in the lemma we have that
for any $s\in[s_{\min},s_{\max}]$,

$\partial_{\hat{y}}\hat{\Gamma}_{-}\big{(}\hat{y}_{*}(s);L_{*}(s)\big{)}=0$

with $L_{*}(s)=\sqrt{\mu a_{*}(s)}$ and $\hat{y}_{*}(s),a_{*}(s)$ defined in
(5.18). In addition, these are the only possible critical points of
$\hat{\Gamma}_{-}$. ∎

We introduce now the boundary values of $\hat{\Gamma}$. We first define, for
$L\in[L_{*}(s_{\max}),L_{\max}]$,

$\hat{\Gamma}_{*}^{-}(L)=\hat{y}_{*}(s)s,\qquad\text{with}\quad
s=s(L)\quad\text{ such that}\quad L=L_{*}(s).$

Then we introduce, for $L\in[L_{\min},L_{\max}]$

$\hat{\Gamma}_{\min}^{-}(L)=\hat{\Gamma}_{-}(\hat{y}_{\min};L),\qquad\hat{\Gamma}_{\max}^{-}(L)=\hat{\Gamma}_{-}\left(\frac{1}{2};L\right)$
(5.19)

that, by Corollary 5.5, satisfy

$\displaystyle\hat{\Gamma}_{\min}^{-}(L)$
$\displaystyle\in\left[\hat{y}_{\min}\frac{2}{5}m,\hat{y}_{\min}m\right]=[0.005433756726029,0.013584391815074]$
(5.20) $\displaystyle\hat{\Gamma}_{\max}^{-}(L)$
$\displaystyle\in\left[\frac{1}{5}m,\frac{1}{2}m\right]=[0.023303027798234,0.058257569495584].$

By Lemma 5.9, $\hat{y}_{*}(s)$ is an increasing function whereas $L_{*}(s)$ is
decreasing. Let $L_{*}^{-1}$ be its inverse. We have that
$\hat{\Gamma}_{*}^{-}(L)=\hat{y}_{*}(L^{-1}_{*}(L))L^{-1}_{*}(L)$ is a
decreasing function. Therefore

$\hat{\Gamma}_{*}^{-}(L)\in[\hat{\Gamma}_{*}^{-}(L_{\max}),\hat{\Gamma}_{*}^{-}(L_{*}(s_{\max}))]$

and, by definition of $s_{\min},s_{\max}$ in Lemma 5.9,

$\hat{\Gamma}_{*}^{-}(L)\in\left[\hat{y}_{*}(s_{\min})s_{\min},\frac{1}{2}s_{\max}\right]=[0.038367525991514,0.048398408167370].$

Next lemma studies the monotonicity properties of the function
$\hat{\Gamma}_{-}$.

######  Lemma 5.10.

For $\hat{y}_{\min}\leq\hat{y}\leq\frac{1}{2}$, the function
$\hat{\Gamma}_{-}(\hat{y};\cdot)$ is strictly decreasing (with respect to $L$)
on $[L_{\min},L_{\max}]$. Therefore, $\hat{\Gamma}_{\min,\max}^{-}(L)$ (see
definition (5.19)) are also decreasing and moreover

$\hat{\Gamma}_{\min}^{-}(L)\in[\hat{\Gamma}_{\min}^{-}(L_{\max}),\hat{\Gamma}_{\min}^{-}(L_{\min})]=[0.013575332545548,0.013584387878826]$

and

$\hat{\Gamma}_{\max}^{-}(L)\in[\hat{\Gamma}_{\max}^{-}(L_{\max}),\hat{\Gamma}_{\max}^{-}(L_{\min})]=[0.027639200647529,0.058244177051364]$

where the values have been computed numerically.

The function $\hat{\Gamma}_{-}(\cdot;L)$ satisfies, for
$\hat{y}_{\min}\leq\hat{y}\leq\frac{1}{2}$, that

  1. 1.

For $L\in[L_{\min},L_{*}(s_{\max})]$, the function $\hat{\Gamma}_{-}(\cdot;L)$
is strictly increasing and moreover
$\hat{\Gamma}_{-}(\hat{y};L)\in[\hat{\Gamma}_{\min}^{-}(L),\hat{\Gamma}_{\max}^{-}(L)]$.

  2. 2.

For $L\in(L_{*}(s_{\max}),L_{\max}]$, $\hat{\Gamma}_{-}(\cdot,L)$ has a
maximum at $\hat{y}_{\max}(L)=\hat{y}_{*}(L_{*}^{-1}(L))$. In this case,

$\hat{\Gamma}_{-}(\hat{y};L)\in[\hat{\Gamma}_{\min}^{-}(L),\hat{\Gamma}_{*}^{-}(L)].$

######  Remark 5.11.

The statement of Lemma 5.10 can be graphically represented as in Figure 9.

Figure 9: For a given value of $L$, the figures represent qualitatively the
function $\hat{\Gamma}_{-}(\hat{y},L)$ (as a function of $\hat{y}$). There are
depicted the three different behaviours depending on the value of $L$. See
also Figure 10 for numerical computations.

###### Proof.

From formula (5.13) of $\partial_{L}\hat{\Gamma}_{-}$ and using that, by Lemma
5.1, $\partial_{s}\mathcal{G}>0$, we have that the sign of
$\partial_{L}\hat{\Gamma}_{-}$, is the same as for
$-\partial_{L}\mathcal{G}_{-}(s_{-},\hat{y},L)$ or, in other words, the same
as $-\mathbf{b}_{-}(s_{-})-\hat{y}^{2}\mathbf{c}_{-}(s_{-})$ with
$s_{-}=s_{-}(\hat{y},L)$. Again, from Lemma 5.1, we have that
$\mathbf{b}_{-},\mathbf{c}_{-}$ are increasing functions and using their
values (5.10) and (5.11) at $s=\frac{2}{5}m$, we conclude that
$\mathbf{b}_{-}(s_{-})+\hat{y}^{2}\mathbf{c}_{-}(s_{-})>0$ and so
$\hat{\Gamma}_{-}(\hat{y};\cdot)$ is decreasing.

By Lemma 5.9, when $L\in[L_{\min},L_{*}(s_{\max})]$ there are no critical
points of $\hat{\Gamma}_{-}$ in the interval
$\left[\hat{y}_{\min},\frac{1}{2}\right)$. That implies that
$\hat{\Gamma}_{-}(\cdot;L)$ is either increasing or decreasing. Notice that,
by definition (5.19) of $\hat{\Gamma}_{\min,\max}^{-}$ and bounds (5.20),

$\hat{\Gamma}_{-}(\hat{y}_{\min};L)=\hat{\Gamma}_{\min}(L)\leq\hat{\Gamma}_{\max}(L)=\hat{\Gamma}_{-}\left(\frac{1}{2};L\right)$

and therefore $\hat{\Gamma}_{-}(\cdot;L)$ is a strictly increasing function in
$[L_{\min},L_{*}(s_{\max}))$.

When $L=L_{*}(s_{\max})$ we have that the corresponding critical point is
$\hat{y}=y_{*}(s_{\max})=\frac{1}{2}$. Therefore, the same argument as before
allows to conclude that also in this case $\hat{\Gamma}_{-}(\cdot;L)$ is
strictly increasing for $\hat{y}_{\min}\leq\hat{y}<\frac{1}{2}$.

Take now $L\in(L_{*}(s_{\max}),L_{\max}]$ and let
$s=L_{*}^{-1}(L)\in[s_{\min},s_{\max}]$ be such that $L=L_{*}(s)$. By Lemma
5.9, we already know that $\hat{y}_{*}(s)$ is a critical point of
$\hat{\Gamma}_{-}(\cdot;L)$ and that
$\hat{\Gamma}_{-}(\hat{y}_{*}(s);L_{*}(s))=\hat{\Gamma}_{*}^{-}(L)$. So we
only need to check that

$\partial^{2}_{\hat{y}}\hat{\Gamma}_{-}(\hat{y}_{*}(L^{-1}_{*}(L));L)=\partial^{2}_{\hat{y}}\hat{\Gamma}_{-}(\hat{y}_{*}(s);L_{*}(s))<0,\qquad
s=L^{-1}_{*}(L).$

We first recall expression (5.15) of $\hat{\Gamma}_{-}$ in (5.15):

$\partial_{\hat{y}}\hat{\Gamma}_{-}(\hat{y};L)=\frac{\mathcal{H}(s_{-}(\hat{y};L),\hat{y};L)}{\partial_{s}\mathcal{G}_{-}(s_{-}(\hat{y};L),\hat{y};L)}.$

Then, using that, by construction, $\mathcal{H}(s,\hat{y}_{*}(s);L_{*}(s))=0$
and $s_{-}(\hat{y}_{*}(s);L_{*}(s))=s$

$\partial_{\hat{y}}^{2}\hat{\Gamma}_{-}(\hat{y}_{*}(s);L_{*}(s))=\frac{\partial_{\hat{y}}\mathcal{H}(s,\hat{y}_{*}(s);L_{*}(s)),}{\partial_{s}\mathcal{G}_{-}(s,y_{*}(s);L_{*}(s))},$

and, since by Lemma 5.1, $\partial_{s}\mathcal{G}>0$, we need to compute the
sign of $\partial_{\hat{y}}\mathcal{H}(s,\hat{y}_{*}(s);L_{*}(s))$. From
(5.16) we have that

$\partial_{\hat{y}}\mathcal{H}(s,\hat{y};L)=2\alpha^{3}L^{4}\rho\big{[}3\overline{\mathbf{b}}(s)\hat{y}^{2}+5\overline{\mathbf{c}}(s)\hat{y}^{4}].$

Then, using that $\mathcal{H}(s,\hat{y}_{*}(s);L_{*}(s))=0$, we deduce that

$\partial_{\hat{y}}\mathcal{H}(s,\hat{y}_{*}(s);L_{*}(s))=\frac{1}{\hat{y}_{*}(s)}\left(-3\overline{\mathbf{a}}(s)+4\alpha^{3}L^{4}\rho\overline{\mathbf{c}}(s)(\hat{y}_{*}(s))^{5}\right).$

On the other hand, from expression (5.17) for $\overline{\mathbf{c}}$,

$\overline{\mathbf{c}}(s)\leq 8U_{2}^{0,0}(5-12s)$

where we have used that $45-60s+5s^{2}+8s^{3}\geq 45-60s\geq 45-60m>0$.
Therefore, using again that $s\in(\frac{2}{5}m,m)$,
$0<\hat{y}_{*}(s)\leq\frac{1}{2}$ and $a_{\min}\leq a\leq a_{\max}$, we obtain

$\displaystyle\hat{y}_{*}(s)\partial_{\hat{y}}\mathcal{H}(s,\hat{y};L)$
$\displaystyle\leq-6s(5s+4)+32\alpha^{3}L^{4}\rho(\hat{y}_{*}(s))^{5}U_{2}^{0,0}(5-12s)$
$\displaystyle\leq-6\cdot\frac{2}{5}m(2m+4)+\frac{a_{\max}^{3}}{a_{\mathrm{M}}^{3}}L_{\max}^{4}\rho
U_{2}^{0,0}\left(5-\frac{24}{5}m\right)$ $\displaystyle=-0.778057059724233<0.$

As a conclusion $\partial_{\hat{y}}\mathcal{H}(s_{-},\hat{y};L)<0$ if
$(\hat{y},L)=(\hat{y}_{*}(s),L_{*}(s))$ and that implies that
$\partial_{\hat{y}}^{2}\hat{\Gamma}_{-}(\hat{y}_{*}(s);L_{*}(s))<0$ so that
$\hat{y}_{*}(s)$ is a maximum of $\hat{\Gamma}(\cdot;L_{*}(s))$. ∎

######  Remark 5.12.

The values of $\hat{\Gamma}_{-}(\hat{y};L)$ can be numerically computed for
any fixed $\hat{y},L$ as $\hat{\Gamma}_{-}(\hat{y};L)=s_{-}(\hat{y};L)\hat{y}$
with $s_{-}$ the zero of the function $\mathcal{G}_{-}(s,\hat{y};L)$. In
Figure 10 we present some representative values of $a$ (recall that
$L=\sqrt{\mu a})$ where we can find the different behaviour described in Lemma
5.10.

|
---|---
|
Figure 10: The function $\hat{\Gamma}_{-}(\hat{y},L)$ for
$a=22000,a_{*}(s_{\max}),26000,29600,\,\mathrm{km}$. Notice that for
$a=22000,a_{*}(s_{\max})$ the function is strictly increasing but when
$a=a_{*}(s_{\max})$ $\hat{\Gamma}_{-}(\hat{y};L_{*}(s_{\max}))$ has a critical
point at $\hat{y}=\frac{1}{2}$. The black dotted line is for
$\hat{\Gamma}_{\max}^{-}(L)$ meanwhile the orange one is for
$\hat{\Gamma}_{*}^{-}(L)$.

######  Remark 5.13.

We first notice that
$\hat{\Gamma}_{*}^{-}(L_{*}(s_{\max}))=\hat{\Gamma}_{\max}^{-}(L_{*}(s_{\max})).$
Then, using that by Lemma 5.10, $\hat{\Gamma}_{\max}^{-}$ is decreasing (in
its variable $L$), we also have that for $L\in(L_{*}(s_{\max}),L_{\max}]$,

$\hat{\Gamma}^{-}_{\max}(L)<\hat{\Gamma}^{-}_{*}(L).$

Moreover, since $s_{\max}=s_{-}\left(\frac{1}{2},L_{*}(s_{\max})\right)$, we
have that $s_{\max}\leq m$. Therefore, using that by (5.14),
$\hat{\Gamma}^{+}_{\max}(L)\geq\frac{1}{2}m$, we have that

$\hat{\Gamma}^{-}_{*}(L)\leq\frac{1}{2}s_{\max}\leq\frac{1}{2}m\leq\hat{\Gamma}^{+}_{\max}(L).$

From (5.20) it is also clear that the constant
$\hat{\Gamma}_{\min}^{+}=\hat{y}_{\min}m$ satisfies that

$\hat{\Gamma}_{\min}^{-}(L)\leq\hat{\Gamma}_{\min}^{+}\leq\hat{\Gamma}_{\max}^{-}(L).$

From this analysis, it is straightforward to deduce the following result about
the existence of equilibrium points of the form $(\pi,y)$ (that is, the
existence results of Theorem 3.5).

######  Proposition 5.14.

Let $L\in[L_{\min},L_{\max}]$. If
$\hat{\Gamma}\in[\hat{\Gamma}_{\min}^{-}(L),\hat{\Gamma}_{\max}^{-}(L)]$,
there exists a unique $\hat{y}_{-}(\hat{\Gamma};L)$ such that system (4.1) has
a fixed point of the form $(\pi,y)=(\pi,L\hat{y}_{-}(\hat{\Gamma};L))$. The
function $\hat{y}_{-}(\cdot;L)$ is strictly increasing.

In addition, for $L\in[L_{\min},L_{*}(s_{\max})]$:

  1. 1.

If
$\hat{\Gamma}\notin[\hat{\Gamma}_{\min}^{-}(L),\hat{\Gamma}_{\max}^{-}(L)]$,
system (4.1) has no fixed points of the form $(\pi,L\hat{y})$ with
$\hat{y}_{\min}\leq\hat{y}\leq\frac{1}{2}$ and
$\hat{\Gamma}\in\left(0,\frac{1}{2}\hat{y}\right)$.

and when $L\in(L_{*}(s_{\max}),L_{\max}]$, we have that

  1. 1.

If $\hat{\Gamma}\in[\hat{\Gamma}_{\max}^{-}(L),\hat{\Gamma}_{*}(L))$, there
exist only two functions $\hat{y}_{-}^{1,2}(\hat{\Gamma};L)$ satisfying

$\hat{y}_{-}^{1}(\hat{\Gamma};L)\leq\hat{y}_{\max}(L)\leq\hat{y}_{-}^{2}(\hat{\Gamma};L)$

with $\hat{y}_{\max}(L)$ defined in Lemma 5.10, such that system (4.1) has two
fixed point of the form $(\pi,y)=(\pi,L\hat{y}_{-}^{1,2}(\hat{\Gamma};L))$. In
addition, $\hat{y}^{1}_{-}(\cdot;L)$ is strictly increasing and
$\hat{y}^{2}_{-}(\cdot;L)$ is strictly decreasing.

  2. 2.

If $\hat{\Gamma}=\hat{\Gamma}_{*}(L)$, there exists only one fixed point
$(\pi,\hat{y}_{\max}(L))$ with $\hat{y}_{\max}(L)$ a decreasing function.

  3. 3.

If $\hat{\Gamma}\notin[\hat{\Gamma}_{\min}^{-}(L),\hat{\Gamma}_{*}(L)]$,
system (4.1) has no fixed points of the form $(\pi,L\hat{y})$ with
$\hat{y}_{\min}\leq\hat{y}\leq\frac{1}{2}$ and
$\hat{\Gamma}\in\left(0,\frac{1}{2}\hat{y}\right)$.

###  5.6 Linearization around the critical points

To complete the proof of Theorem 3.5, we are going to study the character of
the fixed points of the form $(0,y),(\pi,y)$. Using the decomposition of
${\mathcal{B}}_{0,\alpha}$ in (4.4), we write system (4.1) as

$\displaystyle\dot{x}$
$\displaystyle=-g(y;L)\big{[}X_{1}(y,\Gamma;L)+X_{2}(y,\Gamma;L)\cos x]$
$\displaystyle\dot{y}$ $\displaystyle=h(y,\Gamma;L)\sin x$

with

$g(y;L)=\frac{1}{256}\frac{\rho_{0}}{L^{4}y^{6}},\qquad
h(y,\Gamma;L)=\frac{1}{128}\frac{\rho\alpha^{3}}{L^{4}y^{5}}\sqrt{P_{1}(y,\Gamma)}\mathcal{A}_{1}(y,\Gamma;L)$

and

$X_{1}(y,\Gamma;L)=B_{0}(y,\Gamma;L)+\alpha^{3}B_{1}(y,\Gamma;L),\qquad
X_{2}(y,\Gamma;L)=\alpha^{3}\frac{{\mathcal{B}}_{1}(y,\Gamma;L)}{\sqrt{P_{1}(y,\Gamma)}}.$

By Corollary 5.5, we can characterize the fixed points as the sets

$\\{(0,y,\Gamma_{+}(y;L),L)\\},\qquad\\{(\pi,y,\Gamma_{-}(y;L);L)\\}$

with $L\in[L_{\min},L_{\max}]$ and $\hat{y}_{\min}L\leq y<\frac{L}{2}$. They
satisfy $\mathcal{F}_{\pm}(y,\Gamma_{\pm}(y;L);L)=0,$ where
$\mathcal{F}_{\pm}$ are defined in (5.3), the sign $+$ corresponds to $(0,y)$
and $-$ otherwise. Equivalently we have that

$\mathcal{F}_{\pm}(y,\Gamma_{\pm}(y;L);L)=X_{1}(y,\Gamma_{\pm}(y;L);L)\pm
X_{2}(y,\Gamma_{\pm}(y;L);L)=0.$

The variational equation at the fixed points $(0,y,\Gamma_{+}(y;L),L)$ or
$(\pi,y,\Gamma_{-}(y;L),L)$, $\dot{z}=M^{L}(y)z$, is given by

$M^{L}(y):=\left(\begin{array}[]{cc}0&-g(y;L)\big{[}\partial_{y}X_{1}(y,\Gamma_{\pm};L)\pm\partial_{y}X_{2}(y,\Gamma_{\pm};L)\big{]}\\\
\pm h(y,\Gamma_{\pm};L)&0\end{array}\right),$

where $\Gamma_{\pm}=\Gamma_{\pm}(y;L)$. The eigenvalues of $M^{L}$ then
satisfy