# On the Arnold diffusion mechanism in Medium Earth Orbit

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

###### Abstract

Space debris mitigation guidelines represent the most effective method to
preserve the circumterrestrial environment. Among them, end-of-life disposal
solutions play a key role. In this regard, effective strategies should be
conceived not only on the basis of novel technologies, but also following an
advanced theoretical understanding. A growing effort is devoted to exploit
natural perturbations to lead the satellites towards an atmospheric reentry,
reducing the disposal cost, also if departing from high-altitude regions. In
the case of the Medium Earth Orbit region, home of the navigation satellites
(like GPS and Galileo), the main driver is the gravitational perturbation due
to the Moon, that can increase the eccentricity in the long term. In this way,
the pericenter altitude can get into the atmospheric drag domain and the
satellite can eventually reenter.

In this work, we show how an Arnold diffusion mechanism can trigger the
eccentricity growth. Focusing on the case of Galileo, we consider a hierarchy
of Hamiltonian models, assuming that the main perturbations on the motion of
the spacecraft are the oblateness of the Earth and the gravitational
attraction of the Moon. First, the Moon is assumed to lay on the ecliptic
plane and periodic orbits and associated stable and unstable invariant
manifolds are computed for various energy levels, in the neighborhood of a
given resonance. Along each invariant manifold, the eccentricity increases
naturally, achieving its maximum at the first intersection between them. This
growth is, however, not sufficient to achieve reentry. By moving to a more
realistic model, where the inclination of the Moon is taken into account, the
problem becomes non-autonomous and the satellite is able to move along
different energy levels. Under the ansatz of transversality of the stable and
unstable manifolds in the autonomous case, checked numerically, Poincaré-
Melnikov techniques are applied to show how the Arnold diffusion can be
attained, by constructing a sequence of homoclinic orbits that connect
invariant tori at different energy levels on the normally hyperbolic invariant
manifold.

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
    1. 1.1 An Arnold diffusion mechanism
  2. 2 The secular Hamiltonian system and the main result
  3. 3 A good system of coordinates
    1. 3.1 Slow-fast coordinates $(y,x)$
    2. 3.2 Poincaré coordinates $(\eta,\xi)$
  4. 4 The hierarchy of models and dynamical systems tools
    1. 4.1 The hierarchy of models
      1. 4.1.1 Theoretical results for the coplanar and $h$-averaged system
    2. 4.2 Normally hyperbolic invariant manifolds and scattering maps
    3. 4.3 The Arnold diffusion instability mechanism
  5. 5 The coplanar secular Hamiltonian
    1. 5.1 Periodic orbits and their invariant manifolds
      1. 5.1.1 Numerical verification of Ansatz 5.1
    2. 5.2 The invariant cylinder and their invariant manifolds
      1. 5.2.1 The coplanar inner map
      2. 5.2.2 The coplanar scattering maps
  6. 6 Dynamics of the system with $i_{\mathrm{M}}>0$
    1. 6.1 Perturbative analysis of the inner and scattering maps
    2. 6.2 Existence of diffusing orbits
    3. 6.3 Proof of Theorem 6.7
  7. 7 Conclusions
  8. A Expression of the Hamiltonian
    1. A.1 Hamiltonian in slow-fast coordinates $(y,x)$
    2. A.2 Hamiltonian in Poincaré coordinates $(\eta,\xi)$

##  1 Introduction

One of the key actions to mitigate the space debris problem and ensure a
sustainable exploitation of the circumterrestrial environment is to implement
effectively end-of-life procedures. While for the Low Earth Orbit and the
Geosynchronous protected regions, there exist well-defined guidelines [35],
for the Medium Earth Orbit (MEO) region still a discussion is ongoing. This is
due to the fact that the MEO region is very broad (in principle, it goes from
an altitude of 2000 km up to the geostationary altitude) and not yet critical
in terms of population density.

The MEO region is mostly known because it is where the satellites of the
Global Navigation Satellite Systems (GNSS) orbit, namely, GPS, Galileo,
GLONASS and Beidou-M [3]. They cover a range of semi-major axis between about
$25500$ km (GLONASS) and $29600$ km (Galileo) and their nominal inclination is
$55^{\degree}\pm 2^{\degree}$ (GPS and Beidou-M), $56^{\degree}\pm
2^{\degree}$ (Galileo), $65^{\degree}\pm 2^{\degree}$ (GLONASS) (see, e.g.,
[50, 2] and references herein). In the recent ESA’s zero debris policy [20],
the GNSS orbits are defined as “valuable orbits” and the possibility of
extending the protection to this region is introduced. Given the high
altitude, so far at the end of life the GNSS satellites were either left in
the operational orbit or re-orbited by a given amount [3]. The accumulation of
non-operational satellites in a limited region eventually will lead to a
critical situation in terms of potential collisions and thus fragmentation.
For this reason, Jenkin and Gick [36] proposed to dilute the collision
probability by increasing the orbital eccentricity of the satellites at the
end of life. Indeed, the GNSS inclination and altitude are such that in the
long term the third-body perturbation, that is, the gravitational perturbation
exerted by Sun and Moon, could lead to a natural eccentricity growth up to
reentry, if a suitable initial orbital orientation is chosen. This mechanism
was proven numerically by several authors (see, e.g., [50, 48, 2, 4, 32, 44,
45]), but an exhaustive theoretical explanation of the underlying mechanism is
still missing. In particular, the singular resonance hypothesis, that is, an
integrable model, cannot explain the maximum value of eccentricity that can be
observed by considering a full dynamical model, nor the Saros periodicity
detected.

More precisely, following [39] the disturbing function corresponding to the
gravitational perturbation due to a third body can be written as a series
expansion depending on semi-major axis, eccentricity, inclination of both the
satellite and the third body, and a periodic term involving the longitude of
the ascending node, the argument of perigee and the mean anomaly of both the
satellite and the third body. By doubly-averaging the periodic terms (over the
orbital period of the satellite and the orbital period of the Moon), and
considering only the first-order effect, we define the secular
Hamiltonian.111In what follows, this step will be omitted because it can be
found in several past works (see references above).

A possible way to deal with the secular Hamiltonian is to assume that only one
periodic term is dominant at a time, in particular, when its argument is
resonant. Most of the past works222It is beyond the scope of this work to
provide a review of all the past investigations on the subject. Here, we
recall only the works that are relevant to our contribution. focus on such
“isolated resonance hypothesis” and on the eccentricity growth associated with
an inclination-dependent resonance, involving the argument of the pericenter
$\omega$ and the longitude of the ascending node of the satellite $\Omega$.
For quasi-circular orbits at MEO altitudes, and assuming the Galileo
inclination ($\approx 56^{\degree}$), the dominant resonance is
$2\dot{\omega}+\dot{\Omega}\approx 0$333For GLONASS the inclination is
$\approx 63^{\degree}$ and thus the dominant resonance is $2\dot{\omega}$.. As
a first approximation, $\omega$ and $\Omega$ change in time because of the
Earth’s oblateness, but the third body perturbation has also an effect.

The authors in [49, 12] assumed the isolated resonance hypothesis and brought
forward the concept that the eccentricity growth that ensures a reentry is due
to a chaotic behavior, that occurs when two or more resonances overlap. The
Chirikov criterion is mainly proven by detailed Fast Lyapunov Indicators (FLI)
maps that show the location of the chaotic regions. The phase space associated
with each possible resonance of the secular dynamics was investigated also in
[42], but in this case, the authors remarked the role of the Laplace plane.
More recently, the authors in [31, 11, 41] proposed the idea that the high
eccentricity growth is due, not to the chaos generated by overlapping
resonances, rather to the normally hyperbolic invariant manifold (NHIM)
associated with the given resonance and through an Arnold diffusion mechanism.
They showed, with refined FLI maps and an ad hoc Hamiltonian derived for
circular orbits in the neighborhood of the resonance, that what triggers the
phenomenon is the variation of a given integral of motion. They also
emphasized the role of the Laplace plane for the motion on the invariant tori.
Some of the invariant objects “necessary” to achieve Arnold diffusion were
already identified in [11, 41]. These papers also showed, numerically,
hyperbolicity and unstable motions along the resonance.

The goal of the present paper is to explain the mechanism that creates the
drift and how to construct drifting orbits across the neighborhood of an
isolated resonance. In particular, we focus our study on the resonance
$2\omega+\Omega$ with the Galileo satellites values. Nonetheless, our approach
is fairly general and can be adapted for other resonances and values. In other
words, we follow the idea that it is the NHIM the cause that we have to
investigate, and we will show how it is possible to construct homoclinic
connections that “connect” different energy levels on the NHIM associated with
the $2\omega+\Omega$ resonance. Moreover, we will show how to construct orbits
that shadow (follow closely) this sequence of homoclinic orbits. These
shadowing orbits achieve a drift in eccentricty large enough so that the
satellite renters the Earth’s atmosphere.

We will focus only on the third-body effect exerted by the Moon, thus
neglecting the Sun, and we will apply a perturbative approach with respect to
the inclination of the Moon on the ecliptic plane $i_{\mathrm{M}}$ (which has
a real value of $i_{\mathrm{M}}=5.15^{\circ}$). More details on this will be
given in the next subsection. Notice that the model we consider is a rather
simplified version of the initial problem. However, we expect that the outcome
of this work will pave the way for the analysis of a more realistic model,
that can be used to define effective end-of-life solutions.

The paper is structured as follows. In Section 2 we define the secular
Hamiltonian and provide formulas for it. Then, we state the main result of
this paper. In Section 3 we introduce a “good” system of coordinates which
captures the timescales of the model. In Section 4, relying on what we call
hierarchy of models, we explain how the Arnold diffusion mechanism takes place
in the secular Hamiltonian. Finally, we also provide definitions for the main
tools we use: a normally hyperbolic invariant manifold and the associated
scattering maps. In Section 5 we analyze the dynamics for the coplanar model
(that is, taking $i_{\mathrm{M}}=0)$. We use this analysis to describe, in
Section 6, the dynamics for the full secular Hamiltonian with
$i_{\mathrm{M}}>0$ small enough.

###  1.1 An Arnold diffusion mechanism

V. Arnold in 1964 (see [5]) showed, in a cleverly chosen model, that actions
can vary drastically in nearly integrable Hamiltonian systems (of at least 2
and a half degrees of freedom, that is a phase space of at least dimension 5).
Then, he conjectured that such behavior, nowadays called Arnold diffusion,
should be typical. In (nearly integrable) physical models, it is expected that
Arnold diffusion is a fundamental mechanism leading to transport in phase
space. Such transport is achieved by drifting along resonances.

The understanding of Arnold diffusion mechanisms has had outstanding progress
in the last decades, relying on a wide variety of techniques: the original
geometric approach by Arnold, which has been deeply developed in [7, 13, 14,
15, 18, 26, 19, 27], variational methods [9, 8], topological tools [29, 10],
the so-called separatrix map [53, 54] or a combination of different approaches
[6, 37].

As stated before, the goal of this paper is to explain how an Arnold diffusion
mechanism can enhance drift in eccentricity along the $2\omega+\Omega$
resonance. The mechanism we propose relies on geometric tools in the spirit of
the seminal work by Arnold. They are explained in detail in Section 4.2.

Let us mention here just the main ingredients. We show that, along the
resonance, there exists a normally hyperbolic invariant cylinder (see
Definition 4.3 below). The stable and unstable manifolds of this cylinder
intersect transversally. Thus, we can construct trajectories that, following
closely a sequence of homoclinic orbits to the cylinder, achieve a
considerable drift in eccentricity. Such drift allows the satellite, by
flattening its osculating ellipse, to enter the Earth’s atmosphere.

To perform such analysis, we assume that the inclination of the Moon
$i_{\mathrm{M}}$ with respect to the ecliptic plane is small enough. This
allows us to use perturbative arguments to construct the “highways” that lead
to the drift. Note that it is fundamental that $i_{\mathrm{M}}>0$, because in
the coplanar case, namely, when $i_{\mathrm{M}}=0$, there are not enough
dimensions to achieve Arnold diffusion.

Even if our tools rely on the fact that $i_{\mathrm{M}}>0$ is small enough, we
expect that the same mechanism takes place for a realistic value of
$i_{\mathrm{M}}$, because, indeed, it is expected that transport phenomena are
even stronger and more robust in far from integrable Hamiltonians.
Nevertheless, in that setting, one cannot combine analytical and numerical
techniques to describe the cylinder and its invariant manifolds, as we do in
the present paper. Instead, one has to describe them fully numerically. In
this work, numerical simulations are done for a lower dimensional model (what
we call the coplanar Hamiltonian) and can be kept “simple”. Otherwise, it
would require to compute numerically high dimensional objects in the
5-dimensional phase space. Moreover, when $i_{\mathrm{M}}$ increases other
resonances could come to play also an important role and an resonance
overlapping could take place.

Finally, note that, for applications, it is fundamental to know the speed of
such transport mechanisms. For arbitrarily small $i_{\mathrm{M}}>0$ the
mechanism is rather slow (the drifting time is $T\sim 1/i_{\mathrm{M}}$).
However, for a realistic value of $i_{\mathrm{M}}$ it is expected to be
faster. Moreover, one may expect that combining Arnold diffusion mechanism
with maneuvers could speed up considerably the drifting.

## Acknowledgements

E.M. Alessi acknowledges 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 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 The secular Hamiltonian system and the main result

Let us consider a spacecraft that is affected by the gravitational attraction
of the Earth, the perturbation due to the Earth’s oblateness and the lunar
gravitational perturbation. In the geocentric equatorial reference system, the
motion of the spacecraft takes place on an ellipse, that is described by the
orbital elements semi-major axis $a$, eccentricity $e$, inclination $i$,
longitude of the ascending node $\Omega$ and argument of pericenter $\omega$.
The ellipse changes in time due to 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}$, $e_{\mathrm{M}}=0.0549006$,
$i_{\mathrm{M}}=5.15^{\degree}$, while the longitude of the ascending node of
the Moon with respect to the ecliptic plane varies approximately linearly with
time in a period $\mathcal{T}_{\Omega_{\mathrm{M}}}$ of 1 Saros (about
6585.321347 days [51, 46]) 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/\mathcal{T}_{\Omega_{\mathrm{M}}},$ (2.1)

where $\Omega_{{\mathrm{M}},0}$ is the longitude of ascending node of the Moon
at a given epoch.

Since we consider the secular model (the system averaged over the mean
anomalies of both the satellite and the moon), the semi-major axis $a$ is a
constant of motion which, in the case of Galileo, corresponds to $a=29600$ km.
Let us define

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

which characterizes the Earth-satellite distance with respect to the Earth-
Moon distance.

Then, in Delaunay action-angle variables444Recall that in celestial mechanics,
these are the classical canonical variables, defined as  $\displaystyle L$
$\displaystyle=\sqrt{\mu a},$ $\displaystyle\quad l$ $\displaystyle=M,$ (2.2)
$\displaystyle G$ $\displaystyle=L\sqrt{1-e^{2}},$ $\displaystyle\quad g$
$\displaystyle=\omega,$ $\displaystyle H$ $\displaystyle=G\cos i,$
$\displaystyle\quad h$ $\displaystyle=\Omega.$ , the secular dynamics is
described by the Hamiltonian [12]

$\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.3)

where

${\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}$ the mass parameter of the Earth,
and

$\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. The secular perturbative term due to the Moon is instead

$\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}$

The function $U_{2}^{m,\mp s}(\epsilon)$ corresponds to the Giacaglia function
with $\epsilon=23.44^{\degree}$ (see Table 2 in Appendix A.1) being the
obliquity of the ecliptic with respect to the equatorial plane. Also, one has
that

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

where $\tilde{F}_{m,p}(G,H)$ is a function of the Kaula’s inclination
functions $F_{2,m,p}(i)$555Since $H=G\cos i$, one has that
$\tilde{F}_{m,p}(G,H)=F_{2,m,p}(\arccos\frac{H}{G}).$, see [38], and
$\tilde{X}_{p}(L,G)$ is a function of the zero-order Hansen coefficient
$X_{0}^{2,2-2p}(e)$666Since $G=L\sqrt{1-e^{2}}$, one has that
$\tilde{X}_{p}(L,G)=X_{0}^{2,2-2p}\left(\sqrt{1-\frac{G^{2}}{L^{2}}}\right)$.
, see [34], and

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

where $\sigma=\pm s$ defines the relative orientation satellite-Moon. In
addition, the other variables are constants 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.7) $\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 [38]), one obtains that

$\displaystyle
F_{2,s,1}(i_{\mathrm{M}})=\left\\{\begin{array}[]{ll}-\frac{1}{2}+\frac{3}{4}\sin^{2}i_{\mathrm{M}}&\mbox{if
$s=0$,}\\\\[5.0pt] -\frac{3}{2}\sin i_{\mathrm{M}}\cos i_{\mathrm{M}}&\mbox{if
$s=1$,}\\\\[5.0pt] \frac{3}{2}\sin^{2}i_{\mathrm{M}}&\mbox{if
$s=2$.}\end{array}\right.$ (2.11)

We remark that when $i_{\mathrm{M}}=0$, the inclination function
$F_{2,s,1}(0)=0$ for all $s$, except for $s=0$ which is $F_{2,0,1}(0)=-1/2$.
In other words, assuming that the Moon lies on the ecliptic plane, the
Hamiltonian $\mathtt{H}$ in (2.3) 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$.

######  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 $a$ of the orbit of the Gallileo satellites (equal to
29600 km) is the unit of distance and the corresponding orbital period is
2$\pi$.

The main result of the paper is the following. It assumes several ansätze that
are stated below and are verified numerically.

######  Theorem 2.3.

Consider the secular Hamiltonian $\mathtt{H}$ in (2.3) with the parameters
just fixed, take

$\alpha=0.077$ (2.12)

and assume that the Ansätze 5.1, 5.6 and 6.6 are satisfied.

Then, for $i_{\mathrm{M}}>0$ small enough, there exist a time $T>0$ and a
trajectory $z(t)=(L(t),G(t),H(t),g(t),h(t),\Omega_{\mathrm{M}}(t))$ such that

$\mathtt{H}(z(0))=1.7\cdot
10^{-8}\qquad\text{and}\qquad\mathtt{H}(z(T))=1.3\cdot 10^{-6}.$

Moreover, along the energy increase the osculating eccentricity and
inclination satisfy

$e(0)=0\qquad\text{and}\qquad e(T)=0.78,$

and

$i(t)\in[56.06^{\degree},58.09^{\degree}]\qquad\text{for}\qquad t\in[0,T].$

Note that such evolution would not be possible if $i_{\mathrm{M}}=0$ since in
that case the Hamiltonian is autonomous and therefore the energy is a first
integral (recall Remark 2.1).

Concerning the evolution of the eccentricity, we note that it is already
oscillating when $i_{\mathrm{M}}=0$. However, the size of this oscillations
depends on the energy value. Indeed, the orbit that we obtain in Theorem 2.3
is such that at the initial time, assuming $i=56.06^{\degree}$ the
eccentricity is oscillating (approximately) between 0 and 0.35, whereas at the
final time $T$ is oscillating between 0 and 0.78. Therefore, at the initial
time the satellite is far from the Earth’s atmosphere whereas at the final
time the satellite is reentering.

Theorem 2.3 does not provide any estimate on the time $T>0$ needed to achieve
the drift in eccentricity. As mentioned before, more quantitative shadowing
arguments compared to those used in the present paper (see, for instance [52,
47]) should lead to estimates for $T=T(i_{\mathrm{M}})$ of the form

$|T|\leq\frac{C}{i_{\mathrm{M}}},$

for some constant $C>0$ independent of $i_{\mathrm{M}}>0$ (see Remark 4.5
below).

Finally, 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. As well, we omit the dependence of the Hamiltonian on
the variable $L$ since it is a constant of motion.

##  3 A good system of coordinates

Let us start by focusing on the region of the phase space where an orbit
satisfies the “$2g+h$ resonance” in the unperturbed problem, defined by the
Hamiltonian $\tilde{{\mathrm{H}}}_{0}$ in (2.4) (i.e., $\alpha=0$), which is
integrable.

The resonance is the set of points in action space where the condition
$\dot{x}=0$ holds, where $x=2g+h$ is the resonant angle. This is satisfied if
the orbital inclination is equal to $i_{\star}\simeq 56.06^{\degree}$ in the
prograde case or to $i_{\star}\simeq 110.99^{\degree}$ in the retrograde case.
Indeed, according to the equations of motion, the resonance occurs when

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

that is, for all $G\neq 0$ (i.e. $e<1$) and $H=G\cos i_{\star}$ (see (2.2))
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^{\degree}$,

  * •

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

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$ is constant for $i=i_{\star}$,
while it circulates for $\left|i-i_{\star}\right|>0$. 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” rotation of ${\mathcal{O}}(i-i_{\star})$.

###  3.1 Slow-fast coordinates $(y,x)$

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

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

and the symplectic form
${\mathrm{d}}x\wedge{\mathrm{d}}y+{\mathrm{d}}h\wedge{\mathrm{d}}\Gamma.$

Notice that several symplectic transformations are possible, however we prefer
the one such that the resonant action $y$ does not depend on the
inclination777Even if we focus our analysis on the $2g+h$ resonance, our
approach is fairly general. A similar change of coordinates can be considered
for a different resonance. What is important is to choose the angle $x$ as the
one corresponding to the slow dynamics. The other changes of variables are
chosen so that the transformation is symplectic.. Hence, the action-angle
variables $(y,x)$ are associated with the variation in eccentricity.

In slow-fast variables, the Hamiltonian of the full problem 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}}\end{split}$

and

$\begin{split}{\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}$
(3.2)

See Appendix A.1 and, in particular, equation (A.1) for a detailed expansion
of the functions involved in the expression of the Hamiltonian
${\mathrm{H}}_{1}$.

Notice that, in the unperturbed problem ($\alpha=0$), $\Gamma$ and $y$ are
first integrals of the problem, because the Hamiltonian does not depend on the
angles $h$ and $x$, while in the full problem ($\alpha>0$), the phase space is
no more integrable. Moreover, in the integrable case once fixed $a$ and $e$,
the dynamics can be catalogued as a function of the inclination.

In the coordinates just defined, the $2g+h$–resonance becomes the
$x$–resonance. In the unperturbed problem, the action space associated with
the resonance is defined by the condition

$\dot{x}=\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.$

In other words, the resonance can take place on the two lines

$\begin{split}\Gamma=y\frac{-4+\sqrt{21}}{5}\quad\text{and}\quad\Gamma=y\frac{-4-\sqrt{21}}{5},\end{split}$
(3.3)

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

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

The next step is to study the dynamics in the neighborhood of a circular
orbit, that is for small values of $e>0$. However, the slow-fast variables
derived from the Delaunay variables (as happens with the original Delaunay
variables) are singular at $e=0$. In order to overcome this difficulty, we
introduce the set of 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),$
(3.4)

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, that is when
$e\approx 0$.

In this set of coordinates, the Hamiltonian of the full problem can be written
as

$\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}})$
(3.5)

where

$\displaystyle\mathcal{H}_{0}(\eta,\Gamma,\xi)$
$\displaystyle=({\mathrm{H}}_{0}\circ\Upsilon^{\mathrm{Poi}})(\eta,\Gamma,\xi)$
(3.6)
$\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}}$

and

$\mathcal{H}_{1}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})=({\mathrm{H}}_{1}\circ\Upsilon^{\mathrm{Poi}})(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}}).$
(3.7)

See Appendix A.2 and, in particular, equation (A.2) for a detailed expansion
of the functions involved in the expression of Hamiltonian $\mathcal{H}_{1}$.

##  4 The hierarchy of models and dynamical systems tools

We devote this section to explain the geometric framework that will lead to
the drifting orbits. To this end, it is convenient to autonomize the
Hamiltonian. Let us recall equation (2.1) for the time variation of
$\Omega_{\mathrm{M}}\in\mathbb{T}$, namely,

$\Omega_{\mathrm{M}}(t)=\Omega_{{\mathrm{M}},0}+n_{\Omega_{\mathrm{M}}}t.$

We introduce $\Omega_{\mathrm{M}}$ as a variable and define its symplectic
conjugate variable $J$. Then, we define the 3-degree-of-freedom Hamiltonian

$\mathcal{K}(\eta,\Gamma,J,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})=\mathcal{H}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})+n_{\Omega_{\mathrm{M}}}J,$
(4.1)

where $\mathcal{H}$ is the Hamiltonian introduced in (3.5).

We analyze this Hamiltonian in two steps, relying on the smallness of
$i_{\mathrm{M}}$. In Section 4.1, we explain the strategy of this analysis. In
this section, we also introduce the $h$-averaged Hamiltonian. We _will not
rely_ on it to construct our diffusion mechanism. However, we introduce it
since it has been widely studied numerically in literature (see, e.g., [11]).
and is a convenient simplified model to use as a first step in certain
numerical studies of the full problem. In the paper [1], we perform a detailed
analytic study of the $h-$averaged Hamiltonian, and we describe its
equilibrium points and their stability.

In Section 4.2, we introduce the two main tools that will create the “highway”
of unstable orbits: a normally hyperbolic invariant cylinder and the
associated scattering maps.

###  4.1 The hierarchy of models

To construct the 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 ratio between the semi-major axis of
the satellite and the one of the Moon $\alpha$.

In our analysis, we consider $i_{\mathrm{M}}$ arbitrarily small - the
perturbative parameter - and a realistic value for $\alpha$ (see Theorem 2.3).
Still, the smallness of $\alpha$ creates different timescales that may be
taken into account.

##### Main reduction: the Coplanar Model

Since we are assuming that the inclination of the Moon is small, the main
reduction that one can do to have an intermediate model is to take
$i_{\mathrm{M}}=0$. We refer to this model as the Coplanar Model since it
corresponds to assuming that the orbit of the Moon is coplanar to that of the
Earth. When doing this reduction, the Hamiltonian $\mathcal{K}$ in (4.1)
becomes $\Omega_{\mathrm{M}}$ independent, see Remark 2.1.

Starting from (4.1), we define

$\begin{split}\mathcal{K}_{\mathrm{CP}}(\eta,\Gamma,J,\xi,h)=\mathcal{K}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};0)=\mathcal{H}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};0)+n_{\Omega_{\mathrm{M}}}J,\end{split}$
(4.2)

where the subscript $\mathrm{CP}$ stands for coplanar. Since
$\mathcal{K}_{\mathrm{CP}}$ is $\Omega_{\mathrm{M}}$-independent, $J$ is a
first integral. In fact, one can work in the reduced phase space and consider
the 2-degree-of-freedom Hamiltonian

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

This Hamiltonian 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 (3.6) and
$\mathcal{H}_{\mathrm{CP},1}$ is the Hamiltonian $\mathcal{H}_{1}$ in (3.7)
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).$

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

##### A possible further reduction: the $h$-averaged problem

The departing point of our analysis is the coplanar Hamiltonian in (4.3)
above. However, in literature extensive investigations have been based on the
$h$-averaged coplanar model. The reduction to the $h$-averaged model is based
on the fact that, if $\alpha$ is a small parameter, the autonomous Hamiltonian
$\mathcal{H}_{\mathrm{CP}}$ has a timescale separation between the slow and
fast angles, respectively $x$ and $h$. Indeed,

$\dot{\xi},\dot{\eta}\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 another one in which the fast oscillations have
been removed by averaging over 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.$

We refer to this Hamiltonian as the $h$-averaged Hamiltonian. Note that since
the Hamiltonian $\mathcal{H}_{0}$ in (3.6) is $h$-independent,
$\mathcal{H}_{\mathrm{AV}}$ can be written as

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

Notice that $\mathcal{H}_{\mathrm{AV}}$ is a 1-degree-of-freedom Hamiltonian
provided that $\Gamma$ is a first integral of the $h$-averaged system.

####  4.1.1 Theoretical results for the coplanar and $h$-averaged system

The mechanism we use to show the existence of drifting orbits is the Arnold
diffusion mechanism (see Section 4.2 and 4.3 below). A key point in the
analysis is to study the relative position between the invariant manifolds of
the hyperbolic periodic orbits of the coplanar model.

In this work, we use circular periodic orbits, that is, the periodic orbits
located at $(\eta,\xi)=(0,0)$ (see (3.4)). Indeed, let us recall that
$\mathcal{H}_{\mathrm{CP}}(\eta,\Gamma,\xi,h)$ is a $2$-degrees-of-freedom
autonomous Hamiltonian that can be expressed as
$\mathcal{H}_{\mathrm{CP}}=\mathcal{H}_{0}+\alpha^{3}\mathcal{H}_{\mathrm{CP},1}$
with $\mathcal{H}_{0}$ as given in (3.6) and $\mathcal{H}_{\mathrm{CP},1}$ in
(A.3). From the explicit expression of the Hamiltonian obtained in Appendix
A.2, one can easily see that it does not have linear terms with respect to
$(\eta,\xi)=(0,0)$. This implies that $(\eta,\xi)=(0,0)$ is invariant. Then,
to analyze the circular periodic orbits, it is enough to look for periodic
solutions of the $1$-degree-of-freedom Hamiltonian
$\mathcal{H}_{\mathrm{CP}}(0,\Gamma,0,h)$. In Appendix A.2 we provide formulas
for this Hamiltonian (see Lemma A.1).

In [1], we prove the existence of periodic orbits at $\\{\eta=\xi=0\\}$ when
the semi-major axis $a\in[a_{\min},a_{\max}]$ with

$a_{\max}=30000\,\mathrm{km},\qquad a_{\min}=6378.14\,\mathrm{km}\;\text{(to
avoid collision)}.$

We also introduce $L_{\min,\max}=\sqrt{\mu a_{\min,\max}}$.

######  Theorem 4.1.

For any $L\in[L_{\min},L_{\max}]$, we define

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

Then, for any energy level $E\in[\mathbf{E}_{1}(L),\mathbf{E}_{2}(L)]$, there
exists a periodic orbit of the form
$\mathcal{P}_{E}(t)=(0,{\Gamma}_{E}(t),0,h_{E}(t))$ such that, $h_{E}(0)=0$
and, for $t\geq 0$, $\dot{h}_{E}(t)\neq 0$, $\Gamma_{E}(t)\in[0,0.49L]$ and

$\mathcal{H}_{\mathrm{CP}}(0,{\Gamma}_{E}(t),0,h_{E}(t))=E.$

In addition, for the case of Galileo, namely $L=1$ in non-dimensional units
(which corresponds to $a=29600$ km see Remark 2.2), one has that

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

This result ensures the existence of periodic orbits but does not give
information about the character of them. A first theoretical approach for
solving this problem is to study the character of the equilibrium point
$(\eta,\xi)=(0,0)$ for the $h$-averaged system and to consider the coplanar
Hamiltonian as a perturbation with respect to $\alpha$ of
$\mathcal{H}_{\mathrm{AV}}$. Since in this work we deal with a realistic value
of $\alpha=0.077$ (see (2.12)), this perturbative approach is not useful in
our case. However it gives an insight about the scenario we can encounter. In
[1], the following result about $\mathcal{H}_{\mathrm{AV}}$ is proven.

######  Theorem 4.2.

There exist two functions
${\Gamma}_{1},{\Gamma}_{2}:[L_{\min},L_{\max}]\to\left(0,\frac{L}{2}\right)$
such that, for $L\in[L_{\min},L_{\max}]$:

  * •

If either $\Gamma\in(0,\Gamma_{1}(L))$ or
$\Gamma\in\left(\Gamma_{2}(L),\frac{L}{2}\right)$, then $(0,0)$ is a center of
$\mathcal{H}_{\mathrm{AV}}$.

  * •

If $\Gamma\in(\Gamma_{1}(L),\Gamma_{2}(L))$, then $(0,0)$ is a saddle of
$\mathcal{H}_{\mathrm{AV}}$.

  * •

If $\Gamma=\Gamma_{1}(L)$ or $\Gamma=\Gamma_{2}(L)$, the origin is a
degenerated equilibrium point of $\mathcal{H}_{\mathrm{AV}}$.

In addition, $L^{-1}\Gamma_{1}(L)\in\left(0,\frac{m}{2}\right)$ and
$L^{-1}\Gamma_{2}(L)\in\left(\frac{m}{2},\frac{1}{2}\right)$, with
$m:=\frac{-4+\sqrt{21}}{5}$ (i.e. the slope of the prograde resonance line in
(3.3)).

For a given $L$, the values of $\Gamma_{1,2}$ can be numerically computed (as
the zeroes of some appropriate function). For instance, for the case of
Galileo (i.e., $L=1$),

$\Gamma_{1}(1)=0.029613649805289,\qquad\Gamma_{2}(1)=0.084971418151141,$

and the corresponding energy levels (for the $h$-averaged system) are

$\mathbf{E}_{1}^{\mathrm{AV}}(1)=2.072230388690642\cdot
10^{-6},\qquad\text{and}\qquad\mathbf{E}_{2}^{\mathrm{AV}}(2)=-3.473759155836634\cdot
10^{-7},$

respectively. Therefore, for $L=1$ and
$E\in(\mathbf{E}_{2}^{\mathrm{AV}}(1),\mathbf{E}_{1}^{\mathrm{AV}}(1))$, as a
consequence of Theorem 4.2, the origin is a saddle point of the $h$-averaged
system.

Applying perturbation techniques, it can be seen (see [1]) that if
$\alpha^{3}$ is small enough, the circular periodic orbits for the coplanar
Hamiltonian, $\mathcal{H}_{\mathrm{CP}}$, stated in Theorem 4.1, are of saddle
type if they belong to certain energy levels close enough to those of
$\mathcal{H}_{\mathrm{AV}}(0,\Gamma,0)$ with
$\Gamma\in(\Gamma_{1}(L),\Gamma_{2}(L))$. However, as expected, there is no
quantitative information about the maximum value of $\alpha$ such that this
result can be applied. For this reason, we have decided to consider the
coplanar (non-averaged) model to build the Arnold diffusion mechanism. This is
important to have a more realistic numerical computation of the hyperbolic
periodic orbits and the corresponding manifolds (see Section 5.1).

###  4.2 Normally hyperbolic invariant manifolds and scattering maps

We devote this section to explain the main tools that we use to construct the
“instability paths” along which the eccentricity of the satellite drifts.
These are normally hyperbolic invariant manifolds, the homoclinic
intersections of their invariant manifolds and the associated scattering maps.

Such objects can be defined for maps or flows. We will use them both for maps
and flows in Sections 5 and 6 below. Here we provide the definitions for
flows. The ones for maps are analogous.

In this section we denote by $M$ a $\mathcal{C}^{r}$ smooth manifold, by
$X\in\mathcal{C}^{r}(M,TM)$ a vector field on $M$ and by $\varphi^{t}:M\to M$
the associated smooth flow. Let $\Lambda\subset M$ be a compact
$\varphi^{t}$-invariant submanifold, possibly with boundary. By
$\varphi^{t}$-invariant we mean that $X$ is tangent to $\Lambda$, but that
orbits can escape through the boundary (a concept sometimes referred to as
_local_ or _weak_ invariance).

######  Definition 4.3.

We call $\Lambda$ a _normally hyperbolic invariant manifold_ for $\varphi^{t}$
if there is $0<\nu<\vartheta^{-1}$, a positive constant $C$ and an invariant
splitting of the tangent bundle

$T_{\Lambda}M=T\Lambda\oplus E^{\mathrm{s}}\oplus E^{\mathrm{u}}$

with ${\displaystyle T\Lambda=\bigcup_{P\in\Lambda}\big{(}\\{P\\}\times
T_{P}\Lambda\big{)}}$, ${\displaystyle
E^{\mathrm{s},\mathrm{u}}=\bigcup_{P\in\Lambda}\big{(}\\{P\\}\times
E_{P}^{\mathrm{s},\mathrm{u}}\big{)}}$ such that, for all $P\in\Lambda$,

$\begin{array}[]{ll}\|D\varphi^{t}(P)v\|\leq C\vartheta^{|t|}&\mbox{ for all
}t\in\mathbb{R},v\in T_{P}M,\\\ \|D\varphi^{t}(P)v\|\leq C\nu^{t}&\mbox{ for
all }t\geq 0,v\in E_{P}^{\mathrm{s}},\\\ \|D\varphi^{t}(P)v\|\leq
C\nu^{-t}&\mbox{ for all }t\leq 0,v\in E_{P}^{\mathrm{u}}.\end{array}$

Moreover, $\Lambda$ is called an $r$-_normally hyperbolic invariant manifold_
if it is $\mathcal{C}^{r}$ smooth and

$0<\nu<\vartheta^{-r}<1\qquad\text{for}\qquad r\geq 1.$ (4.4)

Fenichel Theory [23, 24, 25] ensures that normally hyperbolic invariant
manifolds are persistent under perturbations. Moreover, they possess stable
and unstable invariant manifolds $W^{\mathrm{s},\mathrm{u}}(\Lambda)\subset M$
defined as follows. The local stable manifold
$W^{\mathrm{s}}_{\mathrm{loc}}(\Lambda)$ is the set of points in a small
neighbourhood of $\Lambda$ whose forward orbits never leave the neighbourhood,
and tend with exponential rate to $\Lambda$. The local unstable manifold
$W^{\mathrm{u}}_{\mathrm{loc}}(\Lambda)$ is the set of points in the
neighbourhood whose backward orbits stay in the neighbourhood and tend
exponentially to $\Lambda$. We then define the stable and unstable manifold of
$\Lambda$ as

$W^{\mathrm{s}}(\Lambda)=\bigcup_{t\geq
0}^{\infty}\varphi^{-t}\left(W^{\mathrm{s}}_{\mathrm{loc}}(\Lambda)\right),\quad
W^{\mathrm{u}}(\Lambda)=\bigcup_{t\geq
0}^{\infty}\varphi^{t}\left(W^{\mathrm{u}}_{\mathrm{loc}}(\Lambda)\right).$

The stable and unstable manifolds of $\Lambda$ are foliated by what is usually
called the strong stable and strong unstable foliations, the leaves of which
we denote by $W^{\mathrm{s},\mathrm{u}}(P)$ for $P\in\Lambda$. For each
$P\in\Lambda$, the leaf $W^{\mathrm{s}}(P)$ (resp. $W^{\mathrm{u}}(P)$) of the
strong stable foliation is tangent at $P$ to $E^{\mathrm{s}}_{P}$ (resp.
$E^{\mathrm{u}}_{P}$). Moreover, the foliations satisfy that
$\varphi^{t}\left(W^{\mathrm{s}}(P)\right)=W^{\mathrm{s}}\left(\varphi^{t}(P)\right)$
and
$\varphi^{t}\left(W^{\mathrm{u}}(P)\right)=W^{\mathrm{u}}\left(\varphi^{t}(P)\right)$
for each $P\in\Lambda$ and $t\in\mathbb{R}$.

Then, one can define the usually called _wave maps_
$\pi^{\mathrm{s},\mathrm{u}}:W^{\mathrm{s},\mathrm{u}}(\Lambda)\to\Lambda$ to
be projections along leaves of the strong stable and strong unstable
foliations. That is to say, if $Q\in W^{\mathrm{s}}(\Lambda)$ then there is a
unique $Q_{+}\in\Lambda$ such that $Q\in W^{\mathrm{s}}(Q_{+})$, and so
$\pi^{\mathrm{s}}(Q)=Q_{+}$. Similarly, if $Q\in W^{\mathrm{u}}(\Lambda)$ then
there is a unique $Q_{-}\in\Lambda$ such that $Q\in W^{\mathrm{u}}(Q_{-})$, in
which case $\pi^{\mathrm{u}}(Q)=Q_{-}$. The points $Q_{\pm}$ satisfy that

$\lim_{t\to\pm\infty}\mathrm{dist}(\varphi^{t}(Q_{\pm}),\varphi^{t}(Q))=0.$

Now, suppose that $Q\in\left(W^{\mathrm{s}}(\Lambda)\pitchfork
W^{\mathrm{u}}(\Lambda)\right)\setminus\Lambda$ is a transverse homoclinic
point such that $Q\in W^{\mathrm{s}}(Q_{+})\cap W^{\mathrm{u}}(Q_{-})$. We say
that the homoclinic intersection at $Q$ is _strongly transverse_ if

$\begin{split}T_{Q}W^{\mathrm{s}}(Q_{+})\oplus
T_{Q}\left(W^{\mathrm{s}}(\Lambda)\cap
W^{\mathrm{u}}(\Lambda)\right)=T_{Q}W^{\mathrm{s}}(\Lambda),\\\
T_{Q}W^{\mathrm{u}}(Q_{-})\oplus T_{Q}\left(W^{\mathrm{s}}(\Lambda)\cap
W^{\mathrm{u}}(\Lambda)\right)=T_{Q}W^{\mathrm{u}}(\Lambda).\end{split}$ (4.5)

In this case, in a small enough neighborhood $\Xi$ of $Q$ in
$W^{\mathrm{s}}(\Lambda)\cap W^{\mathrm{u}}(\Lambda)$, (4.5) holds at each
point of $\Xi$, and the restrictions to $\Xi$ of the holonomy maps,
$\pi^{\mathrm{s},\mathrm{u}}$, are bijections onto their images. We call $\Xi$
a _homoclinic channel_. In such domain, following [16], we define the
scattering map as follows.

######  Definition 4.4.

Let $Q_{-}\in\pi^{\mathrm{u}}\left(\Xi\right)$, let
$Q=\left(\left.\pi^{\mathrm{u}}\right|_{\Xi}\right)^{-1}(Q_{-})$, and let
$Q_{+}=\pi^{\mathrm{s}}(Q)$. The _scattering map_
$S:\pi^{\mathrm{u}}(\Xi)\to\pi^{\mathrm{s}}(\Xi)$ is defined by

$S=\pi^{\mathrm{s}}\circ\left(\pi^{\mathrm{u}}\right)^{-1}:Q_{-}\longmapsto
Q_{+}.$

Assume that both $M$ (the manifold) and $X$ (the vector field) are
$\mathcal{C}^{r}$ with $r\geq 2$ and that (4.4) is satisfied, so that
$\Lambda$ is also $\mathcal{C}^{r}$. Then, the scattering map $S$ is
$\mathcal{C}^{r-1}$ (see [16]). Moreover, if the vector field has a
$\mathcal{C}^{r}$-dependence on some parameters, then the scattering map
depends on a $\mathcal{C}^{r-1}$ on them.

In general, the scattering map is only locally defined, as the transverse
homoclinic intersection of stable and unstable manifolds can be very
complicated. In the present paper, we are able to describe quite precisely the
domains of the scattering maps that we consider.

###  4.3 The Arnold diffusion instability mechanism

Once we have introduced the hierarchy of models in Section 4.1 and the
dynamical systems tools in Section 4.2, we are ready to explain the
instability mechanism that we consider to achieve drift in the eccentricity of
the satellite. Such mechanism fits into what are usually called _a priory
chaotic_ Hamiltonian systems in Arnold diffusion literature [13, 47, 22].

Consider first the coplanar Hamiltonian $\mathcal{H}_{\mathrm{CP}}$ in (4.3),
namely, $i_{\mathrm{M}}=0$. Since $\mathcal{H}_{\mathrm{CP}}$ is autonomous,
its energy is a first integral, and therefore the coordinate $J$ as well for
the Hamiltonian $\mathcal{K}_{\mathrm{CP}}$ in (4.2) (in the extended phase
space). In Section 5 we will see that at each energy level (in a certain
interval), the Hamiltonian $\mathcal{H}_{\mathrm{CP}}$ has a periodic orbit at
the $2g+h$-resonance. In Poincaré variables, this periodic orbit is located at
$\\{\eta=\xi=0\\}$. Moreover, in the considered energy interval, these
periodic orbits are hyperbolic and, thus, they have stable and unstable
invariant manifolds. In addition, we check numerically that they intersect
transversally within each energy level. The obtained homoclinic orbit has a
range of eccentricity between $e=0$ (the periodic orbit) and certain
$e=e_{\mathrm{\max}}$ which depends on the energy level888By energy level we
mean with respect to the coplanar secular Hamiltonian
$\mathcal{H}_{\mathrm{CP}}$.. For energy levels realistic for Galileo
$e_{\mathrm{\max}}\simeq 0.35$ (in particular, its orbit does not enter the
Earth’s atmosphere) whereas for higher energies $e_{\mathrm{\max}}$ reaches
0.78 (which implies that the osculating ellipse of the satellite hits the
Earth atmosphere) or higher.

\begin{overpic}[scale={0.7}]{AddingPerturbation_3D_modified.pdf}
\put(11.0,-4.0){$i_{\mathrm{M}}=0$}
\put(1.0,22.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$\Lambda_{0}$}
\put(27.0,21.0){\footnotesize\color[rgb]{0,0.78515625,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.78515625,0}
\small$W^{\mathrm{s}}(\Lambda_{0})$}
\put(-7.0,12.0){\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\small$J=J_{0}$}
\put(25.0,3.0){\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}
\footnotesize\small$W^{\mathrm{u}}(\Lambda_{0})$} \put(36.0,25.0){\tiny
homoclinic} \put(38.0,23.0){\tiny channel} \put(68.0,-4.0){$i_{\mathrm{M}}>0$}
\put(55.0,22.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$\Lambda_{i_{\mathrm{M}}}$}
\put(81.0,21.0){\color[rgb]{0,0.78515625,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.78515625,0}\footnotesize\small$W^{\mathrm{s}}(\Lambda_{i_{\mathrm{M}}})$}
\put(80.0,2.0){\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\footnotesize\small$W^{\mathrm{u}}(\Lambda_{i_{\mathrm{M}}})$}
\end{overpic} Figure 1: The normally hyperbolic invariant cylinder and its
invariant manifolds. In the left figure $i_{\mathrm{M}}=0$ and therefore $J$
is a first integral. The cylinder is foliated by invariant tori whose
invariant manifolds intersect transversally within $J=\text{constant}$ giving
rise to a homoclinic channel. In the right figure, $i_{\mathrm{M}}>0$ (small
enough). Then, $J$ is not a first integral anymore and one can construct
heteroclinic connections between different tori in the cylinder (diamond
intersections in the picture). See Corollary 5.2 and Theorem 6.1,
respectively, for a detailed description of the notation.

In the extended phase space (i.e., adding the pair $(J,\Omega_{\mathrm{M}})$
and considering the Hamiltonian $\mathcal{K}_{\mathrm{CP}}$ in (4.2)) these
periodic orbits become 2-dimensional tori. If one considers the union of such
tori for an energy interval, one has a normally hyperbolic invariant cylinder
(see Definition 4.3). Analogously, the union of the stable and unstable
manifolds of the periodic orbits (now tori) constitute the invariant manifolds
of the cylinder and their transverse intersections give rise to homoclinic
channels (see Figure 1-left). This allows to define scattering maps (in
suitable domains, detailed in Section 5).

The fact that $J$ is constant when $i_{\mathrm{M}}=0$, implies that solutions
close to the invariant manifolds of the cylinder have a constrained drift in
eccentricity: it ranges from $e=0$ (close to the cylinder) to close to the
already mentioned $e_{\mathrm{max}}$ (close to the “first” homoclinic point),
which depends on $J$.

When $i_{\mathrm{M}}>0$ but still small enough the dynamics changes
drastically. Indeed, in Section 6, we construct solutions of the full
Hamiltonian $\mathcal{K}$ in (4.1) such that the action $J$ undergoes big
changes. Equivalently, we construct solutions of the Hamiltonian $\mathcal{H}$
in (3.5) with energy drift. These solutions travel along the resonance (3.3)
and therefore, while increasing energy along the resonance, they perform
larger homoclinic excursions which lead to larger oscillations in
eccentricity.

The achievement of this drift is given by an Arnold diffusion mechanism, which
relies on the normally hyperbolic invariant manifold and the associated
invariant manifolds and scattering maps (see Section 4.2). Indeed, classical
Fenichel Theory (see [23]) implies that the normally hyperbolic invariant
cylinder of the Hamiltonian $\mathcal{H}_{\mathrm{CP}}$ in (4.3) persists for
$i_{\mathrm{M}}$ small enough and the same happens for the associated stable
and unstable invariant manifolds (see Figure 1-right). Moreover, all the
objects are regular with respect to $i_{\mathrm{M}}$. This implies that the
transverse intersections present in the case $i_{\mathrm{M}}=0$ are also
persistent and the associated scattering maps are also smooth with respect to
$i_{\mathrm{M}}$.

We use the scattering maps to construct what is usually called a pseudo-orbit.
That is, a sequence of concatenated “pieces of orbits”. More precisely, this
pseudo-orbit is formed by pieces of orbits in the cylinder (that is, pieces of
orbits of what is usually called the _inner dynamics_) which are connected by
heteroclinic orbits (see Figure 2-left). Recall that a point in the cylinder
$Q_{-}$ is mapped to a point $Q_{+}$ by the scattering map if there is a
heteroclinic orbit of the Hamiltonian (4.1) between them. Thus, to construct a
pseudo-orbit one has to understand how to “combine” the inner dynamics (the
flow restricted to the cylinder) and the scattering map (in occasions also
called outer map) to achieve drift in $J$.

\begin{overpic}[scale={0.7}]{PseudOrbitShadowing.pdf}
\put(-3.0,12.0){\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}
$J$}
\put(-3.0,25.0){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}$\Lambda_{i_{\mathrm{M}}}$}
\put(27.0,23.0){\color[rgb]{0,0.78515625,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0.78515625,0}\footnotesize\small$W^{\mathrm{s}}(\Lambda_{i_{\mathrm{M}}})$}
\put(27.0,2.0){\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\footnotesize\small$W^{\mathrm{u}}(\Lambda_{i_{\mathrm{M}}})$}
\put(49.0,25.0){ $\Lambda_{i_{\mathrm{M}}}$}
\put(80.0,2.0){\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\small
Shadowing orbit} \end{overpic} Figure 2: The pseudo-orbit (left) and the
shadowing orbit (right): The blue orbit shadows (follows closely) the pseudo-
orbit formed by heteroclinic orbits connecting different points in the
cylinder and pieces of cylinder orbits.

Both the inner and scattering maps are computed numerically relying on
perturbative methods (for $i_{\mathrm{M}}>0$ small enough). Thanks to the
particular form of the Hamiltonian (4.1), these maps have a rather simple
form. Indeed, note that one can reduce the dimension of the model by
considering a Poincaré map associated to the section $\\{h=0\\}$ and by
eliminating the variable $\Gamma$ through energy conservation. Then, the
cylinder becomes two dimensional and a good system of coordinates for it is
$(J,\Omega_{\mathrm{M}})$. In these coordinates, the inner map is of the form

$\mathcal{F}^{\mathrm{in}}_{i_{\mathrm{M}}}:\begin{pmatrix}J\\\
\Omega_{\mathrm{M}}\end{pmatrix}\mapsto\begin{pmatrix}J+i_{\mathrm{M}}\left(A_{1}^{+}(J)e^{i\Omega_{\mathrm{M}}}+A_{1}^{-}(J)e^{-i\Omega_{\mathrm{M}}}\right)+\mathcal{O}(i_{\mathrm{M}}^{2})\\\
\Omega_{\mathrm{M}}+n_{\Omega_{\mathrm{M}}}\mathcal{T}_{0}(J)+\mathcal{O}(i_{\mathrm{M}})\end{pmatrix},$
(4.6)

and one can define two scattering maps (denoted as primary and secondary) of
the form

$\mathcal{F}^{\mathrm{out},*}_{i_{\mathrm{M}}}:\begin{pmatrix}J\\\
\Omega_{\mathrm{M}}\end{pmatrix}\mapsto\begin{pmatrix}J+i_{\mathrm{M}}(B_{1}^{*,+}(J)e^{i\Omega_{\mathrm{M}}}+B_{1}^{*,-}(J)e^{-i\Omega_{\mathrm{M}}})+\mathcal{O}(i_{\mathrm{M}}^{2})\\\
\Omega_{\mathrm{M}}+n_{\Omega_{\mathrm{M}}}\zeta^{*}(J)+\mathcal{O}(i_{\mathrm{M}})\end{pmatrix},$
(4.7)

for $*\in\\{\mathrm{pri},\mathrm{sec}\\}$, associated to two different
homoclinic channels. We compute numerically the functions appearing in the
first orders of these maps. Those of the scattering maps,
$\mathcal{F}^{\mathrm{out},*}_{i_{\mathrm{M}}}$, are given by Melnikov-like
integrals (similar so those in [13, 21]).

Note that the scattering maps are not globally defined since the invariant
manifolds of the cylinder may have tangencies. This is the reason why we must
use two scattering maps: the union of domains of the two scattering maps
contain the whole region of the cylinder we are interested in.

Then, we construct a drifting pseudo-orbit. Note that we cannot use now the
classical results [43, 40], which ensure that such pseudo-orbits exist
provided the inner and scattering maps dynamics share no common invariant
curves, since they require globally defined maps (which is not the case for
the scattering maps in the present paper). Instead, we construct it from a
“generalized transition chain”, that is, a sequence of invariant quasi-
periodic tori of either the inner map or one of the scattering maps, which are
connected by an iteration of the other map. This is done in Section 6.

Once the pseudo-orbit is constructed, the final step is to show that there is
a true orbit of the Hamiltonian (4.1) which “shadows” (i.e. follows closely)
the pseudo-orbit (see Figure 2-right). For this last step we rely on shadowing
results developed in [30]. These results are rather flexible and can be easily
applied to our setting. Unfortunately, they do not provide time estimates.

######  Remark 4.5.

Instead of [30], we could have relied on more quantitative shadowing arguments
(see for instance [52, 47, 29, 10]) which could provide time estimates.
However, they require a more precise knowledge of the inner dynamics and,
therefore, to keep the proof simple, we do not use them. Still, they should be
applicable also to our setting and should lead to instability times
$T=T(i_{\mathrm{M}})$ satisfying

$|T|\leq\frac{C}{i_{\mathrm{M}}},$

for some constant $C>0$ independent of $i_{\mathrm{M}}>0$.

##  5 The coplanar secular Hamiltonian

In this section, we analyze the dynamics of the Hamiltonian
$\mathcal{H}_{\mathrm{CP}}$ (see (4.3)). That is, we consider the secular
Hamiltonian with no inclination of the Moon with respect to the ecliptic plane
(i.e., $i_{\mathrm{M}}=0$). In particular, we analyze certain features of its
dynamics and then we “translate” them to the Hamiltonian
$\mathcal{K}_{\mathrm{CP}}$ (see (4.2)) in the extended phase space.

Since we consider $\alpha=0.077$ as a fixed parameter999Instead of taking it
“small enough”. (see (2.12)), the application of perturbative techniques to
analyze $\mathcal{H}_{\mathrm{CP}}$ is not possible. Instead, in the following
sections, we assume certain ansätze on the Hamiltonian and verify them
numerically.

###  5.1 Periodic orbits and their invariant manifolds

Concerning the secular coplanar Hamiltonian in (4.3) we assume the following
ansatz. It concerns the existence of hyperbolic periodic orbits at circular
motions and the transverse intersections of their stable and unstable
invariant manifolds.

As we claimed in Section 4.1.1, the plane $\\{\eta=\xi=0\\}$ is invariant by
the coplanar Hamiltonian. The ansatz is related to the (circular) periodic
orbits lying on $\\{\eta=\xi=0\\}$.

######  Ansatz 5.1.

The Hamiltonian $\mathcal{H}_{\mathrm{CP}}(\eta,\Gamma,\xi,h)$ given in (4.3)
satisfies the following statements.

  1. 1.

In every energy level $E\in[E_{\min},E_{\max}]=[-2.12\cdot 10^{-7},1.36\cdot
10^{-6}]$, the periodic orbit
${\mathcal{P}}_{E}(t)=(0,\Gamma_{E}(t),0,h_{E}(t))$ provided in Theorem 4.1 is
hyperbolic. Denoting $\mathcal{T}(\mathcal{P}_{E})$ its period, we have that

$n_{\Omega_{\mathrm{M}}}\mathcal{T}(\mathcal{P}_{E})\in[3.9\pi,4.15\pi],\qquad\mathcal{H}_{\mathrm{CP}}(\mathcal{P}_{E}(t))=E.$

Both the periodic orbit and the period depend smoothly on $E$. In addition,
$\mathcal{T}(\mathcal{P}_{E})$ is a strictly increasing function with respect
to $E$.

  2. 2.

For each $E\in[E_{\min},E_{\max}]$, the invariant manifolds
$W^{\mathrm{u}}(\mathcal{P}_{E})$ and $W^{\mathrm{s}}(\mathcal{P}_{E})$
intersect transversally. Let $\mathcal{Q}_{E}(t)$ be any of these homoclinic
orbits.

  3. 3.

For $E\in[E_{\min},E_{\max}]$, the periodic orbit $\mathcal{P}_{E}(t)$ and the
homoclinic orbit $\chi_{E}(t)$ satisfy

$i(\mathcal{Q}_{E}(t)),i(\mathcal{P}_{E}(t))\in[55.70^{\degree},58.18^{\degree}],\qquad\text{for
all }\,t\in\mathbb{R},$

and

$\max_{t\in\mathbb{R}}e(\mathcal{Q}_{E}(t))\geq 0.795,$

where $i(\cdot)$ and $e(\cdot)$ denote the osculating inclination and
eccentricity respectively.

Note that this ansatz assumes the existence of transverse homoclinic orbits to
the periodic orbits but does not provide information on how these homoclinic
orbits depend on $E$. In fact, we will have to consider the intersections
between different branches of the invariant manifolds depending on the energy
levels to avoid homoclinic tangencies.

We devote the next section to verify numerically this ansatz.

####  5.1.1 Numerical verification of Ansatz 5.1

As we claimed in Section 4.1.1, the existence of the circular periodic orbits
located at $\\{\eta=\xi=0\\}$ is guaranteed by Theorem 4.1. Considering non-
dimensional units (see Remark 2.2), some examples of these periodic orbits
computed for the system given by $\mathcal{H}_{\mathrm{CP}}$ are shown in
Figure 3. Note that Galileo corresponds to $\mathcal{H}_{\mathrm{CP}}=1.7\cdot
10^{-8}$ (assuming $i=56.06^{\degree})$.

Figure 3: Examples of periodic orbits for $i_{M}=0$, non-dimensional units.
The colorbar reports the value of ${\mathcal{H}}_{\mathrm{CP}}$. Here it is
shown the behavior of a range of energies larger than the one we are
interested in.

Concerning the stability, the eigenvalues of the monodromy matrix of a given
periodic orbit are of the type
$\big{\\{}1,1,e^{\mathcal{T}\lambda},e^{-\mathcal{T}\lambda}\big{\\}}$, with
$\mathcal{T}$ the period of the given orbit, in case of orbits with hyperbolic
nature, or $\big{\\{}1,1,i\nu,-i\nu\big{\\}}$ in case of orbits with elliptic
nature. We will consider only the former case. In Figure 4, we show the
behavior of the eigenvalue greater than 1 for the hyperbolic periodic orbits,
as a function of ${\mathcal{H}}_{\mathrm{CP}}$.

Figure 4: The value of the eigenvalue greater than 1, $\lambda$, as a function
of the ${\mathcal{H}}_{\mathrm{CP}}$ for the hyperbolic periodic orbits.

In Figure 5, we show the period of the hyperbolic periodic orbits and the
value of $n_{\Omega_{\mathrm{M}}}\mathcal{T}$ as a function of
${\mathcal{H}}_{\mathrm{CP}}$. Recall that
$\Omega_{\mathrm{M}}(t)=\Omega_{{\mathrm{M}},0}+n_{\Omega_{\mathrm{M}}}t$ (see
(2.1)), therefore the value $n_{\Omega_{\mathrm{M}}}\mathcal{T}$ indicates the
ratio between the periods of the variables $h$ and $\Omega_{\mathrm{M}}$.
Notice that at ${\mathcal{H}}_{\mathrm{CP}}\approx 4.4472\cdot
10^{-7}\in[E_{\min},E_{\max}]$ (non-dimensional units), the period is such
that $n_{\Omega_{\mathrm{M}}}\mathcal{T}=4\pi$ which corresponds to the
resonance $2h-\Omega_{\mathrm{M}}$ (that is,
$\mathcal{T}=2\mathcal{T}_{\Omega_{\mathrm{M}}}$). This double resonance was
already observed in [11] and will be commented further at the end of the
paper.

Figure 5: The period (left) $\mathcal{T}$ of the hyperbolic periodic orbits as
a function of ${\mathcal{H}}_{\mathrm{CP}}$. On the right,
$n_{\mathrm{Saros}}\mathcal{T}$ (where $n_{\mathrm{Saros}}\equiv
n_{\Omega_{\mathrm{M}}}$) as a function of ${\mathcal{H}}_{\mathrm{CP}}$.

Next, we analyze the dynamics of the stable and the unstable manifolds of the
circular periodic orbits and we look for transverse intersections between them
at each energy level (Item 2 in the ansatz). Note that the Hamiltonian
$\mathcal{H}_{\mathrm{CP}}$ is reversible with respect to the involutions

$\Phi^{\mathrm{h}}(\eta,\Gamma,\xi,h)=(-\eta,\Gamma,\xi,-h)\quad\text{and}\quad{\Phi}^{\mathrm{v}}(\eta,\Gamma,\xi,h)=(\eta,\Gamma,-\xi,-h),$

where $\mathrm{h}$ and $\mathrm{v}$ stand for horizontal and vertical
respectively, referring to the projection of the symmetry axis onto the
$(\xi,\eta)$ plane, which are $\\{\eta=0,\,h=0\\}$, $\\{\eta=0,\,h=\pi\\}$ and
$\\{\xi=0,\,h=0\\}$, $\\{\xi=0,\,h=\pi\\}$ respectively. These symmetries
simplify the numerical verification of the ansatz since one can easily see
that the intersections between $W^{\mathrm{u}}(\mathcal{P}_{E})$ and
$W^{\mathrm{s}}(\mathcal{P}_{E})$ have points at these symmetry axes in all
the considered energy levels. However, these intersections may not be
transverse. Indeed, the angles between the invariant manifolds at the symmetry
axes depend analytically on the energy and they may vanish for a discrete set
of values of the energy.

In Figure 6, on the left, we show the behavior of the unstable invariant
manifolds associated with the periodic orbit on the Poincaré section
$\\{h=0\\}$, together with the corresponding eccentricity growth (colorbar).
Each curve corresponds to a different energy levels. In the same figure on the
right, we show the behavior of the invariant manifolds on the Poincaré section
$\\{h=0\\}$ (purple) and $\\{h=\pi\\}$ (green), for a same value of energy.

Figure 6: Left: the unstable manifolds on the Poincaré section $\\{h=0\\}$ for
different energy levels. The colorbar depicts the evolution of the
eccentricity along these unstable manifolds. Right: the behavior of the
unstable manifolds computed assuming as Poincaré map $\\{h=0\\}$ (purple) and
$\\{h=\pi\\}$ (green), for a given energy level.

All the computations of the hyperbolic manifolds have been performed following
the procedure explained in [22]. The first intersection between the stable and
the unstable manifolds takes place at $\eta=0$ ($\xi\neq 0$), after about 30
iterations on the Poincaré map of the fundamental domain assuming an initial
displacement of $10^{-8}$ in non-dimensional units from the periodic orbit
along the eigendirections. This choice corresponds to an error of $10^{-11}$
following [22].

To verify Item 2 in the ansatz we must show that, at every energy level, the
circular periodic orbit has a transverse homoclinic. We first consider the
homoclinic point at $\\{\eta=0,\,h=0\\}$ with $\xi>0$. The corresponding value
of $\xi$ and the maximum eccentricity growth achieved are shown in Figure 7.
The invariant manifolds are transverse at these points for all energies within
the considered energy range except for 6 values. (These values will be showed
later in the first column of Table 1, see also Figure 8).

Figure 7: Left: first point of intersection $\xi_{0}$ between negative branch
of the stable and the unstable manifolds of a given periodic orbit at the
crossing with the $\\{\eta=0,\,h=0\\}$ axis, as a function of
${\mathcal{H}}_{\mathrm{CP}}$ (non-dimensional units). Right: the
corresponding maximum eccentricity achievable along the hyperbolic manifolds
as a function of ${\mathcal{H}}_{\mathrm{CP}}$ (non-dimensional units).

For these 6 energy levels we must look for other homoclinic points and check
the corresponding transversality. Note that one cannot consider the homoclinic
points at $\\{\eta=0,\,h=0\\}$ with $\xi<0$ (see Figure 6) since they have the
same tangencies as the first branch considered as they are
$\Phi^{\mathrm{v}}$-symmetric. The first choice is to look for the homoclinic
points at the symmetry axis $\\{\eta=0,\,h=\pi\\}$ (with either $\xi>0$ or
$\xi<0$, since they are $\Phi^{\mathrm{v}}$-symmetric). However, the energy
values for which these homoclinic points develop tangencies are very close to
those of the other branch considered. This is shown in Table 1 and Figure 8,
where we show the splitting angle computed at the first point of intersection
between the negative branch of the stable and the unstable manifolds at
$\\{\eta=0,\,h=0\\}$ and $\\{\eta=0,\,h=\pi\\}$ as a function of the given
energy level. One should expect that even if they look to be very close, the
tangency energy values of the two branches are disjoint. However, to check
this would require a numerical analysis with a much higher precision.

$\mathcal{H}_{\mathrm{CP}}$ | $\Delta\mathcal{H}_{\mathrm{CP}}$ | $\Delta\mathcal{H}_{\mathrm{CP}}/\mathcal{H}_{\mathrm{CP}}$ | $\phi$
---|---|---|---
1.34294e-6 | 1.2e-11 | 8.9e-6 | 0.788
1.23642e-6 | 1.5e-10 | 1.2e-4 | 0.997
1.09175e-6 | 6.1e-12 | 5.6e-6 | 1.183
8.9030e-7 | 6.6e-12 | 7.4e-6 | 1.276
6.0662e-7 | 2.5e-11 | 4.1e-5 | 1.554
2.1005e-7 | 7.5e-11 | 3.5e-4 | 1.935
Table 1: The value of $\mathcal{H}_{\mathrm{CP}}$ corresponding to a non-
transverse intersection at $\\{\eta=0,\,h=0\\}$ and the absolute and relative
difference with respect to the value of $\mathcal{H}_{\mathrm{CP}}$
corresponding to a non-transverse intersection at $\\{\eta=0,\,h=\pi\\}$. On
the last column, the splitting angle $\phi$ (rad) computed at
$\\{\xi=0,\,h=0\\}$ at the same energy level.  Figure 8: Splitting angle
between the negative branch of the stable and the unstable manifolds computed
at $\\{\eta=0,\,h=0\\}$ (purple) and $\\{\eta=0,\,h=\pi\\}$ axis (green), as a
function of ${\mathcal{H}}_{\mathrm{CP}}$ (non-dimensional units). The
tangencies of both branches are very close to each other and beyond the
numerical accuracy that we consider.

Instead, to keep the numerical analysis “simple”, we consider a different
branch: we look for homoclinic points at the symmetry axis
$\\{\xi=0,\,h=0\\}$. We refer to these homoclinic points as secondary since
their orbits approach twice the saddle, in contrast to the first branch, which
we denote by primary. The corresponding splitting angle $\phi$ is shown in
Table 1 (last column) and an example of such intersection is given in Figure
9.

Figure 9: Homoclinic intersection (green point) computed at
$\\{\xi=0,\,h=0\\}$ for ${\mathcal{H}}_{\mathrm{CP}}=2.1005\cdot 10^{-7}$
(non-dimensional units). In red the negative branch of the unstable manifold,
in blue the negative branch of the stable manifold.

Once we have computed the periodic orbits and the associated transverse
homoclinics, Item 3 of the ansatz can be easily verified by integrating
numerically the flow of the coplanar system.

###  5.2 The invariant cylinder and their invariant manifolds

In this section we analyze how the periodic orbits provided by Ansatz 5.1
(Item 1) give rise to a normally hyperbolic invariant cylinder and, by means
of Item 2 of Ansatz 5.1, we analyze the associated stable and unstable
invariant manifolds and their intersections. Then, we study the inner dynamics
(dynamics restricted to the cylinder) and the scattering maps of the coplanar
secular model (see Section 4.2).

Our aim is to study the complete problem as a perturbation of the coplanar
one. To this end, we consider the Hamiltonian $\mathcal{K}_{\mathrm{CP}}$ (see
(4.2)) defined in the extended phase space
$(\eta,\Gamma,J,\xi,h,\Omega_{\mathrm{M}})$. In other words, we keep the
conjugated variables $(J,\Omega_{\mathrm{M}})$ even if $\Omega_{\mathrm{M}}$
is a cyclic variable (see (2.1)). As a result, the periodic orbits considered
in Ansatz 5.1 become $2$-dimensional invariant tori. Indeed, we first notice
that the energy level $\mathcal{K}_{\mathrm{CP}}=0$ in the extended phase
space corresponds to $\mathcal{H}_{\mathrm{CP}}=E$ with
$E=-n_{\Omega_{\mathrm{M}}}J$. In addition, since $\dot{h}_{E}(t)\neq 0$,
$t=h_{E}^{-1}(h)$. Then, from functions $h^{-1}_{E}(t)$ and $\Gamma_{E}(t)$,
we can obtain a parametrization of these tori in terms of
$(h,\Omega_{\mathrm{M}})$ of the form

$\Big{(}0,\Gamma_{J}^{\mathrm{ext}}(h),J,0,h,\Omega_{\mathrm{M}}\Big{)},\qquad\text{for
}(h,\Omega_{\mathrm{M}})\in\mathbb{T}^{2},$

for

$J\in[J_{\min},J_{\max}]=\left[-\frac{E_{\max}}{n_{\Omega_{\mathrm{M}}}},-\frac{E_{\min}}{n_{\Omega_{\mathrm{M}}}}\right],$

with $E_{\max}$, $E_{\min}$ as given in Ansatz 5.1. Therefore, the union of
these $2$-dimensional invariant tori form a normally hyperbolic invariant
$3$-dimensional manifold $\Lambda_{0}$, diffeomorphic to a cylinder
$\mathbb{T}^{2}\times[J_{\min},J_{\max}]$. Applying the implicit function
theorem with respect to the energy $J$, one can see that the cylinder is
analytic (by Ansatz 5.1, the periodic orbits are hyperbolic, thus non-
degenerate). We summarize this fact in the following lemma, which is a direct
consequence of Ansatz 5.1 (see also Figure 1-left).

######  Corollary 5.2.

Assume Ansatz 5.1. Then, the Hamiltonian $\mathcal{K}_{\mathrm{CP}}$ (see
(4.2)) has an analytic normally hyperbolic invariant $3$-dimensional cylinder
$\Lambda_{0}$, which is foliated by $2$-dimensional invariant tori on the
constant invariant hyperplanes $J=J_{0}\in[J_{\min},J_{\max}]$.

Moreover, the cylinder $\Lambda_{0}$ has 4-dimensional invariant manifolds,
denoted by $W^{\mathrm{u}}(\Lambda_{0})$ and $W^{\mathrm{s}}(\Lambda_{0})$,
which, for every $J\in[J_{\min},J_{\max}]$, intersect transversally either at
the symmetry axis $\\{\eta=0,h=0\\}$ with $\xi>0$ or at $\\{\xi=0,h=0\\}$ with
$\eta>0$.

From now on, we will call primary homoclinic orbits (and denote them with a
superindex $\mathrm{pri}$) those which intersect the symmetry axis
$\\{\eta=0,h=0\\}$ with $\xi>0$ and we will call secondary homoclinic orbits
(and denote them with a superindex $\mathrm{sec}$) those which intersect the
symmetry axis $\\{\xi=0,h=0\\}$ with $\eta>0$. Indeed, the first ones only get
close to the periodic orbit once while the second ones approach it twice.

To analyze the dynamics on the cylinder and the scattering map associated with
the transverse intersections of its invariant manifolds it is more convenient
to consider a Poincaré map. Therefore, we define a global Poincaré section and
work with maps to reduce the dimension by one.

There are two natural choices of section, $\\{h=0\\}$ or
$\\{\Omega_{\mathrm{M}}=0\\}$, since both variables satisfy $\dot{h}\neq 0$
(see Theorem 4.1) and $\dot{\Omega}_{\mathrm{M}}=n_{\Omega_{\mathrm{M}}}\neq
0$. It would seem more natural to choose the section
$\\{\Omega_{\mathrm{M}}=0\\}$, however since we aim to prove drift on the
energy $J$, it is more convenient to consider the section

$\Sigma=\\{h=0\\}$

and the induced Poincaré map

$\Pi_{0}:\Sigma\to\Sigma\qquad$ (5.1)

Notice that this Poincaré section was already used in the numerical study in
Section 5.1. In addition, we denote the intersection of the cylinder
$\Lambda_{0}$ with the section $\Sigma$ as

$\widetilde{\Lambda}_{0}=\Lambda_{0}\cap\Sigma.$

Notice that $\widetilde{\Lambda}_{0}$ is a 2-dimensional normally hyperbolic
manifold for the Poincaré map $\Pi_{0}$ with 3-dimensional stable and unstable
manifolds, which we denote $W^{\mathrm{u}}(\widetilde{\Lambda}_{0})$ and
$W^{\mathrm{s}}(\widetilde{\Lambda}_{0})$ for $j=1,2$. They intersect at
transverse homoclinic orbits which are just the those analyzed in Corollary
5.2 intersected with the section $\Sigma$.

To fix notation, we denote by $\mathtt{D}^{*}_{0}$, with
$*\in\\{\mathrm{pri},\mathrm{sec}\\}$, the subset of $[J_{\min},J_{\max}]$ for
which there exist primary/secondary homoclinic orbits to the corresponding
invariant torus in the cylinder. By Corollary 5.2,
$\mathtt{D}^{\mathrm{pri}}_{0}\cup\mathtt{D}_{0}^{\mathrm{sec}}=[J_{\min},J_{\max}]$.

######  Corollary 5.3.

Assume Ansatz 5.1. The Poincaré map $\Pi_{0}$ defined in (5.1) and induced by
Hamiltonian $\mathcal{K}_{\mathrm{CP}}$ (see (4.2)) has a normally hyperbolic
$2$-dimensional cylinder $\widetilde{\Lambda}_{0}$ foliated by invariant
curves. In addition, there exists an analytic function
$\mathcal{G}_{0}:[J_{\min},J_{\max}]\times\mathbb{T}\to\mathbb{R}^{3}\times\mathbb{T}^{3}$,

$\mathcal{G}_{0}(J,\Omega_{\mathrm{M}})=\big{(}0,{\mathcal{G}}_{0}^{\Gamma}(J),J,0,0,\Omega_{\mathrm{M}}\big{)},$

that parametrizes $\widetilde{\Lambda}_{0}$,

$\widetilde{\Lambda}_{0}=\\{\mathcal{G}_{0}(J,\Omega_{\mathrm{M}})\,:\,(J,\Omega_{\mathrm{M}})\in[J_{\min},J_{\max}]\times\mathbb{T}\\}.$

Moreover, for $*\in\\{\mathrm{pri},\mathrm{sec}\\}$, within the hypersurface
$J=J_{0}$ with $J_{0}\in\mathtt{D}^{*}_{0}$, the invariant manifolds
$W^{\mathrm{u}}(\widetilde{\Lambda}_{0})$ and
$W^{\mathrm{s}}(\widetilde{\Lambda}_{0})$ intersect transversally at the
symmetry axis of the involutions $\Phi^{\mathrm{h}}$ (for $*=\mathrm{pri}$)
and $\Phi^{\mathrm{v}}$ (for $*=\mathrm{sec}$).

Let $\Xi_{0}^{*}$ denote these transverse intersections on the symmetry axes.
Then, there exists an analytic function
$\mathcal{C}^{*}_{0}:\mathtt{D}^{*}_{0}\times\mathbb{T}\to\mathbb{R}^{3}\times\mathbb{T}^{3},$
such that

$\displaystyle\mathcal{C}^{\mathrm{pri}}_{0}(J,\Omega_{\mathrm{M}})$
$\displaystyle=\Big{(}0,\mathcal{C}_{0}^{\Gamma,\mathrm{pri}}(J),J,\mathcal{C}_{0}^{\xi,\mathrm{pri}}(J),0,\Omega_{\mathrm{M}}\Big{)},$
$\displaystyle\mathcal{C}^{\mathrm{sec}}_{0}(J,\Omega_{\mathrm{M}})$
$\displaystyle=\Big{(}\mathcal{C}_{0}^{\eta,\mathrm{sec}}(J),\mathcal{C}_{0}^{\Gamma,\mathrm{sec}}(J),J,0,0,\Omega_{\mathrm{M}}\Big{)},$

that parametrizes $\Xi_{0}^{*}$. That is,

$\Xi^{*}_{0}=\\{\mathcal{C}^{*}_{0}(J,\Omega_{\mathrm{M}})\,:\,(J,\Omega_{\mathrm{M}})\in[J_{\min},J_{\max}]\times\mathbb{T}\\}.$

The subscript $0$ in the previous structures and parameterizations indicates
that we are dealing with the coplanar case, i.e., $i_{\mathrm{M}}=0$.

Notice that Corollary 5.2 gives global coordinates for the cylinder
$\widetilde{\Lambda}_{0}$. Moreover, these coordinates are symplectic with
respect to the canonical symplectic form $d\Omega_{\mathrm{M}}\wedge dJ$.
Indeed, Corollary 5.3 implies that at the cylinder $\widetilde{\Lambda}_{0}$
one has $\eta=\xi=h=0$ and $\Gamma=\mathcal{G}_{0}^{\Gamma}(J)$. Then, the
pullback of the canonical form $d\xi\wedge d\eta+dh\wedge
d\Gamma+d\Omega_{\mathrm{M}}\wedge dJ$ to $\widetilde{\Lambda}_{0}$ is just
$d\Omega_{\mathrm{M}}\wedge dJ$.

Next, we consider the inner and the two scattering maps (for the primary and
secondary homoclinic channels). The inner map is defined in the whole cylinder
$\widetilde{\Lambda}_{0}$ whereas the scattering maps are defined in open
domains of the cylinder where the associated homoclinic channels are
transverse. Since $J$ is conserved by the inner and scattering maps, these
maps are integrable and the variables $(J,\Omega_{\mathrm{M}})$ are action-
angle coordinates. In these variables, it will be easier to later understand
the effect of the inclination $i_{\mathrm{M}}$ on the inner and scattering
maps.

####  5.2.1 The coplanar inner map

In this section we study the inner map, as introduced in (4.6), restricted to
the normally hyperbolic manifold $\widetilde{\Lambda}_{0}$. In particular, the
inner map
$\mathcal{F}_{0}^{\mathrm{in}}:\widetilde{\Lambda}_{0}\to\widetilde{\Lambda}_{0}$
is defined as the Poincaré map $\Pi_{0}$ in (5.1) restricted to the symplectic
invariant manifold $\widetilde{\Lambda}_{0}$. We express
$\mathcal{F}_{0}^{\mathrm{in}}$ using the global coordinates
$(J,\Omega_{\mathrm{M}})$ of $\widetilde{\Lambda}_{0}$. Since $J$ is an
integral of motion and $\dot{\Omega}_{{\mathrm{M}}}=n_{\Omega_{\mathrm{M}}}$,
the inner map has the form of

$\mathcal{F}_{0}^{\mathrm{in}}:\begin{pmatrix}J\\\
\Omega_{\mathrm{M}}\end{pmatrix}\mapsto\begin{pmatrix}J\\\
\Omega_{\mathrm{M}}+n_{\Omega_{\mathrm{M}}}\mathcal{T}_{0}(J)\end{pmatrix},$
(5.2)

where $\mathcal{T}_{0}(J)$ is the period of the periodic orbit obtained in
Ansatz 5.1 on the corresponding energy surface $E=-n_{\Omega_{\mathrm{M}}}J$.

The following lemma is a direct consequence of Ansatz 5.1 (see Figure 5).

######  Lemma 5.4.

Assume Ansatz 5.1. The analytic symplectic inner map
$\mathcal{F}_{0}^{\mathrm{in}}$ defined in (5.2) is twist, that is

$\partial_{J}\mathcal{T}_{0}(J)<0\quad\text{for}\quad
J\in[J_{\min},J_{\max}].$

In addition, there exists a unique $J_{\mathrm{res}}\in[J_{\min},J_{\max}]$
such that

$n_{\Omega_{\mathrm{M}}}\mathcal{T}_{0}(J_{\mathrm{res}})=4\pi.$

Moreover,

$n_{\Omega_{\mathrm{M}}}\mathcal{T}_{0}(J)\not\in
2\pi\mathbb{Z}\qquad\text{for all }\quad
J\in[J_{\min},J_{\max}]\setminus\\{J_{\mathrm{res}}\\}.$

This lemma is crucial to later, in Section 6.3, construct the heteroclinic
connections between different tori of the perturbed cylinder (with
$i_{\mathrm{M}}>0$).

####  5.2.2 The coplanar scattering maps

In this section we study the scattering maps, as introduced in (4.7),
associated to the Poincaré map $\Pi_{0}$ which sends $\widetilde{\Lambda}_{0}$
to itself by means of the homoclinic channels $\Xi_{0}^{\mathrm{pri}}$ and
$\Xi_{0}^{\mathrm{sec}}$ (see Corollary 5.3). For each homoclinic channel, we
denote

$\mathcal{F}^{\mathrm{out},\mathrm{pri}}_{0}:\widetilde{\Lambda}_{0}\to\widetilde{\Lambda}_{0},\qquad\mathcal{F}^{\mathrm{out},\mathrm{sec}}_{0}:\widetilde{\Lambda}_{0}\to\widetilde{\Lambda}_{0}.$

Notice the abuse of notation since the scattering maps are only defined
provided $J\in\mathtt{D}_{0}^{\mathrm{pri}}$ or
$J\in\mathtt{D}_{0}^{\mathrm{sec}}$, respectively, and not in the whole
cylinder $\widetilde{\Lambda}_{0}$.

The scattering map is always exact symplectic, see [16]. So, since $J$ is
preserved, the scattering maps must be of the form

$\mathcal{F}^{\mathrm{out},*}_{0}:\begin{pmatrix}J\\\
\Omega_{\mathrm{M}}\end{pmatrix}\mapsto\begin{pmatrix}J\\\
\Omega_{\mathrm{M}}+n_{\Omega_{\mathrm{M}}}\zeta^{*}(J)\end{pmatrix},\quad\text{for}\quad*\in\\{\mathrm{pri},\mathrm{sec}\\}.$
(5.3)

The functions $\zeta^{*}$ are usually called phase shift. Indeed the
homoclinic orbits in $\Xi^{*}_{0}$ are homoclinic to a periodic orbit.
However, they are asymptotic to different trajectories (that is, different
phases) in the periodic orbit.

To compute the phase shifts $\zeta^{*}$, it is convenient to deal with a flow
instead of the associated Poincaré map, because then one can rely on Melnikov
Theory. However, the outer map induced by the flow associated to the
Hamiltonian (4.2) does not preserve the section $\\{h=0\\}$ and, therefore,
one cannot deduce the outer map associated with the Poincaré map $\Pi_{0}$
from that of the flow. Then, we reparameterize the flow so that its scattering
maps preserve this section. This reparameterization corresponds to identifying
the variable $h$ with time $t$ and is given by

$\displaystyle\frac{d\xi}{ds}$
$\displaystyle=\frac{\partial_{\eta}\mathcal{H}_{\mathrm{CP}}(\eta,\Gamma,\xi,h)}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}(\eta,\Gamma,\xi,h)},$
$\displaystyle\frac{d\eta}{ds}$
$\displaystyle=-\frac{\partial_{\xi}\mathcal{H}_{\mathrm{CP}}(\eta,\Gamma,\xi,h)}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}(\eta,\Gamma,\xi,h)},$
(5.4) $\displaystyle\frac{dh}{ds}$ $\displaystyle=1,$
$\displaystyle\frac{d\Gamma}{ds}$
$\displaystyle=-\frac{\partial_{h}\mathcal{H}_{\mathrm{CP}}(\eta,\Gamma,\xi,h)}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}(\eta,\Gamma,\xi,h)},$
$\displaystyle\frac{d\Omega_{\mathrm{M}}}{ds}$
$\displaystyle=\frac{n_{\Omega_{\mathrm{M}}}}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}(\eta,\Gamma,\xi,h)},$
$\displaystyle\frac{dJ}{ds}$ $\displaystyle=0.$

Note that the right hand side of the system of equations (5.4) does not depend
on $\Omega_{\mathrm{M}}$. Let
$\widehat{\varphi}_{\mathrm{CP}}\\{s,(\eta,\Gamma,J,\xi,h,\Omega_{\mathrm{M}})\\}$
be the flow associated to equation (5.4) and
$\varphi_{\mathrm{CP}}\\{s,(\eta,\Gamma,\xi,h)\\}$ be the flow associated with
its $(\eta,\Gamma,\xi,h)$ components. Componentwise, for the reduced and
extended phase space, it can be written as

$\displaystyle\varphi_{\mathrm{CP}}\\{s,(\eta,\Gamma,\xi,h)\\}=\big{(}$
$\displaystyle\varphi^{\eta}\\{s,(\eta,\Gamma,\xi,h)\\},\varphi^{\Gamma}\\{s,(\eta,\Gamma,\xi,h)\\},$
(5.5) $\displaystyle\varphi^{\xi}\\{s,(\eta,\Gamma,\xi,h)\\},h+s\big{)},$
$\displaystyle\widehat{\varphi}_{\mathrm{CP}}\\{s,(\eta,\Gamma,J,\xi,h,\Omega_{\mathrm{M}})\\}=\big{(}$
$\displaystyle\varphi^{\eta}\\{s,(\eta,\Gamma,\xi,h)\\},\varphi^{\Gamma}\\{s,(\eta,\Gamma,\xi,h)\\},J,$
$\displaystyle\varphi^{\xi}\\{s,(\eta,\Gamma,\xi,h)\\},h+s,\Omega_{\mathrm{M}}+\varphi^{\Omega_{\mathrm{M}}}\\{s,(\eta,\Gamma,\xi,h)\\}\big{)}.$

In addition, we denote the trajectories departing from points either in the
normally hyperbolic cylinder or in one of the homoclinic channels,
respectively, as

$\displaystyle\gamma_{J}(s)$
$\displaystyle=\varphi_{\mathrm{CP}}\\{s,(0,\mathcal{G}_{0}^{\Gamma}(J),0,0)\\},$
$\displaystyle\chi^{\mathrm{pri}}_{J}(s)$
$\displaystyle=\varphi_{\mathrm{CP}}\\{s,(0,\mathcal{C}_{0}^{\Gamma,\mathrm{pri}}(J),\mathcal{C}_{0}^{\xi,\mathrm{pri}}(J),0)\\},$
(5.6) $\displaystyle\chi^{\mathrm{sec}}_{J}(s)$
$\displaystyle=\varphi_{\mathrm{CP}}\\{s,(\mathcal{C}_{0}^{\eta,\mathrm{sec}}(J),\mathcal{C}_{0}^{\Gamma,\mathrm{sec}}(J),0,0)\\}.$

######  Lemma 5.5.

Assume Ansatz 5.1. The functions $\zeta^{*}(J)$ involved in the definition of
the outer maps in (5.3) are given by

$\zeta^{*}(J)=\zeta^{*}_{+}(J)-\zeta^{*}_{-}(J).$

where

$\zeta^{*}_{+}(J)=-\zeta^{*}_{-}(J)=\lim_{N\to+\infty}\sum_{l=0}^{N-1}\left(\int_{2\pi
l}^{2\pi(l+1)}\frac{ds}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}\circ\chi^{*}_{J}(s)}-\mathcal{T}_{0}(J)\right).$

In Figure 10-left, we plot the function $\zeta^{\mathrm{pri}}(J)$ for the case
under study. This figure prompts us to assume the following ansatz, which will
be used in Section 6.3. Roughly speaking it will ensure that the inner and
outer maps do not have simultaneous resonances at their first orders in
$i_{\mathrm{M}}$.

######  Ansatz 5.6.

For $J=J_{\mathrm{res}}$, the value introduced in Lemma 5.4, the function
$\zeta^{\mathrm{pri}}$ satisfies that

$\zeta^{\mathrm{pri}}(J_{\mathrm{res}})\not\in\pi\mathbb{Z}\qquad\text{and}\qquad\partial_{J}{\zeta^{\mathrm{pri}}}(J_{\mathrm{res}})\neq
0,$

Figure 10: As a function of ${\mathcal{H}}_{\mathrm{CP}}\equiv-
n_{\Omega_{\mathrm{M}}}J$, we show on the left
$n_{\mathrm{Saros}}\zeta^{\mathrm{pri}}(J)$. Here $n_{\mathrm{Saros}}\equiv
n_{\Omega_{\mathrm{M}}}$. Note that at
${\mathcal{H}}_{\mathrm{CP}}=4.4472\times 10^{-7}$ corresponding to the
existence of a secondary resonance (see Fig. 5 on the right),
$n_{\mathrm{Saros}}\zeta^{\mathrm{pri}}(J)/\pi\approx 5.61$. On the right,
$n_{\mathrm{Saros}}\zeta^{\mathrm{sec}}(J)/2\pi$ for the non-transverse values
(see Table 1). The integral of Lemma 5.5 has been computed by means of the
function qags of the quadpack Fortran package.

Note that in Figure 10-left we only depict the horizontal lines at
$\zeta^{\mathrm{pri}}=2\pi k$ instead of $\zeta=\pi k$ to have a clear
picture. In Figure 10-right, we plot the functions $\zeta^{\mathrm{sec}}(J)$
for the values of $J\in[J_{\min},J_{\max}]$ where the primary homoclinic
channel is not transverse (see Table 1).

###### Proof of Lemma 5.5.

Since the flow component
$(\widehat{\varphi}_{\mathrm{CP}})^{\Omega_{\mathrm{M}}}$ is of the form
$\Omega_{\mathrm{M}}+\varphi^{\Omega_{\mathrm{M}}}$ with
$\varphi^{\Omega_{\mathrm{M}}}$ independent of $\Omega_{\mathrm{M}}$, its
behavior is given by

$\varphi^{\Omega_{\mathrm{M}}}\\{s,(\eta,\Gamma,\xi,h)\\}=n_{\Omega_{\mathrm{M}}}\int_{0}^{s}\frac{ds}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}\circ\varphi_{\mathrm{CP}}\\{s,(\eta,\Gamma,\xi,h)\\}}.$
(5.7)

First, we obtain an integral expression for the period $\mathcal{T}_{0}(J)$
given in (5.2). Indeed, since the inner map is just the $2\pi$-time map of the
flow $\widehat{\varphi}_{\mathrm{CP}}$ for the components
$(J,\Omega_{\mathrm{M}})$, by Corollary 5.3 one has that,

$\mathcal{T}_{0}(J)=\varphi^{\Omega_{\mathrm{M}}}\\{2\pi,(0,\mathcal{G}_{0}^{\Gamma}(J),0,0)\\}$

and, by (5.6) and (5.7),

$\displaystyle\mathcal{T}_{0}(J)$
$\displaystyle=\int_{0}^{2\pi}\frac{ds}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}\circ\gamma_{J}(s)}.$
(5.8)

Let us now consider the primary scattering map. One has the homoclinic point

$\mathcal{C}_{0}^{\mathrm{pri}}(J,\Omega_{\mathrm{M}})=(0,\mathcal{C}_{0}^{\Gamma,\mathrm{pri}}(J),J,\mathcal{C}_{0}^{\xi,\mathrm{pri}}(J),0,\Omega_{\mathrm{M}})\in
W^{\mathrm{u},1}({\Lambda}_{0})\cap
W^{\mathrm{s},1}({\Lambda}_{0})\cap\\{h=0\\},$

(see Corollaries 5.2 and 5.3). Since the $(\eta,\Gamma,\xi,h)$ components are
independent of $\Omega_{\mathrm{M}}$, this point is forward asymptotic (in the
reparametrized time) to a point

$\displaystyle
Q_{+}=\big{(}0,\mathcal{G}_{0}^{\Gamma}(J),J,0,0,\Omega_{\mathrm{M}}+n_{\Omega_{\mathrm{M}}}\zeta^{\mathrm{pri}}_{+}(J)\big{)},$

and backward asymptotic to a point

$\displaystyle
Q_{-}=\big{(}0,\mathcal{G}_{0}^{\Gamma}(J),J,0,0,\Omega_{\mathrm{M}}+n_{\Omega_{\mathrm{M}}}\zeta^{\mathrm{pri}}_{-}(J)\big{)}.$

The scattering map corresponds to the application $Q_{-}\mapsto Q_{+}$ (see
Definition 4.4). Then, by the expression of the map
$\mathcal{F}^{\mathrm{out},\mathrm{pri}}_{0}$ in (5.3), one has that

$\zeta^{\mathrm{pri}}(J)=\zeta^{\mathrm{pri}}_{+}(J)-\zeta^{\mathrm{pri}}_{-}(J).$

By the definition of $\zeta_{+}^{\mathrm{pri}}(J)$, using definitions (5.6) of
$\chi_{J}^{\mathrm{pri}},\gamma_{J}$ and applying (5.7), one obtains

$\displaystyle\zeta^{\mathrm{pri}}_{+}(J)$
$\displaystyle=\frac{1}{n_{\Omega_{\mathrm{M}}}}\lim_{s\to+\infty}\left(\varphi^{\Omega_{\mathrm{M}}}\\{s,(0,\mathcal{C}_{0}^{\Gamma,\mathrm{pri}}(J),\mathcal{C}_{0}^{\xi,\mathrm{pri}},0)\\}-\varphi^{\Omega_{\mathrm{M}}}\\{s,(0,\mathcal{G}_{0}^{\Gamma}(J),0,0)\\}\right)$
(5.9)
$\displaystyle=\lim_{s\to+\infty}\int_{0}^{s}\left(\frac{d\sigma}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}\circ\chi^{\mathrm{pri}}_{J}(\sigma)}-\frac{d\sigma}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}\circ\gamma_{J}(\sigma)}\right).$

Since system (5.4) is $2\pi$-periodic in $s$, due to the identification of $s$
with $h$, it is more convenient to write these integrals as

$\displaystyle\zeta^{\mathrm{pri}}_{+}(J)$
$\displaystyle=\lim_{N\to+\infty}\int_{0}^{2\pi
N}\left(\frac{ds}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}\circ\chi^{\mathrm{pri}}_{J}(s)}-\frac{ds}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}\circ\gamma_{J}(s)}\right),$

and, taking into account (5.8), one has that

$\displaystyle\zeta^{\mathrm{pri}}_{+}(J)$
$\displaystyle=\lim_{N\to+\infty}\sum_{l=0}^{N-1}\int_{2\pi
l}^{2\pi(l+1)}\left(\frac{ds}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}\circ\chi^{\mathrm{pri}}_{J}(s)}-\mathcal{T}_{0}(J)\right).$

Let us recall that $\mathcal{H}_{\mathrm{CP}}$ is reversible with respect to
$\Phi^{\mathrm{h}}(\eta,\Gamma,\xi,h)=(-\eta,\Gamma,\xi,-h)$ (see (5.1.1)).
Therefore, the flow satisfies
$\varphi_{\mathrm{CP}}\\{s,(0,\Gamma,\xi,0)\\}=\Phi^{\mathrm{h}}\circ\varphi_{\mathrm{CP}}\\{-s,(0,\Gamma,\xi,0)\\}$.
Then, one can see that
$\varphi^{\Omega_{\mathrm{M}}}\\{s,(0,\Gamma,\xi,0)\\}=-\varphi^{\Omega_{\mathrm{M}}}\\{-s,(0,\Gamma,\xi,0)\\}$.
As a result,

$\displaystyle\zeta^{\mathrm{pri}}_{-}(J)$
$\displaystyle=\frac{1}{n_{\Omega_{\mathrm{M}}}}\lim_{s\to-\infty}\left(\varphi^{\Omega_{\mathrm{M}}}\\{s,(0,\mathcal{C}_{0}^{\Gamma,\mathrm{pri}}(J),\mathcal{C}_{0}^{\xi,\mathrm{pri}},0)\\}-\varphi^{\Omega_{\mathrm{M}}}\\{s,(0,\mathcal{G}_{0}^{\Gamma}(J),0,0)\\}\right)$
(5.10) $\displaystyle=-\zeta^{\mathrm{pri}}_{+}(J).$

We proceed analogously for the secondary scattering map. ∎

##  6 Dynamics of the system with $i_{\mathrm{M}}>0$

In this section, we consider the Hamiltonian $\mathcal{K}$ in (4.1) and assume
$0<i_{\mathrm{M}}\ll 1$. Since we want to compare its dynamics with that of
$\mathcal{K}_{\mathrm{CP}}$ (see also (4.2)), we write $\mathcal{K}$ as

$\displaystyle\mathcal{K}(\eta,\Gamma,J,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})$
$\displaystyle=\mathcal{H}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})+n_{\Omega_{\mathrm{M}}}J$
(6.1)
$\displaystyle=\mathcal{H}_{\mathrm{CP}}(\eta,\Gamma,\xi,h)+\alpha^{3}i_{\mathrm{M}}\mathcal{R}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})+n_{\Omega_{\mathrm{M}}}J,$

where $\mathcal{H}_{\mathrm{CP}}$ has been introduced in (4.3). See Appendix
A.2 for the whole expression of $\mathcal{R}$.

###  6.1 Perturbative analysis of the inner and scattering maps

The first step is to compute perturbatively the inner and outer maps. To this
end, we apply Poincaré-Melnikov techniques. As done in the previous section
(see (5.4)), we reparameterize the flow so that the variable $h$ becomes
“time”. Namely,

$\displaystyle\frac{d\xi}{ds}$
$\displaystyle=\frac{\partial_{\eta}\mathcal{H}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})}{\partial_{\Gamma}\mathcal{H}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})},$
$\displaystyle\frac{d\eta}{ds}$
$\displaystyle=-\frac{\partial_{\xi}\mathcal{H}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})}{\partial_{\Gamma}\mathcal{H}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})},$
(6.2) $\displaystyle\frac{dh}{ds}$ $\displaystyle=1,$
$\displaystyle\frac{d\Gamma}{ds}$
$\displaystyle=-\frac{\partial_{h}\mathcal{H}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})}{\partial_{\Gamma}\mathcal{H}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})},$
$\displaystyle\frac{d\Omega_{\mathrm{M}}}{ds}$
$\displaystyle=\frac{n_{\Omega_{\mathrm{M}}}}{\partial_{\Gamma}\mathcal{H}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})}$
$\displaystyle\frac{dJ}{ds}$
$\displaystyle=-\alpha^{3}i_{\mathrm{M}}\frac{\partial_{\Omega_{\mathrm{M}}}\mathcal{R}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})}{\partial_{\Gamma}\mathcal{H}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})}.$

This system is $\mathcal{O}(i_{\mathrm{M}})$-close to (5.4). Analogously, we
consider the Poincaré map associated with this system,

$\displaystyle\Pi_{i_{\mathrm{M}}}:\Sigma\to\Sigma,\quad\Sigma=\\{h=0\\},$
(6.3)

which is a $\mathcal{O}(i_{\mathrm{M}})$-perturbation of the map $\Pi_{0}$
given in (5.1).

We study the system given by Hamiltonian $\mathcal{K}$ (see (6.1)) as a
perturbation of the one given by $\mathcal{K}_{\mathrm{CP}}$ (see (4.2)). By
classical Fenichel Theory (see [23]) and Corollary 5.2, one has that the flow
associated to the Hamiltonian $\mathcal{K}$ has a normally hyperbolic
invariant $3$-dimensional cylinder which is $i_{\mathrm{M}}$-close to the
cylinder $\Lambda_{0}$. Analogously, by Corollary 5.3, the Poincaré map
$\Pi_{i_{\mathrm{M}}}$ has a normally hyperbolic invariant $2$-dimensional
cylinder at section $\Sigma$.

Recall that, also by Corollary 5.3, the cylinder $\widetilde{\Lambda}_{0}$
possesses two homoclinic channels which we denote by $\mathrm{pri}$ and
$\mathrm{sec}$, which are defined for
$J\in\mathtt{D}^{*}_{0}\subset[J_{\min},J_{\max}]$ for
$*\in\\{\mathrm{pri},\mathrm{sec}\\}$. Moreover, the sets $\mathtt{D}^{*}_{0}$
satisfy that
$\mathtt{D}_{0}^{\mathrm{pri}}\cup\mathtt{D}_{0}^{\mathrm{sec}}=[J_{\min},J_{\max}]$.

For the full model, with $i_{\mathrm{M}}>0$ small enough, we still have
transverse homoclinic connections but in slightly smaller domains. In order to
characterize them, we define

$\mathtt{D}^{*}_{\delta}=\left\\{J\in\mathtt{D}^{*}_{0}\,:\,\mathrm{dist}(J,\partial\mathtt{D}^{*}_{0})>\delta\right\\},\quad\text{
for }*\in\\{\mathrm{pri},\mathrm{sec}\\}.$ (6.4)

Analogously, we also define

$\mathtt{D}_{\delta}=[J_{\min}+\delta,J_{\max}-\delta].$

Next theorem, which is a direct consequence of the classical Fenichel theory
and Corollary 5.3, summarizes the perturbative context.

######  Theorem 6.1.

Assume Ansatz 5.1. For any $\delta>0$ and $r\geq 4$, there exists
$i_{\mathrm{M}}^{0}>0$ such that, for any
$i_{\mathrm{M}}\in(0,i_{\mathrm{M}}^{0})$, the map $\Pi_{i_{\mathrm{M}}}$
introduced in (6.3) has the following properties.

  1. 1.

It has a $\mathcal{C}^{r}$ normally hyperbolic invariant manifold
$\widetilde{\Lambda}_{i_{\mathrm{M}},\delta}$, which is $i_{\mathrm{M}}$-close
in the $\mathcal{C}^{1}$-topology to $\widetilde{\Lambda}_{0}$. In addition,
there exists a function
$\mathcal{G}_{i_{\mathrm{M}}}:\mathtt{D}_{\delta}\times\mathbb{T}\to\mathbb{R}^{3}\times\mathbb{T}^{3}$
of the form

$\mathcal{G}_{i_{\mathrm{M}}}(J,\Omega_{\mathrm{M}})=\big{(}\mathcal{G}_{i_{\mathrm{M}}}^{\eta}(J,\Omega_{\mathrm{M}}),{\mathcal{G}}_{i_{\mathrm{M}}}^{\Gamma}(J,\Omega_{\mathrm{M}}),J,\mathcal{G}_{i_{\mathrm{M}}}^{\xi}(J,\Omega_{\mathrm{M}}),0,\Omega_{\mathrm{M}}\big{)},$

which parameterizes the cylinder $\widetilde{\Lambda}_{i_{\mathrm{M}},\delta}$
as a graph, that is

$\widetilde{\Lambda}_{i_{\mathrm{M}},\delta}=\\{\mathcal{G}_{i_{\mathrm{M}}}(J,\Omega_{\mathrm{M}})\,:\,(J,\Omega_{\mathrm{M}})\in\mathtt{D}_{\delta}\times\mathbb{T}\\}.$

  2. 2.

There exists a nonvanishing $\mathcal{C}^{r-1}$ function,
$a_{i_{\mathrm{M}}}$, such that the pull back of the canonical form
$d\xi\wedge d\eta+dh\wedge d\Gamma+d\Omega_{\mathrm{M}}\wedge dJ$ onto the
cylinder $\widetilde{\Lambda}_{i_{\mathrm{M}},\delta}$ is of the form

$a_{i_{\mathrm{M}}}(J,\Omega_{\mathrm{M}})d\Omega_{\mathrm{M}}\wedge dJ.$
(6.5)

  3. 3.

The homoclinic channels obtained in Corollary 5.3 are persistent. That is,
there exist $\mathcal{C}^{r}$ functions
$\mathcal{C}^{*}_{i_{\mathrm{M}}}:\mathtt{D}^{*}_{\delta}\times\mathbb{T}\to\mathbb{R}^{3}\times\mathbb{T}^{3},$

$\mathcal{C}^{*}_{i_{\mathrm{M}}}(J,\Omega_{\mathrm{M}})=\Big{(}\mathcal{C}^{\eta,*}_{i_{\mathrm{M}}}(J,\Omega_{\mathrm{M}}),\mathcal{C}_{i_{\mathrm{M}}}^{\Gamma,*}(J,\Omega_{\mathrm{M}}),J,\mathcal{C}^{\xi,*}_{i_{\mathrm{M}}}(J,\Omega_{\mathrm{M}}),0,\Omega_{\mathrm{M}}\Big{)},$

with $*\in\\{\mathrm{pri},\mathrm{sec}\\}$, such that the manifolds

$\Xi^{*}_{i_{\mathrm{M}},\delta}=\\{\mathcal{C}^{*}_{i_{\mathrm{M}}}(J,\Omega_{\mathrm{M}})\,:\,(J,\Omega_{\mathrm{M}})\in\mathtt{D}^{*}_{\delta}\times\mathbb{T}\\},$

belong to the transverse intersections between the invariant manifolds
$W^{\mathrm{u}}(\widetilde{\Lambda}_{i_{\mathrm{M}},\delta})$ and
$W^{\mathrm{s}}(\widetilde{\Lambda}_{i_{\mathrm{M}},\delta})$. Moreover, the
functions $\mathcal{C}^{*}_{i_{\mathrm{M}}}$ are $\mathcal{C}^{r}$ regular
with respect to $i_{\mathrm{M}}$ and are $\mathcal{O}(i_{\mathrm{M}})$ to the
function $\mathcal{C}^{*}_{0}$ obtained in Corollary 5.3.

We can define the inner and scattering maps in the invariant cylinder
$\widetilde{\Lambda}_{i_{\mathrm{M}},\delta}$ given in Theorem 6.1 as we have
done in Lemmas 5.4 and 5.5, respectively, for the coplanar case
($i_{\mathrm{M}}=0$). We also compute for them first order expansions. To this
end, we consider the following definition.

######  Definition 6.2.

Let $f$ be a $2\pi$-periodic function in $\Omega_{\mathrm{M}}$ and denote by
$f^{[k]}$ its $k$-th Fourier coefficient. Then, we define the set

$\mathcal{N}_{\Omega_{\mathrm{M}}}(f)=\\{k\in\mathbb{Z}\,:\,f^{[k]}\not\equiv
0\\}.$

Let us recall that
$\mathcal{H}=\mathcal{H}_{\mathrm{CP}}+\alpha^{3}i_{\mathrm{M}}\mathcal{R}$
(see (6.1)) and that $\mathcal{H}_{{}_{\mathrm{CP}}}$, the coplanar
Hamiltonian, is independent of $\Omega_{\mathrm{M}}$. However, the remainder
$\mathcal{R}$ has the following harmonic structure

$\displaystyle\mathcal{R}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}};i_{\mathrm{M}})=$
$\displaystyle\,\mathcal{R}_{1}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}})+i_{\mathrm{M}}\mathcal{R}_{2}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}})+\mathcal{O}(i_{\mathrm{M}}^{2}),$

with

$\mathcal{N}_{\Omega_{\mathrm{M}}}(\mathcal{R}_{1})=\\{\pm
1\\}\qquad\text{and}\qquad\mathcal{N}_{\Omega_{\mathrm{M}}}(\mathcal{R}_{2})=\\{0,\pm
2\\}.$

This is a consequence of the particular formulas of $\mathcal{R}$ in Appendix
A.2 (see equations (A.4) and (A.5)). To fix notation, we write
$\mathcal{R}_{1}$ as

$\mathcal{R}_{1}(\eta,\Gamma,\xi,h,\Omega_{\mathrm{M}})=e^{i\Omega_{\mathrm{M}}}\mathcal{R}_{1}^{+}(\eta,\Gamma,\xi,h)+e^{-i\Omega_{\mathrm{M}}}\mathcal{R}_{1}^{-}(\eta,\Gamma,\xi,h).$

Note that $\mathcal{R}_{1}^{-}=\overline{\mathcal{R}_{1}^{+}}$. We analyze the
harmonic structure of all the objects involved in the definition of the inner
and scattering maps. The first step is to study the asymptotic expansion with
respect to $i_{\mathrm{M}}$ of the flow associated to the vector field (6.2).

######  Lemma 6.3.

Let $\psi^{i_{\mathrm{M}}}\\{s,(\eta,\Gamma,J,\xi,h,\Omega_{\mathrm{M}})\\}$
be the flow associated with the vector field in (6.2) induced by the
Hamiltonian $\mathcal{K}$ in (6.1). Then, for $i_{\mathrm{M}}>0$ small enough,
it has an expansion

$\displaystyle\psi^{i_{\mathrm{M}}}\\{s,(\eta,\Gamma,J,\xi,h,\Omega_{\mathrm{M}})\\}=$
$\displaystyle\,\widehat{\varphi}_{\mathrm{CP}}\\{s,(\eta,\Gamma,J,\xi,h,\Omega_{\mathrm{M}})\\}+i_{\mathrm{M}}\psi_{1}\\{s,(\eta,\Gamma,J,\xi,h,\Omega_{\mathrm{M}})\\}$
$\displaystyle+i_{\mathrm{M}}^{2}\psi_{2}\\{s,(\eta,\Gamma,J,\xi,h,\Omega_{\mathrm{M}})\\}+\mathcal{O}(i_{\mathrm{M}}^{3}),$

where $\widehat{\varphi}_{\mathrm{CP}}$ has been introduced in (5.5) and the
remaining functions satisfy that

$\displaystyle\mathcal{N}_{\Omega_{\mathrm{M}}}(\psi_{1}\\{s,(\eta,\Gamma,J,\xi,h,\Omega_{\mathrm{M}})\\})$
$\displaystyle=\\{\pm 1\\},$
$\displaystyle\mathcal{N}_{\Omega_{\mathrm{M}}}(\psi_{2}\\{s,(\eta,\Gamma,J,\xi,h,\Omega_{\mathrm{M}})\\})$
$\displaystyle=\\{0,\pm 2\\}.$

The proof of this lemma is analogous to that of Lemma 3.6 in [22], and relies
on classical perturbation theory jointly with the special form of (6.2).

######  Lemma 6.4.

Assume Ansatz 5.1. Let $\mathcal{G}_{i_{\mathrm{M}}}$, $a_{i_{\mathrm{M}}}$
and $\mathcal{C}_{i_{\mathrm{M}}}^{*}$, for
$*\in\\{\mathrm{pri},\mathrm{sec}\\}$, be as given in Theorem 6.1. Then, for
$i_{\mathrm{M}}>0$ small enough, they have an asymptotic expansion of the form

$\displaystyle\mathcal{G}_{i_{\mathrm{M}}}(J,\Omega_{\mathrm{M}})$
$\displaystyle=\mathcal{G}_{0}(J,\Omega_{\mathrm{M}})+i_{\mathrm{M}}\mathcal{G}_{1}(J,\Omega_{\mathrm{M}})+i_{\mathrm{M}}^{2}\mathcal{G}_{2}(J,\Omega_{\mathrm{M}})+\mathcal{O}(i_{\mathrm{M}}^{3}),$
$\displaystyle a_{i_{\mathrm{M}}}(J,\Omega_{\mathrm{M}})$
$\displaystyle=1+i_{\mathrm{M}}a_{1}(J,\Omega_{\mathrm{M}})+i_{\mathrm{M}}^{2}a_{2}(J,\Omega_{\mathrm{M}})+\mathcal{O}(i_{\mathrm{M}}^{3}),$
$\displaystyle\mathcal{C}_{i_{\mathrm{M}}}^{*}(J,\Omega_{\mathrm{M}})$
$\displaystyle=\mathcal{C}_{0}^{*}(J,\Omega_{\mathrm{M}})+i_{\mathrm{M}}\mathcal{C}_{1}^{*}(J,\Omega_{\mathrm{M}})+i_{\mathrm{M}}^{2}\mathcal{C}_{2}^{*}(J,\Omega_{\mathrm{M}})+\mathcal{O}(i_{\mathrm{M}}^{3}),$

where $\mathcal{G}_{0}$ and $\mathcal{C}_{0}^{*}$ have been introduced in
Corollary 5.3 and the remaining functions satisfy that

$\displaystyle\mathcal{N}_{\Omega_{\mathrm{M}}}(\mathcal{G}_{1})=\\{\pm
1\\},\qquad\mathcal{N}_{\Omega_{\mathrm{M}}}(\mathcal{G}_{2})=\\{0,\pm 2\\},$
$\displaystyle\mathcal{N}_{\Omega_{\mathrm{M}}}(a_{1})=\\{\pm
1\\},\qquad\mathcal{N}_{\Omega_{\mathrm{M}}}(a_{2})=\\{0,\pm 2\\},$
$\displaystyle\mathcal{N}_{\Omega_{\mathrm{M}}}(\mathcal{C}_{1}^{*})=\\{\pm
1\\},\qquad\mathcal{N}_{\Omega_{\mathrm{M}}}(\mathcal{C}_{2}^{*})=\\{0,\pm
2\\},$

and that, for $i=1,2$,

$\displaystyle\mathcal{G}_{i}(J,\Omega_{\mathrm{M}})$
$\displaystyle=\big{(}\mathcal{G}_{i}^{\eta}(J,\Omega_{\mathrm{M}}),\mathcal{G}_{i}^{\Gamma}(J,\Omega_{\mathrm{M}}),0,\mathcal{G}_{i}^{\xi}(J,\Omega_{\mathrm{M}}),0,0\big{)},$
$\displaystyle\mathcal{C}_{i}^{*}(J,\Omega_{\mathrm{M}})$
$\displaystyle=\big{(}\mathcal{C}_{i}^{\eta,*}(J,\Omega_{\mathrm{M}}),\mathcal{C}_{i}^{\Gamma,*}(J,\Omega_{\mathrm{M}}),0,\mathcal{C}_{i}^{\xi,*}(J,\Omega_{\mathrm{M}}),0,0\big{)}.$

Lemmas 6.3 and 6.4 allow to compute asymptotic expansions in
$i_{\mathrm{M}}>0$ for the inner and scattering maps of the inclined model.
First, we introduce notation (see (5.6)) for the evolution of the coplanar
flow for points on the normally hyperbolic invariant cylinder
$\widetilde{\Lambda}_{0}$ and on the homoclinic channels $\Xi^{*}_{0}$ (see
Corollary 5.3):

$\displaystyle\gamma_{J}(s)$
$\displaystyle=\varphi_{\mathrm{CP}}\\{s,(0,\mathcal{G}_{0}^{\Gamma}(J),0,0)\\},$
$\displaystyle\chi^{\mathrm{pri}}_{J}(s)$
$\displaystyle=\varphi_{\mathrm{CP}}\\{s,(0,\mathcal{C}_{0}^{\Gamma,\mathrm{pri}}(J),\mathcal{C}_{0}^{\xi,\mathrm{pri}}(J),0)\\},$
(6.6) $\displaystyle\chi^{\mathrm{sec}}_{J}(s)$
$\displaystyle=\varphi_{\mathrm{CP}}\\{s,(\mathcal{C}_{0}^{\eta,\mathrm{sec}}(J),\mathcal{C}_{0}^{\Gamma,\mathrm{sec}}(J),0,0)\\},$
$\displaystyle\widetilde{\gamma}_{J}(s)$
$\displaystyle=\varphi^{\Omega_{\mathrm{M}}}\\{s,(0,\mathcal{G}_{0}^{\Gamma}(J),0,0)\\},$
$\displaystyle\widetilde{\chi}^{\mathrm{pri}}_{J}(s)$
$\displaystyle=\varphi^{\Omega_{\mathrm{M}}}\\{s,(0,\mathcal{C}_{0}^{\Gamma,\mathrm{pri}}(J),\mathcal{C}_{0}^{\xi,\mathrm{pri}}(J),0)\\},$
$\displaystyle\widetilde{\chi}^{\mathrm{sec}}_{J}(s)$
$\displaystyle=\varphi^{\Omega_{\mathrm{M}}}\\{s,(\mathcal{C}_{0}^{\eta,\mathrm{sec}}(J),\mathcal{C}_{0}^{\Gamma,\mathrm{sec}}(J),0,0)\\},$

where $\varphi_{\mathrm{CP}}$ and $\varphi^{\Omega_{\mathrm{M}}}$ have been
introduced in (5.5).

######  Theorem 6.5.

Assume Ansatz 5.1. Fix $\delta>0$ and $i_{\mathrm{M}}>0$ small enough. The
normally hyperbolic manifold $\widetilde{\Lambda}_{i_{\mathrm{M}},\delta}$
given in Theorem 6.1 of the map $\Pi_{i_{\mathrm{M}}}$ (see (6.3)) has
associated inner and outer maps which are $\mathcal{C}^{r}$ (also with respect
to $i_{\mathrm{M}}$) and are of the following form.

  * •

The inner map is of the form

$\mathcal{F}^{\mathrm{in}}_{i_{\mathrm{M}}}\begin{pmatrix}J\\\
\Omega_{M}\end{pmatrix}=\begin{pmatrix}J+i_{\mathrm{M}}A_{1}(J,\Omega_{\mathrm{M}})+i_{\mathrm{M}}^{2}A_{2}(J,\Omega_{\mathrm{M}})+\mathcal{O}(i_{\mathrm{M}}^{3})\\\
\Omega_{M}+n_{\Omega_{\mathrm{M}}}\left\\{\mathcal{T}_{0}(J)+i_{\mathrm{M}}\mathcal{T}_{1}(J,\Omega_{\mathrm{M}})+i_{\mathrm{M}}^{2}\mathcal{T}_{2}(J,\Omega_{\mathrm{M}})+\mathcal{O}(i_{\mathrm{M}}^{3})\right\\}\end{pmatrix},$

where $(J,\Omega_{\mathrm{M}})\in\mathtt{D}_{\delta}\times\mathbb{T}$, the
function $\mathcal{T}_{0}$ has been introduced in (5.2) and the functions
$A_{1},A_{2},\mathcal{T}_{1}$ and $\mathcal{T}_{2}$ satisfy that

$\displaystyle\mathcal{N}_{\Omega_{\mathrm{M}}}(A_{1}),\mathcal{N}_{\Omega_{\mathrm{M}}}(\mathcal{T}_{1})$
$\displaystyle=\\{\pm 1\\},$
$\displaystyle\mathcal{N}_{\Omega_{\mathrm{M}}}(A_{2}),\mathcal{N}_{\Omega_{\mathrm{M}}}(\mathcal{T}_{2})$
$\displaystyle=\\{0,\pm 2\\}.$

Moreover, $A_{1}$ is of the form

$\displaystyle A_{1}(J,\Omega_{\mathrm{M}})$
$\displaystyle=A_{1}^{+}(J)e^{i\Omega_{\mathrm{M}}}+A_{1}^{-}(J)e^{-i\Omega_{\mathrm{M}}},$

where

$\displaystyle A_{1}^{\pm}(J)$ $\displaystyle=\mp
i\alpha^{3}\int_{0}^{2\pi}\frac{\mathcal{R}_{1}^{\pm}(\gamma_{J}(s))}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}(\gamma_{J}(s))}e^{\pm
i\widetilde{\gamma}_{J}(s)}ds.$

  * •

For $*\in\\{\mathrm{pri},\mathrm{sec}\\}$, the scattering map is of the form

$\mathcal{F}^{\mathrm{out},*}_{i_{\mathrm{M}}}\begin{pmatrix}J\\\
\Omega_{M}\end{pmatrix}=\begin{pmatrix}J+i_{\mathrm{M}}B_{1}^{*}(J,\Omega_{\mathrm{M}})+i_{\mathrm{M}}^{2}B_{2}^{*}(J,\Omega_{\mathrm{M}})+\mathcal{O}(i_{\mathrm{M}}^{3})\\\
\Omega_{M}+n_{\Omega_{\mathrm{M}}}\left\\{\zeta^{*}(J)+i_{\mathrm{M}}D_{1}^{*}(J,\Omega_{\mathrm{M}})+i_{\mathrm{M}}^{2}D_{2}^{*}(J,\Omega_{\mathrm{M}})+\mathcal{O}(i_{\mathrm{M}}^{3})\right\\}\end{pmatrix},$

where $(J,\Omega_{\mathrm{M}})\in\mathtt{D}^{*}_{\delta}\times\mathbb{T}$, the
function $\zeta^{*}(J)=\zeta^{*}_{+}(J)-\zeta^{*}_{-}(J)$ has been introduced
in Lemma 5.5 and the functions $B^{*}_{1},B^{*}_{2},D^{*}_{1}$ and $D^{*}_{2}$
satisfy that

$\displaystyle\mathcal{N}_{\Omega_{\mathrm{M}}}(B_{1}^{*}),\mathcal{N}_{\Omega_{\mathrm{M}}}(D_{1}^{*})$
$\displaystyle=\\{\pm 1\\},$
$\displaystyle\mathcal{N}_{\Omega_{\mathrm{M}}}(B_{2}^{*}),\mathcal{N}_{\Omega_{\mathrm{M}}}(D_{2}^{*})$
$\displaystyle=\\{0,\pm 2\\}.$

Moreover, the functions $B_{1}^{*}$ are of the form

$\displaystyle
B_{1}^{*}(J,\Omega_{\mathrm{M}})=B_{1}^{*,+}(J)e^{i\Omega_{\mathrm{M}}}+B_{1}^{*,-}(J)e^{-i\Omega_{\mathrm{M}}},$

where

$\displaystyle B_{1}^{*,\pm}(J)=\pm i\alpha^{3}\lim_{s\to-\infty}\Bigg{[}$
$\displaystyle\int_{0}^{s}\frac{\mathcal{R}_{1}^{\pm}(\chi_{J}^{*}(\sigma))}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}(\chi^{*}_{J}(\sigma))}e^{\pm
i(\widetilde{\chi}^{*}_{J}(\sigma)+n_{\Omega_{\mathrm{M}}}\zeta_{+}^{*}(J))}d\sigma$
$\displaystyle-\int_{0}^{s}\frac{\mathcal{R}_{1}^{\pm}(\gamma_{J}(\sigma))}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}(\gamma_{J}(\sigma))}e^{\pm
i\widetilde{\gamma}_{J}(\sigma)}d\sigma\Bigg{]}$ $\displaystyle\mp
i\alpha^{3}\lim_{s\to+\infty}\Bigg{[}$
$\displaystyle\int_{0}^{s}\frac{\mathcal{R}_{1}^{\pm}(\chi_{J}^{*}(\sigma))}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}(\chi^{*}_{J}(\sigma))}e^{\pm
i(\widetilde{\chi}^{*}_{J}(\sigma)+n_{\Omega_{\mathrm{M}}}\zeta_{+}^{*}(J))}d\sigma$
$\displaystyle-\int_{0}^{s}\frac{\mathcal{R}_{1}^{\pm}(\gamma_{J}(\sigma))}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}(\gamma_{J}(\sigma))}e^{\pm
i(\widetilde{\gamma}_{J}(\sigma)+2n_{\Omega_{\mathrm{M}}}\zeta_{+}^{*}(J))}d\sigma\Bigg{]}.$

###### Proof.

First, we focus on the inner map. Its regularity is a consequence of the
regularity of the original flow and the regularity of the invariant cylinder
(which is a consequence of Fenichel Theory). Notice that, since it is the
dynamics of the Poincaré map $\Pi_{i_{\mathrm{M}}}$ (defined in (6.3))
restricted to $\widetilde{\Lambda}_{i_{\mathrm{M}},\delta}$, it satisfies the
following homological equation

$\displaystyle{\Pi}_{i_{\mathrm{M}}}\circ\mathcal{G}_{i_{\mathrm{M}}}=\mathcal{G}_{i_{\mathrm{M}}}\circ\mathcal{F}^{\mathrm{in}}_{i_{\mathrm{M}}},$
(6.7)

where $\mathcal{G}_{i_{\mathrm{M}}}$ has been introduced in Theorem 6.1.
Notice that $\Pi_{i_{\mathrm{M}}}$ is defined by

$\displaystyle\Pi_{i_{\mathrm{M}}}(\eta,\Gamma,J,\xi,0,\Omega_{\mathrm{M}})=\psi^{i_{\mathrm{M}}}\\{2\pi,(\eta,\Gamma,J,\xi,0,\Omega_{\mathrm{M}})\\},$

where $\psi^{i_{\mathrm{M}}}$ is the flow associated to the vector field
(6.2). Therefore, by Lemma 6.3, one can consider the following expansion

$\displaystyle\Pi_{i_{\mathrm{M}}}={\Pi}_{0}+i_{\mathrm{M}}\Pi_{1}+i_{\mathrm{M}}^{2}\Pi_{2}+\mathcal{O}(i_{\mathrm{M}}^{3}),$

where,

$\displaystyle\Pi_{0}(\eta,\Gamma,J,\xi,0,\Omega_{\mathrm{M}})$
$\displaystyle=\widehat{\varphi}_{\mathrm{CP}}\\{2\pi,(\eta,\Gamma,J,\xi,0,\Omega_{\mathrm{M}})\\},$
$\displaystyle\Pi_{k}(\eta,\Gamma,J,\xi,0,\Omega_{\mathrm{M}})$
$\displaystyle=\psi_{k}\\{2\pi,(\eta,\Gamma,J,\xi,0,\Omega_{\mathrm{M}})\\},\quad
k=1,2.$

Let us consider the following expansions of the inner map

$\displaystyle\mathcal{F}^{\mathrm{in}}_{i_{\mathrm{M}}}$
$\displaystyle={\mathcal{F}}_{0}^{\mathrm{in}}+i_{\mathrm{M}}\mathcal{F}^{\mathrm{in}}_{1}+i_{\mathrm{M}}^{2}\mathcal{F}^{\mathrm{in}}_{2}+\mathcal{O}(i_{\mathrm{M}}^{3}).$

By (6.7) and Lemma 6.4, the map $\mathcal{F}_{0}^{\mathrm{in}}$ is as given in
(5.2). Moreover,

$\Pi_{1}\circ\mathcal{G}_{0}+(D\Pi_{0}\circ\mathcal{G}_{0})\mathcal{G}_{1}=\mathcal{G}_{1}\circ\mathcal{F}_{0}^{\mathrm{in}}+(D\mathcal{G}_{0}\circ\mathcal{F}_{0}^{\mathrm{in}})\mathcal{F}_{1}^{\mathrm{in}}.$

and

$\displaystyle\Pi_{2}\circ\mathcal{G}_{0}+(D\Pi_{1}\circ\mathcal{G}_{0})\mathcal{G}_{1}+\frac{1}{2}(D^{2}\Pi_{0}\circ\mathcal{G}_{0})\mathcal{G}_{1}^{\otimes
2}+(D\Pi_{0}\circ\mathcal{G}_{0})\mathcal{G}_{2}$ (6.8)
$\displaystyle=\mathcal{G}_{2}\circ\mathcal{F}_{0}^{\mathrm{in}}+(D\mathcal{G}_{1}\circ\mathcal{F}_{0}^{\mathrm{in}})\mathcal{F}_{1}^{\mathrm{in}}+\frac{1}{2}(D^{2}\mathcal{G}_{0}\circ\mathcal{F}_{0}^{\mathrm{in}})+(D\mathcal{G}_{0}\circ\mathcal{F}_{0}^{\mathrm{in}})\mathcal{F}_{2}^{\mathrm{in}}.$

Then, by Corollary 5.3 and Lemma 6.4,

$\pi_{\diamond}\\{\mathcal{F}_{1}^{\mathrm{in}}\\}=\pi_{\diamond}\\{\Pi_{1}\circ\mathcal{G}_{0}+(D\Pi_{0}\circ\mathcal{G}_{0})\mathcal{G}_{1}\\}\\},\quad\diamond=J,\Omega_{\mathrm{M}},$
(6.9)

where $\pi_{J}$ and $\pi_{\Omega_{\mathrm{M}}}$ denote the projections onto
the coordinates $J$ and $\Omega_{\mathrm{M}}$, respectively. Since
$\mathcal{N}_{\Omega_{\mathrm{M}}}(\mathcal{G}_{1}),\mathcal{N}_{\Omega_{\mathrm{M}}}(\Pi_{1})=\\{\pm
1\\}$ one has that
$\mathcal{N}_{\Omega_{\mathrm{M}}}(\mathcal{F}_{1}^{\mathrm{in}})=\\{\pm
1\\}.$ Moreover, using analogous arguments and equation (6.8), one can obtain
that
$\mathcal{N}_{\Omega_{\mathrm{M}}}(\mathcal{F}_{2}^{\mathrm{in}})=\\{0,\pm
2\\}.$

Now it only remains to compute the formula for $A_{1}(J,\Omega_{\mathrm{M}})$.
Indeed, by (6.9) and taking into account that
$\partial_{J}\Pi_{0}=(0,0,1,0,0,0)$, we have that

$\displaystyle
A_{1}(J,\Omega_{\mathrm{M}})=\pi_{J}\\{\mathcal{F}_{1}^{\mathrm{in}}(J,\Omega_{\mathrm{M}})\\}=\pi_{J}\\{\Pi_{1}\circ\mathcal{G}_{0}\\}=\pi_{J}\,{\psi_{1}\\{2\pi,\mathcal{G}_{0}(J,\Omega_{\mathrm{M}})\\}}.$

Since $\psi^{i_{\mathrm{M}}}$ is the flow associated to the vector field
(6.2), we can apply the fundamental theorem of calculus and the expression of
$\mathcal{R}$ in (A.2) to obtain

$\displaystyle\pi_{J}\,\psi^{i_{\mathrm{M}}}\\{2\pi,\mathcal{G}_{0}(J,\Omega_{\mathrm{M}})\\}$
$\displaystyle=-\alpha^{3}i_{\mathrm{M}}\int_{0}^{2\pi}\frac{\partial_{\Omega_{\mathrm{M}}}\mathcal{R}(\psi^{i_{\mathrm{M}}}\\{s,\mathcal{G}_{0}(J,\Omega_{\mathrm{M}})\\};i_{\mathrm{M}})}{\partial_{\Gamma}\mathcal{H}(\psi^{i_{\mathrm{M}}}\\{s,\mathcal{G}_{0}(J,\Omega_{\mathrm{M}})\\};i_{\mathrm{M}})}ds$
$\displaystyle=-\alpha^{3}i_{\mathrm{M}}\int_{0}^{2\pi}\frac{\partial_{\Omega_{\mathrm{M}}}\mathcal{R}_{1}(\widehat{\varphi}_{\mathrm{CP}}\\{s,\mathcal{G}_{0}(J,\Omega_{\mathrm{M}})\\})}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}(\varphi_{\mathrm{CP}}\\{s,\mathcal{G}_{0}(J,\Omega_{\mathrm{M}})\\})}ds+\mathcal{O}(i_{\mathrm{M}}^{2}).$

Then, by the expression of $\mathcal{R}_{1}$ in (A.4) and of functions (6.6),
one has that

$\displaystyle A_{1}(J,\Omega_{\mathrm{M}})=$
$\displaystyle\,-i\alpha^{3}\int_{0}^{2\pi}\frac{\mathcal{R}_{1}^{+}(\gamma_{J}(s))e^{i(\Omega_{\mathrm{M}}+\widetilde{\gamma}_{J}(s))}}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}(\gamma_{J}(s))}ds$
$\displaystyle+i\alpha^{3}\int_{0}^{2\pi}\frac{\mathcal{R}_{1}^{-}(\gamma_{J}(s))e^{-i(\Omega_{\mathrm{M}}+\widetilde{\gamma}_{J}(s))}}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}(\gamma_{J}(s))}ds.$

Next, we compute the scattering maps for a fixed
$*=\\{\mathrm{pri},\mathrm{sec}\\}$ (see Definition 4.4). Its regularity is
proven in [17] (for regular vector fields with regular normally hyperbolic
invariant manifolds).

Let us consider points
$\mathcal{G}_{i_{\mathrm{M}}}(J^{+},\Omega_{{\mathrm{M}}}^{+}),\mathcal{G}_{i_{\mathrm{M}}}(J^{-},\Omega_{{\mathrm{M}}}^{-})\in\widetilde{\Lambda}_{i_{\mathrm{M}},\delta}$
and
$\mathcal{C}_{i_{\mathrm{M}}}^{*}(J^{0},\Omega_{{\mathrm{M}}}^{0})\in\Xi_{i_{\mathrm{M}},\delta}^{*}$
(see Theorem 6.1) such that the unstable fiber of
$\mathcal{G}_{i_{\mathrm{M}}}(J^{-},\Omega_{{\mathrm{M}}}^{-})$ intersects the
stable fiber of
$\mathcal{G}_{i_{\mathrm{M}}}(J^{+},\Omega_{{\mathrm{M}}}^{+})$ at the point
$\mathcal{C}_{i_{\mathrm{M}}}^{*}(J^{0},\Omega_{{\mathrm{M}}}^{0})$ in the
homoclinic channel $\Xi_{i_{\mathrm{M}},\delta}^{*}$. This implies that

$\lim_{s\to\pm\infty}\big{|}\psi^{i_{\mathrm{M}}}\\{s,\mathcal{G}_{i_{\mathrm{M}}}(J^{\pm},\Omega_{{\mathrm{M}}}^{\pm})\\}-\psi^{i_{\mathrm{M}}}\\{s,\mathcal{C}_{i_{\mathrm{M}}}^{*}(J^{0},\Omega_{{\mathrm{M}}}^{0})\\}\big{|}=0,$
(6.10)

where $\psi^{i_{\mathrm{M}}}$ is the flow associated to the vector field
(6.2). In other words, the scattering map satisfies that
$(J^{+},\Omega_{{\mathrm{M}}}^{+})=\mathcal{F}^{\mathrm{out},*}_{i_{\mathrm{M}}}(J^{-},\Omega_{{\mathrm{M}}}^{-})$.

The analysis of the harmonic structure can be done as for the inner map. We
show now how to compute the first order of the $J$-component. Recall that
$\pi_{J}\mathcal{G}_{i_{\mathrm{M}}}(J,\Omega_{\mathrm{M}})=\pi_{J}\mathcal{C}_{i_{\mathrm{M}}}^{*}(J,\Omega_{\mathrm{M}})=J$.
By the fundamental theorem of calculus and using (6.10), one has that

$\displaystyle J^{0}-J^{\pm}=\lim_{s\to{\pm}\infty}\Bigg{[}$
$\displaystyle\int_{s}^{0}\pi_{J}\partial_{s}\psi^{i_{\mathrm{M}}}\\{\sigma,\mathcal{C}_{i_{\mathrm{M}}}^{*}(J^{0},\Omega_{{\mathrm{M}}}^{0})\\}d\sigma$
$\displaystyle-\int_{s}^{0}\pi_{J}\partial_{s}\psi^{i_{\mathrm{M}}}\\{\sigma,\mathcal{G}_{i_{\mathrm{M}}}(J^{\pm},\Omega_{{\mathrm{M}}}^{\pm})\\}d\sigma\Bigg{]}.$

Then, by (6.2),

$\displaystyle
J^{0}-J^{\pm}=\alpha^{3}i_{\mathrm{M}}\lim_{s\to\pm\infty}\Bigg{(}$
$\displaystyle\int_{0}^{s}\frac{\partial_{\Omega_{\mathrm{M}}}\mathcal{R}(\psi^{i_{\mathrm{M}}}\\{\sigma,\mathcal{C}_{i_{\mathrm{M}}}^{*}(J^{0},\Omega_{{\mathrm{M}}}^{0}\\};i_{\mathrm{M}})}{\partial_{\Gamma}\mathcal{H}(\psi^{i_{\mathrm{M}}}\\{\sigma,\mathcal{C}_{i_{\mathrm{M}}}^{*}(J^{0},\Omega_{{\mathrm{M}}}^{0})\\};i_{\mathrm{M}})}d\sigma$
$\displaystyle-\int_{0}^{s}\frac{\partial_{\Omega_{\mathrm{M}}}\mathcal{R}(\psi^{i_{\mathrm{M}}}\\{\sigma,\mathcal{G}_{i_{\mathrm{M}}}(J^{\pm},\Omega_{{\mathrm{M}}}^{\pm})\\};i_{\mathrm{M}})}{\partial_{\Gamma}\mathcal{H}(\psi^{i_{\mathrm{M}}}\\{\sigma,\mathcal{G}_{i_{\mathrm{M}}}(J^{\pm},\Omega_{{\mathrm{M}}}^{\pm})\\};i_{\mathrm{M}})}d\sigma\Bigg{)}$

and, by Lemmas 6.3, 6.4 and the expression of $\mathcal{R}$ in (A.2), we
obtain

$\displaystyle
J^{0}-J^{\pm}=\alpha^{3}i_{\mathrm{M}}\lim_{s\to\pm\infty}\Bigg{(}$
$\displaystyle\int_{0}^{s}\frac{\partial_{\Omega_{\mathrm{M}}}\mathcal{R}_{1}(\widehat{\varphi}_{\mathrm{CP}}\\{\sigma,\mathcal{C}_{0}^{*}(J^{0},\Omega_{{\mathrm{M}}}^{0})\\})}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}(\widehat{\varphi}_{\mathrm{CP}}\\{\sigma,\mathcal{C}_{0}^{*}(J^{0},\Omega_{{\mathrm{M}}}^{0})\\})}d\sigma$
$\displaystyle-\int_{0}^{s}\frac{\partial_{\Omega_{\mathrm{M}}}\mathcal{R}_{1}(\widehat{\varphi}_{\mathrm{CP}}\\{\sigma,\mathcal{G}_{0}(J^{\pm},\Omega_{{\mathrm{M}}}^{\pm})\\})}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}(\widehat{\varphi}_{\mathrm{CP}}\\{\sigma,\mathcal{G}_{0}(J^{\pm},\Omega_{{\mathrm{M}}}^{\pm})\\})}d\sigma\Bigg{)}+\mathcal{O}(i_{\mathrm{M}}^{2}).$

Finally, by the expression of $\mathcal{R}_{1}$ in (A.4), of
$\widehat{\varphi}_{\mathrm{CP}}$ in (5.5) and of the functions (6.6)

$\displaystyle J^{0}$
$\displaystyle-J^{\pm}=i\alpha^{3}i_{\mathrm{M}}\lim_{s\to\pm\infty}$
$\displaystyle\Bigg{(}\int_{0}^{s}\frac{\mathcal{R}_{1}^{+}(\chi^{*}_{J^{0}}(\sigma))e^{i\Omega_{{\mathrm{M}}}^{0}+i\widetilde{\chi}_{J^{0}}^{*}(\sigma)}-\mathcal{R}_{1}^{-}(\chi^{*}_{J^{0}}(\sigma))e^{-i\Omega_{{\mathrm{M}}}^{0}-i\widetilde{\chi}_{J^{0}}^{*}(\sigma)}}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}(\chi^{*}_{J^{0}}(\sigma))}d\sigma$
$\displaystyle-\int_{0}^{s}\frac{\mathcal{R}_{1}^{+}(\gamma_{J^{\pm}}(\sigma))e^{i\Omega_{{\mathrm{M}}}^{\pm}+i\widetilde{\gamma}_{J^{\pm}}(\sigma)}-\mathcal{R}_{1}^{-}(\gamma_{J^{\pm}}(\sigma))e^{-i\Omega_{{\mathrm{M}}}^{\pm}-i\widetilde{\gamma}_{J^{\pm}}(\sigma)}}{\partial_{\Gamma}\mathcal{H}_{\mathrm{CP}}(\gamma_{J^{\pm}}(\sigma))}d\sigma\Bigg{)}+\mathcal{O}(i_{\mathrm{M}}^{2}).$

Notice that, by Lemma 5.5, (5.9) and (5.10), one has that

$\displaystyle\Omega_{{\mathrm{M}}}^{\pm}$
$\displaystyle=\Omega_{{\mathrm{M}}}^{0}+n_{\Omega_{\mathrm{M}}}\zeta^{*}_{\pm}(J^{0})+\mathcal{O}(i_{\mathrm{M}}),$

with $\zeta_{+}^{*}=-\zeta_{-}^{*}$. Therefore, since
$J^{\pm}=J^{0}+\mathcal{O}(i_{\mathrm{M}})$, one has that

$\displaystyle\Omega_{{\mathrm{M}}}^{0}$
$\displaystyle=\Omega_{{\mathrm{M}}}^{-}+n_{\Omega_{\mathrm{M}}}\zeta^{*}_{+}(J^{-})+\mathcal{O}(i_{\mathrm{M}}),$
$\displaystyle\Omega_{{\mathrm{M}}}^{+}$
$\displaystyle=\Omega_{{\mathrm{M}}}^{-}+2n_{\Omega_{\mathrm{M}}}\zeta^{*}_{+}(J^{-})+\mathcal{O}(i_{\mathrm{M}})$

and, as a consequence,

$\displaystyle J^{+}=$ $\displaystyle\,J^{-}+(J^{0}-J^{-})-(J^{0}-J^{+})$
$\displaystyle=$
$\displaystyle\,J^{-}+i_{\mathrm{M}}\left(B_{1}^{*,+}(J^{-})e^{i\Omega_{{\mathrm{M}}}^{-}}+B_{1}^{*,-}(J^{-})e^{-i\Omega_{{\mathrm{M}}}^{-}}\right)+\mathcal{O}(i_{\mathrm{M}}^{2}),$

with $B_{1}^{*,\pm}$ are the functions defined in the statement of the result.
∎

###  6.2 Existence of diffusing orbits

Once we have computed the first orders in $i_{\mathrm{M}}$ of both the inner
and the scattering maps (see Theorem 6.5), the next step is to construct a
“drifting pseudo-orbit”. By a pseudo-orbit, we mean a sequence of points in
the cylinder obtained by applying successively iterations of the inner map
$\mathcal{F}_{i_{\mathrm{M}}}^{\mathrm{in}}$ and the scattering maps
$\mathcal{F}_{i_{\mathrm{M}}}^{\mathrm{out},\mathrm{pri}}$ and
$\mathcal{F}_{i_{\mathrm{M}}}^{\mathrm{out},\mathrm{sec}}$. By drifting, we
mean that we look for a pseudo-orbit such that its initial condition is close
to the bottom of the cylinder, i.e $J\sim J_{\min}$, and that eventually it
hits a neighborhood of its top, i.e., $J\sim J_{\max}$.

To construct this pseudo-orbit we rely on the following ansatz, which we
verify numerically. See Figures 11, 12 and 13.

Figure 11: The components of $A_{1}^{+}$ as a function of
${\mathcal{H}}_{\mathrm{CP}}$ (non-dimensional units).

Figure 12:  For $*\in\\{\mathrm{pri},\mathrm{sec}\\}$, we plot the real
functions $B_{1,\cos}^{*}=B_{1}^{+}+B_{1}^{-}$ and
$B_{1,\sin}^{*}=i(B_{1}^{+}-B_{1}^{-})$ as a function of
${\mathcal{H}}_{\mathrm{CP}}$ (non-dimensional units). Left, the
$*=\mathrm{pri}$ case. Right, the $*=\mathrm{sec}$ case for the non-transverse
values (see Table 1).

######  Ansatz 6.6.

For any $\delta>0$ and $*\in\\{\mathrm{pri},\mathrm{sec}\\}$, the following
functions of $J$

$\displaystyle f^{*}_{\pm}(J)$ $\displaystyle=\left(e^{\pm
in_{\Omega_{\mathrm{M}}}\mathcal{T}_{0}(J)}-1\right)B_{1}^{*,\pm}(J)-\left(e^{\pm
in_{\Omega_{\mathrm{M}}}\zeta^{*}(J)}-1\right)A_{1}^{\pm}(J),$

do not vanish in the domain $\mathtt{D}_{\delta}^{*}$.

The following theorem ensures the existence of a pseudo-orbit. Its proof is
deferred to Section 6.3.

Figure 13: As a function of ${\mathcal{H}}_{\mathrm{CP}}\equiv-J$, we show the
norm of $f^{\mathrm{pri}}_{+}(J)$ in Ansatz 6.6 on the left. The minimum value
computed is 0.0014295 and corresponds to
${\mathcal{H}}_{\mathrm{CP}}=4.81143\times 10^{-7}$. On the right, the same
norm computed for the non-transverse cases (see Table 1). The corresponding
integrals have been computed by means of the function qags of the quadpack
Fortran package.

######  Theorem 6.7.

Assume Ansätze 5.1, 5.6, 6.6. Fix $\nu>0$ and $\delta>0$. Then, if
$i_{\mathrm{M}}>0$ is small enough, there exist
$z_{0}\in\widetilde{\Lambda}_{i_{\mathrm{M}},\delta}$, $N\in\mathbb{N}$ and
$l_{m}\in\\{\mathrm{pri},\mathrm{sec}\\}$, $j_{m},k_{m}\in\mathbb{N}$ for
$m=1\dots N$ such that the sequence

$z_{m}=\left(\mathcal{F}_{i_{\mathrm{M}}}^{\mathrm{in}}\right)^{j_{m}}\circ\mathcal{F}_{i_{\mathrm{M}}}^{\mathrm{out},l_{m}}\circ\left(\mathcal{F}_{i_{\mathrm{M}}}^{\mathrm{in}}\right)^{k_{m}}(z_{m-1})$

satisfies $z_{m}\in\widetilde{\Lambda}_{i_{\mathrm{M}},\delta}$ for all
$m=1\ldots N$ and the initial and final points $z_{0}$ and $z_{N}$ satisfy

$|\pi_{J}z_{0}-J_{\min}|\leq\nu\quad\text{and}\quad|\pi_{J}z_{N}-J_{\max}|\leq\nu,$

where $\pi_{J}$ denotes the projection onto the $J$ component.

Moreover, in the regions

$\mathcal{J}_{\min}=(J_{\min},J_{\min}+\nu)\times\mathbb{T},\qquad\mathcal{J}_{\max}=(J_{\max}-\nu,J_{\max})\times\mathbb{T},$

the inner map $\mathcal{F}_{i_{\mathrm{M}}}^{\mathrm{in}}$ has invariant tori
which are a graph over $h$ and whose dynamics are conjugated to rigid quasi-
periodic rotation.

Once we have a pseudo-orbit, the final step to prove Theorem 2.3 is to
construct a true orbit of the Poincaré map $\Pi_{i_{\mathrm{M}}}$ which
shadows the pseudo-orbit, that is, which visits small neighborhoods of the
points $z_{m}$ of the pseudo-orbit given by Theorem 6.7. This shadowing will
rely on the following theorem by M. Gidea, R. de la Llave and T. M. Seara in
[30].

######  Theorem 6.8 (Lemma 3.2 in [30]).

Assume that $f:\mathcal{M}\to\mathcal{M}$ is a $\mathcal{C}^{r}$ map with
$r\geq 4$, $\Lambda\subset\mathcal{M}$ is a normally hyperbolic invariant
manifold, $\Gamma_{j}$, $j=1\ldots N$ are homoclinic channels and
$\sigma^{j}:D_{1}^{j}\to D_{2}^{j}$ are the associated scattering maps, where
$D_{1,2}^{j}\subset\Lambda$ are open sets. Assume that $\Lambda$ and
$\Gamma_{j}$ are compact manifolds (possibly with boundary).

Then, for every $\nu>0$, there exist functions
$n_{i}^{*}:\mathbb{N}^{i}\to\mathbb{N}$ and functions
$m_{i}^{*}:\mathbb{N}^{2i+1}\times\mathbb{N}^{i+1}\to\mathbb{N}$ such that for
every pseudo-orbit $\\{y_{i}\\}_{y\geq 0}$ in $\Lambda$ of the form

$y_{i+1}=f^{m_{i}}\circ\sigma^{\alpha_{i}}\circ f^{n_{i}}(y_{i}),$

with $n_{i}\geq n_{i}^{*}(\alpha_{0},\ldots,\alpha_{i-1})$, $m_{i}\geq
m_{i}^{*}(n_{0},\ldots,n_{i},m_{0},\ldots.m_{i-1},\alpha_{0},\ldots,\alpha_{i})$
and $\alpha_{i}\in\\{1,\ldots,N\\}$, there exists an orbit $\\{z_{i}\\}_{i\geq
0}$ of $f$ in $\mathcal{M}$ such that, for all $i\geq 0$,

$z_{i+1}=f^{m_{i}+n_{i}}(z_{i})\quad\text{and}\quad d(z_{i},y_{i})<\nu.$

Theorem 2.3 is a direct consequence of the following lemma.

######  Lemma 6.9.

Assume Ansätze 5.1, 5.6, 6.6 and fix $\nu>0$ small. Then, for
$i_{\mathrm{M}}>0$ small enough, there exist a point $p_{0}\in\Sigma$ (see
(6.3)) and $k_{1}\ldots k_{N}$ such that the iterates

$p_{j}=\Pi_{i_{\mathrm{M}}}^{k_{j}}(p_{0})$

satisfy that

$\left|p_{j}-\mathcal{G}_{i_{\mathrm{M}}}(z_{j})\right|\leq 2\nu,\qquad
j=1\ldots N,$

where $\\{z_{j}\\}_{i=0}^{N}$ is the pseudo-orbit obtained in Theorem 6.7 and
$\mathcal{G}_{i_{\mathrm{M}}}$ is the parameterization of the cyilinder given
in Lemma 6.1.

Moreover,

$|\pi_{J}p_{0}-J_{\min}|\leq
2\nu\qquad\text{and}\qquad|\pi_{J}p_{N}-J_{\max}|\leq 2\nu.$

###### Proof.

To prove Lemma 6.9 we have to apply Theorem 6.8. To this end, we denote by
$\Lambda$ a domain of the cylinder
$\widetilde{\Lambda}_{i_{\mathrm{M}},\delta}$ delimited by one of the “‘top”
and “bottom” invariant tori given by Theorem 6.7. This makes $\Lambda$ compact
and invariant. Note that the projection onto $J$ of the homoclinic channels
$\Xi^{*}_{i_{\mathrm{M}},\delta}$ obtained in Theorem 6.1 covers the whole
interval $[J_{\min}+\delta,J_{\max}-\delta]$ and that one can choose compact
sets inside the homoclinic channels that still cover the same interval.

Then, to apply Theorem 6.8, it only remains to ensure that the amount of
iterates of the inner maps is large enough to fit the hypotheses of the
theorem. We follow ideas developed in [30].

Since $\Lambda$ is compact and invariant by $\Pi_{i_{\mathrm{M}}}$ and
$\mathcal{F}_{i_{\mathrm{M}}}^{\mathrm{in}}$ is an area preserving map, we can
use the Poincaré Recurrence Theorem to construct a pseudo-orbit which is
arbitrarily close to that of Theorem 6.7 and satisfies the hypotheses of
Theorem 6.8. Indeed, the Poincaré Recurrence Theorem assures that there exists
$\tilde{z}_{0}$ and $\tilde{k}_{1}$ big enough such that
$\widetilde{Q}_{0}:=\big{(}\mathcal{F}^{\mathrm{in}}_{i_{\mathrm{M}}}\big{)}^{\tilde{k}_{1}}(\tilde{z}_{0})$
is as close as necessary to