# Computer Assisted Proof of Drift Orbits Along Normally Hyperbolic Manifolds
II: Application to the Restricted Three Body Problem

Maciej J. Capiński 111Partially supported by the NCN grant
2018/29/B/ST1/00109. Faculty of Applied Mathematics, AGH University of Science
and Technology, al. Mickiewicza 30, 30-059 Kraków, Poland.  Natalia Wodka
222Partially supported by the NCN grant 2018/29/B/ST1/00109.

###### Abstract

We present a computer assisted proof or diffusion in the Planar Elliptic
Restricted Three Body Problem. We treat the elliptic problem as a perturbation
of the circular problem, where the perturbation parameter is the eccentricity
of the primaries. The unperturbed system preserves energy, and we show that
for sufficiently small perturbations we have orbits with explicit energy
changes, independent from the size of the perturbation. The result is based on
shadowing of orbits along transversal intersections of stable/unstable
manifolds of a normally hyperbolic cylinder.

###### keywords:

Normally hyperbolic manifold , Arnold diffusion , scattering map , topological
shadowing , computer assisted proof , three body problem

###### MSC:

[2010] 37J25 , 37J40

##  1 Introduction

Autonomous Hamiltonian systems have an integral of motion given by the
Hamiltonian. We shall refer to this integral as the energy. Our question is
whether for an arbitrarily small perturbation we can have orbits, which change
in energy by a prescribed constant independent form the size of the
perturbation.

In our case, the system of interest is the Planar Restricted Three Body
Problem that describes the motion of a small massless particle under the
gravitational pull of two large bodies. They rotate on Keplerian orbits and we
call them primaries. The motion of the massless particle is on the same plane
as the primaries. When the Keplerian orbit is circular, then in the coordinate
frame which rotates with the primaries we have an autonomous Hamiltonian
system; the Planar Circular Restricted Three Body Problem (PCR3BP). When the
Keplerian orbits are elliptic, then we deal with the Elliptic Problem
(PER3BP), which is a time dependent Hamiltonian system.

The system we study has a normally hyperbolic invariant manifold (NHIM) prior
to the perturbation, with stable and unstable manifolds which intersect
transversally. This means that we consider an ‘apriori-chaotic’ system. (This
is a much simpler problem than starting with a fully integrable system; as in
the famous example of Arnold [1].) Our diffusion mechanism is based on
establishing the existence of trajectories that shadow the intersections of
the stable and unstable manifolds and change energy under the influence of the
perturbation.

Our result is based on the geometric and shadowing results of Delshams, Gidea,
de la Llave and Seara. At the core of the mechanism is the scattering map [2];
that is a map which for a point on the NHIM assigns another point from the
NHIM, if their unstable and stable fibres intersect in a nontrivial way. The
benefit of studying a scattering map is twofold. Firstly, one can exploit
perturbative techniques to study the effect of the parameter on the scattering
map. This is typically done by using Melnikov integrals in the case of
continuous systems, or Melnikov sums in the case of discrete systems [2, 3].
Once the effects of the perturbation on the scattering map is established, the
second benefit is that there exist true orbits of the system, which will
shadow the ‘pseudo orbits’ of the scattering map [4]. This means that if one
can prove existence of ‘pseudo orbits’ of the scattering map which have
macroscopic changes in energy, then this ensures the existence of true
trajectories of the system that will do the same.

The strategy for proving diffusion using this scheme is described in [2, 3, 4,
5]. For our proof we will use the results from [6], which provides a
formulation for these methods, that can be implemented to obtain a computer
assisted proof. The main feature of [6] is that the assumptions needed for the
diffusion mechanism can be validated in a finite number of steps by checking
various bounds on the properties of the system.

In our case the components that need to be established are the following. We
establish bounds on the family of Lyapunov orbits in the PCR3BP, which
constitute our NHIM prior to the perturbation. Next we establish bounds for
the local stable/unstable manifolds of the NHIM, and prove that they intersect
transversally. We ensure that the scattering map is well defined and obtain
explicit bounds on its dynamics. We also validate twist conditions, which
ensure that the NHIM after perturbation contains a Cantor set of KAM tori.
Then we can apply the theorems from [6] to prove diffusion. All the above
steps can be established with the assistance of rigorous, interval arithmetic
based estimates, performed with the aid of a computer.

For proofs of diffusion in the PER3BP a computer assisted approach is not
strictly necessary. The paper [7] provides an analytic proof of diffusion in
the PER3BP. This result also used the scattering map theory and shadowing as
the mechanism. To obtain the analytic proof though [7] required that the mass
of one of the primaries is sufficiently small, and that the angular momentum
of the massless particle is sufficiently large. The difference with the result
from the current paper is that we can work with the explicit mass of Jupiter
and with the explicit energy of a massless particle, which corresponds to that
of the comet Oterma. As in [7] we have to assume that the eccentricity of the
system is sufficiently small.

A recent result which works with explicit masses, explicit energy and explicit
interval of eccentricities that starts from zero and reaches a physical value
is given in [8]. This is based on a computer assisted construction of
‘correctly aligned windows’. The difference is that in the current paper we
only need to establish the intersections of the manifolds for the unperturbed
system, and check some explicit conditions from [6]. The shadowing of orbits
is then automatically taken care of by [4], and we do not need to carry out
the explicit construction of windows as in [8]. Our results are weaker though
than [8].

The geometric setup for our method follows very closely that from [9, 10].
There are some small differences though, the main being that we use Poincaré
maps and compute finite Melnikov sums, instead of Melnikov integrals. Our
original plan was to follow directly the setup from [10], but we have found
that computing sharp bounds on integrals along trajectories is more difficult
than computing bounds on sums along discrete trajectories of Poincaré maps.

Our tool for the computer assisted implementation is the CAPD333Computer
Assisted Proofs in Dynamics: http://capd.ii.uj.edu.pl library [11]. The tools
which we use are quite standard: The existence of the family of Lyapunov
orbits as well as the proof of their homoclinic orbits is done by exploiting
symmetry properties of the PCR3BP, combined with parallel shooting and the
Krawczyk method. To establish bounds on the stable and unstable manifolds and
the transversality of their intersections needed to ensure that the scattering
map is well defined we use cones. The assumptions of theorems [6] leading to
diffusion can then be checked by computing sums along finite fragments of
homoclinic orbits.

The paper is organised as follows. Section 2 contains preliminaries, which
include the Krawczyk method, introduction to the scattering map and the
shadowing results for scattering maps, as well as a short introduction to the
restricted three body problem. In section 3 we state the main result of the
paper, which is written in Theorem 12. Section 4 contains the proof of the
main theorem. Some of the technical issues are delegated to the Appendix.

##  2 Preliminaries

###  2.1 Notations

For a set $U$ in a topological space we will denote its interior by
$\mathrm{int}U$, its closure by $\bar{U}$ and its boundary by $\partial U$.

We shall denote identity by $Id$. For a point $p$ expressed in some
coordinates $p=\left(x,y\right)$ we shall denote by $\pi_{x}\mathbf{\ }$and
$\pi_{y}$ the projections onto the given coordinates, i.e. $\pi_{x}p=x$ and
$\pi_{y}p=y$. For $p=\left(p_{1},\ldots,p_{n}\right)$ we will write $\pi_{i}$
to denote the projection onto the $i$-th coordinate, i.e. $\pi_{i}p=p_{i}$.

We write $\left\|\cdot\right\|_{\max}$ for the maximum norm in
$\mathbb{R}^{n}$, and for a matrix $A$ write $\left\|A\right\|_{\max}$ for the
matrix norm induced by $\left\|\cdot\right\|_{\max}$.

We shall denote a $k$-dimensional torus as $\mathbb{T}^{k}$, with the
convention that $\mathbb{T}:=\mathbb{R}\,\mathrm{mod}\,2\pi$.

###  2.2 Krawczyk method

We refer to a cartesian product of intervals as an interval set. For a set
$U\subset\mathbb{R}^{n}$ we shall denote by $\left[U\right]$ an interval set
such that $U\subset\left[U\right].$ We refer to such set as an interval
enclosure of $U$. An interval enclosure is not unique. In our applications we
shall consider interval enclosures of objects of interest; namely: fixed
points, periodic orbits, homoclinic orbits and invariant manifolds. The
smaller the interval enclosure is, the more accurate is our bound on the
object of interest. We shall refer to a matrix whose coefficients are
intervals as an interval matrix.

For a $C^{1}$ function $F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ we denote
by $\left[DF\left(X\right)\right]\subset\mathbb{R}^{n\times n}$ the interval
matrix enclosure of the derivatives of $F$ on $X$, namely we consider

$\left[DF\left(X\right)\right]=\left\\{A=\left(a_{i,j}\right):a_{ij}\in\left[\inf_{x\in
X}\frac{\partial f_{i}}{\partial x_{j}}\left(x\right),\sup_{x\in
X}\frac{\partial f_{i}}{\partial x_{j}}\left(x\right)\right]\right\\}.$

Below theorem, known as the Krawczyk method, can be used to establish bounds
on zeros of functions.

###### Theorem 1

[12]Let $F:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ be a $C^{1}$ function. Let
$X\subset\mathbb{R}^{n}$ be an interval set, let $x\in X$, let
$C\in\mathbb{R}^{n\times n}$ be a linear isomorphism and let

$K\left(x,X,F\right):=x-CF\left(x\right)+\left(Id-C\left[DF\left(X\right)\right]\right)\left(X-x\right).$

If

$K\left(x,X,F\right)\subset\mathrm{int}X,$

then there exists a unique point $x^{\ast}\in X$ for which

$F\left(x^{\ast}\right)=0.$

###  2.3 Normally hyperbolic invariant manifolds and the scattering map

In this section we recall the results which we use for our diffusion
mechanism. We follow the setup from [6], which is based on the scattering map
theory described in [2], and the diffusion mechanism from [4], which is the
main tool for our proof.

###### Definition 2

Let $M$ be a smooth $n$-dimensional manifold, and let $f:M\to M$ be a $C^{r}$
diffeomorphism, with $r>1$. Let $\Lambda\subset M$ be a compact manifold
without boundary, invariant under $f$, i.e., $f(\Lambda)=\Lambda$. We say that
$\Lambda$ is a normally hyperbolic invariant manifold (with symmetric rates)
if there exists a constant $C>0,$ rates $0<\lambda<\mu^{-1}<1$ and a
$Tf$-invariant splitting for every $x\in\Lambda$

$T_{x}M=E_{x}^{u}\oplus E_{x}^{s}\oplus T_{x}\Lambda$

such that

$\displaystyle v$ $\displaystyle\in
E_{x}^{u}\Leftrightarrow\left\|Df^{k}(x)v\right\|\leq
C\lambda^{-k}\left\|v\right\|,\quad k\leq 0,$ (1) $\displaystyle v$
$\displaystyle\in E_{x}^{s}\Leftrightarrow\left\|Df^{k}(x)v\right\|\leq
C\lambda^{k}\left\|v\right\|,\quad k\geq 0,$ (2) $\displaystyle v$
$\displaystyle\in T_{x}\Lambda\Rightarrow\left\|Df^{k}(x)v\right\|\leq
C\mu^{|k|}\left\|v\right\|,\quad k\in\mathbb{Z}.$ (3)

Let $d\left(x,\Lambda\right)$ stand for the distance between a point $x$ and
the manifold $\Lambda$. Given a normally hyperbolic invariant manifold and a
suitable small tubular neighbourhood $U\subset M$ of $\Lambda$ one defines its
local unstable and local stable manifold [13] as

$\displaystyle W_{\Lambda}^{u}\left(f,U\right)$ $\displaystyle=\left\\{y\in
M\,|\,f^{k}(y)\in U,d\left(f^{k}(y),\Lambda\right)\leq
C_{y}\lambda^{\left|k\right|},\,k\leq 0\right\\},$ $\displaystyle
W_{\Lambda}^{s}\left(f,U\right)$ $\displaystyle=\left\\{y\in M\,|\,f^{k}(y)\in
U,d\left(f^{k}(y),\Lambda\right)\leq C_{y}\lambda^{k},\,k\geq 0\right\\},$

where $C_{y}$ is a positive constant, which can depend on $y$. We define the
(global) unstable and stable manifolds as

$W_{\Lambda}^{u}\left(f\right)=\bigcup_{n\geq
0}f^{n}\left(W_{\Lambda}^{u}\left(f,U\right)\right),\qquad
W_{\Lambda}^{s}\left(f\right)=\bigcup_{n\geq
0}f^{-n}\left(W_{\Lambda}^{s}\left(f,U\right)\right).$

The manifolds $W_{\Lambda}^{u}\left(f,U\right)$,
$W_{\Lambda}^{s}\left(f,U\right)$, $W_{\Lambda}^{u}\left(f\right)$ and
$W_{\Lambda}^{s}\left(f\right)$ are foliated by

$\displaystyle W_{x}^{u}\left(f,U\right)$ $\displaystyle=\left\\{y\in
M\,|\,f^{k}(y)\in U,d(f^{k}(y),f^{k}(x))\leq
C_{x,y}\lambda^{\left|k\right|},\,k\leq 0\right\\},$ $\displaystyle
W_{x}^{s}\left(f,U\right)$ $\displaystyle=\left\\{y\in M\,|\,f^{k}(y)\in
U,d(f^{k}(y),f^{k}(x))\leq C_{x,y}\lambda^{k},\,k\geq 0\right\\},$

where $x\in\Lambda$ and $C_{x,y}$ is a positive constant, which can depend on
$x$ and $y$,

$W_{x}^{u}\left(f\right)=\bigcup_{n\geq
0}f^{n}\left(W_{f^{-n}\left(x\right)}^{u}\left(f,U\right)\right),\qquad
W_{x}^{s}\left(f\right)=\bigcup_{n\geq
0}f^{-n}\left(W_{f^{n}\left(x\right)}^{s}\left(f,U\right)\right).$

Let

$l<\min\left\\{r,\frac{|\log\lambda|}{\log\mu}\right\\}.$ (4)

The manifold $\Lambda$ is $C^{l}$ smooth, the manifolds
$W_{\Lambda}^{u}\left(f\right),W_{\Lambda}^{s}\left(f\right)$ are $C^{l-1}$
and $W_{x}^{u}\left(f\right)$, $W_{x}^{s}\left(f\right)$ are $C^{r}$ [2].
Normally hyperbolic manifolds, as well as their stable and unstable manifolds
and their fibres persist under small perturbations [13].

From now let $(M,\omega)$ be a smooth symplectic manifold. Let us assume that
$\Lambda\subset M$ is a normally hyperbolic invariant manifold for a $C^{r}$
symplectic map $f:M\to M$, where $r>1$. We assume that $\Lambda$ is even
dimensional and symplectic with the symplectic form $\omega|_{\Lambda}$, and
that $f|_{\Lambda}$ is symplectic on $\Lambda$. We define two maps,

$\displaystyle\Omega_{+}$
$\displaystyle:W_{\Lambda}^{s}(f)\rightarrow\Lambda,$
$\displaystyle\Omega_{-}$
$\displaystyle:W_{\Lambda}^{u}\left(f\right)\rightarrow\Lambda,$

where $\Omega_{+}(x)=x_{+}$ iff $x\in W_{x_{+}}^{s}\left(f\right)$, and
$\Omega_{-}(x)=x_{-}$ iff $x\in W_{x_{-}}^{u}\left(f\right).$ These are
referred to as the wave maps.

###### Definition 3

We say that a manifold $\Gamma\subset W_{\Lambda}^{u}\left(f\right)\cap
W_{\Lambda}^{s}\left(f\right)$ is a homoclinic channel for $\Lambda$ if the
following conditions hold:

  * (i)

for every $x\in\Gamma$

$\displaystyle
T_{x}W_{\Lambda}^{s}\left(f\right)+T_{x}W_{\Lambda}^{u}\left(f\right)$
$\displaystyle=T_{x}M,$ (5) $\displaystyle
T_{x}W_{\Lambda}^{s}\left(f\right)\cap T_{x}W_{\Lambda}^{u}\left(f\right)$
$\displaystyle=T_{x}\Gamma,$ (6)

  * (ii)

the fibres of $\Lambda$ intersect $\Gamma$ transversally in the following
sense

$\displaystyle T_{x}\Gamma\oplus T_{x}W_{x_{+}}^{s}\left(f\right)$
$\displaystyle=T_{x}W_{\Lambda}^{s}\left(f\right),$ (7) $\displaystyle
T_{x}\Gamma\oplus T_{x}W_{x_{-}}^{u}\left(f\right)$
$\displaystyle=T_{x}W_{\Lambda}^{u}\left(f\right),$ (8)

for every $x\in\Gamma$,

  * (iii)

the wave maps $(\Omega_{\pm})|_{\Gamma}:\Gamma\rightarrow\Lambda$ are
diffeomorphisms onto their image.

###### Definition 4

Assume that $\Gamma$ is a homoclinic channel for $\Lambda$ and let

$\Omega_{\pm}^{\Gamma}:=\left(\Omega_{\pm}\right)|_{\Gamma}.$

We define a scattering map $\sigma^{\Gamma}$ for the homoclinic channel
$\Gamma$ as

$\sigma^{\Gamma}:=\Omega_{+}^{\Gamma}\circ\left(\Omega_{-}^{\Gamma}\right)^{-1}:\Omega_{-}^{\Gamma}\left(\Gamma\right)\rightarrow\Omega_{+}^{\Gamma}\left(\Gamma\right).$

We have the following theorem, which is the basis for the diffusion mechanism
from the subsequent section.

###### Theorem 5

[4] Assume that $f:M\rightarrow M$ is a sufficiently smooth map,
$\Lambda\subset M$ is a normally hyperbolic invariant manifold with stable and
unstable manifolds which intersect transversally along a homoclinic channel
$\Gamma\subset M,$ and $\sigma$ is the scattering map associated to $\Gamma$.

Let $m_{1},\ldots,m_{l}\in\mathbb{N}$ be a fixed sequence of integers. Let
$\left\\{x_{i}\right\\}_{i=0,\ldots,l}$ be a finite pseudo-orbit in $\Lambda$,
that is a sequence of points in $\Lambda$ of the form

$x_{i+1}=f^{m_{i}}\circ\sigma^{\Gamma}\left(x_{i}\right),\qquad
i=0,\ldots,l-1,\,l\geq 1.$ (9)

Then for every $\delta>0$ there exists an orbit
$\left\\{z_{i}\right\\}_{i=0,\ldots,l}$ of $f$ in $M$, with
$z_{i+1}=f^{k_{i}}\left(z_{i}\right)$ for some $k_{i}>0$, such that
$d\left(z_{i},x_{i}\right)<\delta$ for all $i=0,\ldots,l$.

###### Remark 6

The result can be immediately extended to the case where we have a finite
number of scattering maps $\sigma_{1},\ldots,\sigma_{L}$ to shadow

$x_{i+1}=f^{m_{i}}\circ\sigma_{\alpha_{i}}\left(x_{i}\right),\qquad
i=0,\ldots,l-1,\,l\geq 1,$

for two prescribed sequences $m_{1},\ldots,m_{l}\in\mathbb{N}$ and
$\alpha_{1},\ldots,\alpha_{l}\in\left\\{1,\ldots,L\right\\}$; see [4, Theorem
3.7].

###### Remark 7

Theorem 5 has been simplified taking advantage of our particular setting. See
comments in [6, Remark 6].

###  2.4 Diffusion mechanism

In this section we recall the results from [6]. These are based on [4], but
were adapted in [6] to allow for computer assisted validation of the required
assumptions.

Consider the following family of Hamiltonian systems
$H_{\varepsilon}:\mathbb{R}^{4}\times\mathbb{T}\rightarrow\mathbb{R}$, that
depends smoothly on the parameter $\varepsilon\in\mathbb{R}$, which generates
the following ODE in the extended phase space

$\displaystyle x^{\prime}$
$\displaystyle=J\nabla_{x}H_{\varepsilon}\left(x,t\right),$ (10)
$\displaystyle t^{\prime}$ $\displaystyle=1,$

where

$J=\left(\begin{array}[c]{ll}0&Id\\\
-Id&0\end{array}\right),\qquad\text{for}\qquad
Id=\left(\begin{array}[c]{ll}1&0\\\ 0&1\end{array}\right).$

Let $\Phi_{t}^{\varepsilon}$ be the flow of (10). We refer to
$H_{\varepsilon=0}$ as the unperturbed system. We make the following important
assumption

$H_{\varepsilon=0}\left(x,t\right)=H\left(x\right),$

where $H:\mathbb{R}^{4}\rightarrow\mathbb{R}$; in other words, we assume that
the unperturbed system is autonomous and hence $H$ is a constant of motion.

Consider a Poincaré section $\Sigma$ in $\mathbb{R}^{4}$ for the system
$x^{\prime}=J\nabla H\left(x\right)$, and consider a section
$\tilde{\Sigma}=\Sigma\times\mathbb{T}$ for the perturbed system (10) in the
extended phase space. Let us consider a Pioncaré map

$f_{\varepsilon}:\tilde{\Sigma}\rightarrow\tilde{\Sigma}.$ (11)

We assume that for $\varepsilon=0$ the map $f_{\varepsilon=0}$ has a normally
hyperbolic invariant manifold $\Lambda_{0}\subset\tilde{\Sigma}$.

The coordinates on $\tilde{\Sigma}$ can be identified with
$\mathbb{R}^{3}\times\mathbb{T}$. We will write $\theta\in\mathbb{T}$ for the
extended phase space (time) coordinate. We assume that one of the remaining
coordinates on $\tilde{\Sigma}$ is the energy, expressed by the Hamiltonian
$H$. We will write $I$ for this coordinate. Since for $\varepsilon=0$ the
energy is preserved we have

$\pi_{I}f_{0}\left(z\right)=\pi_{I}z\qquad\text{for all }z\in\tilde{\Sigma}.$
(12)

We assume that $\Lambda_{0}$ is parameterised by $\theta,I$. For the remaining
two coordinates on $\tilde{\Sigma}$ we will write $u$ and $s$ and equip
$\tilde{\Sigma}$ with the standard symplectic from $\omega=du\wedge
ds+d\theta\wedge dI$. Due to the Hamiltonian nature of the system, the map
$f_{\varepsilon}$ is symplectic. We assume that

$\Lambda_{0}=\left\\{\left(u,s,I,\theta\right):u=s=0\right\\},$

is a normally hyperbolic invariant manifold for $f_{0}$. The coordinates $u$
and $s$ stand for the ‘unstable’ and ‘stable’ coordinates of $\Lambda_{0}$,
respectively. We assume also that for sufficiently small $\varepsilon$ the
manifold $\Lambda_{0}$ is perturbed to $\Lambda_{\varepsilon}$, which contains
a Cantor set of one dimensional KAM invariant tori.

Let
$g:\mathbb{R}\times\mathbb{R}^{3}\times\mathbb{T}\rightarrow\mathbb{R}^{3}\times\mathbb{T}$
be defined as

$g\left(\varepsilon,x\right):=\frac{1}{\varepsilon}\left(f_{\varepsilon}(x)-f_{0}(x)\right).$

Then

$f_{\varepsilon}(x)=f_{0}(x)+\varepsilon g(\varepsilon,x).$

The following theorem provides conditions under which for any sufficiently
small $\varepsilon>0$ there exists a point $x_{\varepsilon}$ and a number of
iterates $n_{\varepsilon}$ for which

$\pi_{I}\left(f_{\varepsilon}^{n_{\varepsilon}}\left(x_{\varepsilon}\right)-x_{\varepsilon}\right)>1.$
(13)

We first give a definition and follow with the statement of the theorem.

Figure 1: The setting for Theorem 9. [6]

###### Definition 8

Consider the topology on
$\Lambda_{0}\cap\left\\{I\in\left[0,1\right]\right\\}$ induced by
$\Lambda_{0}$. We say that an open set
$S\subset\Lambda_{0}\cap\left\\{I\in\left[0,1\right]\right\\}$ is a strip in
$\Lambda_{0}$ iff

$S\cap\left\\{z\in\Lambda_{0}:\pi_{I}z=\iota\right\\}\neq\emptyset\qquad\text{for
any }\iota\in\left[0,1\right].$

###### Theorem 9

[6] Assume that there is a neighborhood $U$ of $\Lambda_{0}$ and a positive
constant $L_{g}$ such that for every $z\in\Lambda_{0},x_{u}\in
W_{z}^{u}(f_{0},U),x_{s}\in W_{z}^{s}(f_{0},U),$

$\begin{array}[c]{c}\left|\pi_{I}\left(g(0,x_{u})-g\left(0,z\right)\right)\right|\leq
L_{g}\left\|x_{u}-z\right\|,\vskip 6.0pt plus 2.0pt minus 2.0pt\\\
\left|\pi_{I}\left(g(0,x_{s})-g\left(0,z\right)\right)\right|\leq
L_{g}\left\|x_{s}-z\right\|.\end{array}$ (14)

Assume also that there exist positive constants $C,\lambda$, where
$\lambda\in(0,1)$, such that for every $z\in\Lambda_{0}$ and every $x_{u}\in
W_{z}^{u}(f_{0},U),x_{s}\in W_{z}^{s}(f_{0},U)$ we have

$\begin{array}[c]{lll}\left\|f_{0}^{n}\left(z\right)-f_{0}^{n}\left(x_{u}\right)\right\|<C\lambda^{\left|n\right|}&&\text{for
all }n\leq 0,\vskip 6.0pt plus 2.0pt minus 2.0pt\\\
\left\|f_{0}^{n}\left(z\right)-f_{0}^{n}\left(x_{s}\right)\right\|<C\lambda^{n}&&\text{for
all }n\geq 0.\end{array}$ (15)

Assume also that for $\varepsilon=0$ we have the sequence of scattering maps
$\sigma_{\alpha}:\mathrm{dom}\left(\sigma_{\alpha}\right)\rightarrow\Lambda_{0}$
for $\alpha=1,\ldots,L$. Let $S^{+}\subset\Lambda_{0}$ be a strip. Assume that
for every $z\in\overline{S^{+}}$

  1. 1.

there exists an $\alpha\in\left\\{1,\ldots,L\right\\}$ for which
$z\in\mathrm{dom}\left(\sigma_{\alpha}\right)$, and there exists a constant
$m\in\mathbb{N}$ such that

$f_{0}^{m}\circ\sigma_{\alpha}\left(z\right)\in S^{+},$ (16)

  2. 2.

there exists a point $x\in W_{z}^{u}\left(f_{0},U\right)\cap
W_{\sigma_{\alpha}\left(z\right)}^{s}\left(f_{0}\right)$ such that
$f_{0}^{m}\left(x\right)\in
W_{f_{0}^{m}\left(\sigma_{\alpha}\left(z\right)\right)}^{s}(f_{0},U)$ and

$\sum_{j=0}^{m-1}\pi_{I}g\left(0,f_{0}^{j}\left(x\right)\right)-\frac{1+\lambda}{1-\lambda}L_{g}C>0.$
(17)

(The choice of $m$ and $\alpha$ can depend on $z$.)

Then for sufficiently small $\varepsilon>0$ there exists an $x_{\varepsilon}$
and $n_{\varepsilon}>0$ such that

$\pi_{I}\left(f_{\varepsilon}^{n_{\varepsilon}}\left(x_{\varepsilon}\right)-x_{\varepsilon}\right)>1.$

The following theorem can be used to establish orbits of the perturbed system,
whose $I$ coordinate decreases.

###### Theorem 10

[6]Assume that conditions (14) and (15) are satisfied, and that for
$\varepsilon=0$ we have the sequence of scattering maps
$\sigma_{\alpha}:\mathrm{dom}\left(\sigma_{\alpha}\right)\rightarrow\Lambda_{0}$
for $\alpha=1,\ldots,L$. Let $S^{-}\subset\Lambda_{0}$ be a strip. Assume that
for every $z\in\overline{S^{-}}$

  1. 1.

there exists an $\alpha\in\left\\{1,\ldots,L\right\\}$ for which
$z\in\mathrm{dom}\left(\sigma_{\alpha}\right)$, and there exists a constant
$m\in\mathbb{N}$ such that

$f_{0}^{m}\circ\sigma_{\alpha}\left(z\right)\in S^{-},$

  2. 2.

there exists a point $x\in W_{z}^{u}\left(f_{0},U\right)\cap
W_{\sigma_{\alpha}\left(z\right)}^{s}\left(f_{0}\right)$ such that
$f_{0}^{m}\left(x\right)\in
W_{f_{0}^{m}\left(\sigma_{\alpha}\left(z\right)\right)}^{s}(f_{0},U)$ and

$\sum_{j=0}^{m-1}\pi_{I}g\left(0,f_{0}^{j}\left(x\right)\right)+\frac{1+\lambda}{1-\lambda}L_{g}C<0.$

(The choice of $m$ and $\alpha$ can depend on $z$.)

Then for sufficiently small $\varepsilon>0$ there exists an $x_{\varepsilon}$
and $n_{\varepsilon}>0$ such that

$\pi_{I}\left(x_{\varepsilon}-f_{\varepsilon}^{n_{\varepsilon}}\left(x_{\varepsilon}\right)\right)>1.$

By combining the two strips $S^{+}$ and $S^{-}$ we obtain shadowing of any
prescribed finite sequence of actions.

###### Theorem 11

[6] Assume that two strips $S^{+}$ and $S^{-}$ satisfy assumptions of Theorems
9 and 10, respectively. If in addition

  1. 1.

for every $z\in\overline{S^{+}}$ there exists an $n$ (which can depend on $z$)
such that $f_{0}^{n}(z)\in S^{-}$, and

  2. 2.

for every $z\in\overline{S^{-}}$ there exists an $n$ (which can depend on $z$)
such that $f_{0}^{n}(z)\in S^{+}$,

then there exists an $\mathcal{M}$ such that for any given finite sequence
$\\{I_{k}\\}_{k=0}^{N}$ and any given $\delta>0$, for sufficiently small
$\varepsilon$ there exists an orbit of $f_{\varepsilon}$ which
$\varepsilon\mathcal{M}$-shadows the actions $I_{k}$; i.e. there exists a
point $z_{0}^{\varepsilon}$ and a sequence of integers
$n_{1}^{\varepsilon}\leq n_{2}^{\varepsilon}\leq\ldots\leq
n_{N}^{\varepsilon}$ such that

$\left\|\pi_{I}f_{\varepsilon}^{n_{k}^{\varepsilon}}\left(z_{0}^{\varepsilon}\right)-I_{k}\right\|<\varepsilon\mathcal{M}.$

###  2.5 The Planar Elliptic Restricted Three Body Problem

The Planar Elliptic Restricted Three Body Problem (PER3BP) describes the
motion of a massless particle (e.g., an asteroid or a spaceship), under the
gravitational pull of two large bodies, which we call primaries. The primaries
rotate in a plane along Keplerian elliptical orbits with eccentricity
$\varepsilon$, while the massless particle moves in the same plane and has no
influence on the orbits of the primaries. We use normalized units, in which
the masses of the primaries are $\mu$ and $1-\mu$. We consider a frame of
‘pulsating’ coordinates that rotates together with the primaries, making their
position fixed on the horizontal axis [14]. The motion of the massless
particle is described via the Hamiltonian
$H_{\varepsilon}:\mathbb{R}^{4}\times\mathbb{T\rightarrow R}$

$\begin{split}H_{\varepsilon}\left(X,Y,P_{X},P_{Y},\theta\right)&=\frac{\left(P_{X}+Y\right)^{2}+\left(P_{Y}-X\right)^{2}}{2}-\frac{\Omega\left(X,Y\right)}{1+\varepsilon\cos(\theta)},\\\
\Omega\left(X,Y\right)&=\frac{1}{2}\left(X^{2}+Y^{2}\right)+\frac{\left(1-\mu\right)}{r_{1}}+\frac{\mu}{r_{2}},\\\
r_{1}^{2}&=\left(X-\mu\right)^{2}+Y^{2},\\\
r_{2}^{2}&=\left(X-\mu+1\right)^{2}+Y^{2}.\end{split}$ (18)

The corresponding Hamilton equations are:

$\begin{array}[c]{lll}\displaystyle\frac{dX}{d\theta}=\frac{\partial
H_{\varepsilon}}{\partial P_{X}},\vskip 3.0pt plus 1.0pt minus
1.0pt&&\displaystyle\frac{dP_{X}}{d\theta}=-\frac{\partial
H_{\varepsilon}}{\partial X},\\\
\displaystyle\frac{dY}{d\theta}=\frac{\partial H_{\varepsilon}}{\partial
P_{Y}},&&\displaystyle\frac{dP_{Y}}{d\theta}=-\frac{\partial
H_{\varepsilon}}{\partial Y},\end{array}$ (19)

where $X,Y\in\mathbb{R}$ are the position coordinates of the massless
particle, and $P_{X},P_{Y}\in\mathbb{R}$ are the associated linear momenta.
The variable $\theta\in\mathbb{T}$ is the true anomaly of the Keplerian orbits
of the primaries, where $\mathbb{T}$ denotes the $1$-dimensional torus. The
system is non-autonomous, thus we consider it in the extended phase space, of
dimension $5$, with $\theta$ as an independent variable. We use the notation
$\Phi_{t}^{\varepsilon}$ to denote the flow of (19) in the extended phase
space, which includes $\theta\in\mathbb{T}$, i.e.

$\Phi_{t}^{\varepsilon}:\mathbb{R}^{4}\times\mathbb{T}^{1}\rightarrow\mathbb{R}^{4}\times\mathbb{T}^{1}.$

When $\varepsilon=0$ the primaries rotate around the center of mass along
circular orbits. The PER3BP becomes the Planar Circular Restricted Three Body
Problem (PCR3BP). We shall use the notation $\Phi_{t}$ for the flow of the
PCR3BP. Since the ODE of the PCR3BP is autonomous, this flow is not in the
extended phase space, i.e. for a given $t\in\mathbb{R}$,

$\Phi_{t}:\mathbb{R}^{4}\rightarrow\mathbb{R}^{4}.$

We naturally have the relation
$\Phi_{t}\left(x\right)=\pi_{x}\Phi_{t}^{0}\left(x,\theta\right)$ (for
$\varepsilon=0$ the right hand side does not depend on the choice of
$\theta$).

In the PCR3BP the system has the following time reversing symmetry

$\mathcal{R}\circ\Phi_{t}=\Phi_{-t}\circ\mathcal{R},$ (20)

where

$\mathcal{R}\left(X,Y,P_{X},P_{Y}\right)=\left(X,-Y,-P_{X},P_{Y}\right).$

We say that an orbit $x\left(t\right)$ is $\mathcal{R}$-symmetric if

$x(-t)=\mathcal{R}\left(x(t)\right).$

Note that for such orbits we need to have
$\pi_{Y}x\left(0\right)=\pi_{P_{X}}x\left(0\right)=0$.

The Jacobi integral $-2\cdot H_{0}$ is a preserved coordinate for the
unperturbed system with $\varepsilon=0.$ We shall refer to $H_{0}$ as the
energy. For a point $x\in\mathbb{R}^{4}$ we shall write

$I(x):=H_{0}\left(x\right).$

##  3 Statement of the main result

The PCR3BP has five libration fixed points. Three of these are on $\\{Y=0\\}$;
we refer to them as collinear and denote as $L_{1},L_{2},L_{3}$. We use the
convention, in which in the $X,Y$ coordinates the Jupiter, with mass $\mu$
lies to the left of the origin at $\left(\mu-1,0\right)$, and the Sun, with
mass $1-\mu$ is to the right of the origin at $\left(\mu,0\right)$. One of the
collinear libration points, which we denote as $L_{1}$ lies between the Sun
and Jupiter. Around this fixed point we have a family of Lyapunov periodic
orbits. Each orbit has different energy. The Lyapunov orbits are
$\mathcal{R}$-symmetric and we can choose points belonging to the orbits of
the form $\left(X,0,0,P_{X}\right)$, with suitable choices of $X$ and $P_{X}$.
The $P_{X}$ depends on the choice of $X$, so we shall write
$P_{X}=P_{X}\left(X\right)$.

The Lyapunov orbits considered by us will be parameterised by $X$. From now
on, we will write $\mathcal{X}\in\mathbb{R}$ to denote the value
$X=\mathcal{X}$, which determines the choice of a point on the Lyapunov orbit
$\left(\mathcal{X},0,0,P_{X}\left(\mathcal{X}\right)\right)$. This will allow
us to make the distinction between $\mathcal{X}$, which are values
parameterising the orbits and the coordinate $X$. We will be interested in the
family of Lyapunov orbits with

$\mathcal{X}\in\left[-0.95,-0.95+10^{-9}\right].$ (21)

The energy of such orbits, measured by the Jacobi integral, is approximately
$3.03$, which is the energy of the comet Oterma [15] that was observed in the
Jupiter-Sun system. We choose to work with this particular energy, since it
has a physical meaning of a known celestial body, but we could easily work
with a different $\mathcal{X}$-interval than (21), since there is nothing
special about the energy of Oterma. In fact, choosing a different energy level
can make some technical aspects of the computer assisted proof easier. (We
elaborate on this in Remark 46.) We fix an energy level of some physical
object to show that the method can be applied in a concrete setting.

For each $\mathcal{X}$ from (21) the corresponding Lyapunov orbit has a
different energy $I(\mathcal{X})$. We will show that for sufficiently small
$\varepsilon>0$, that is for the PER3BP, we can visit the arbitrary
$I(\mathcal{X})$ for $\mathcal{X}$ from (21). We shall refer to these as
diffusing orbits. Formally, the main result of this paper is as follows.

###### Theorem 12

(Main theorem) Let $I(\mathcal{X})$ denote the energy of a Lyapunov orbit
starting from a given Lyapunov orbit, i.e.

$I\left(\mathcal{X}\right)=H\left(\mathcal{X},0,0,P_{X}(\mathcal{X})\right).$

There exists a constant $\mathcal{M}>0$ such that for an arbitrary (finite)
sequence $\mathcal{X}_{1},\ldots,\mathcal{X}_{N}$ from the interval (21) and
for sufficiently small $\varepsilon>0$, there exists a sequence
$t_{1}^{\varepsilon},\ldots,t_{n}^{\varepsilon}$ and a point $x^{\varepsilon}$
such that

$\left\|H\left(\Phi_{t_{i}^{\varepsilon}}^{\varepsilon}(x^{\varepsilon})\right)-I\left(\mathcal{X}_{i}\right)\right\|<\varepsilon\mathcal{M}\qquad\qquad\text{for
}i=1,\ldots,N.$

###### Remark 13

The interval (21) is very narrow. The actual physical distance between the two
points on the Lyapunov orbits on the section $\\{Y=0\\}$, which correspond to
endpoints of (21) is about $1$ km. This means that the diffusion established
by us is over a very narrow range, but it is not completely negligible.

###### Remark 14

The proof was performed with computer-assisted tools, and for the interval
(21) it took 17 minutes, running on a single thread on a standard laptop. Our
computer assisted proof could be performed on other $\mathcal{X}$-intervals,
which combined together would lead to diffusion over longer distances. Such
proof could be performed by parallel computations on a cluster. We are more
interested in the proof of concept rather than in obtaining long intervals by
brute force, so we have not performed such validation.

##  4 Proof of the main result

The proof is performed in the following steps:

  1. 1.

establish the existence of the family of Lyapunov orbits, which will form the
NHIM $\Lambda_{0}$,

  2. 2.

establish bounds for the local stable/unstable manifolds of $\Lambda_{0}$,

  3. 3.

prove that the stable and unstable manifolds of $\Lambda_{0}$ intersect
transversally, and that we have a homoclinic channel along the intersection,

  4. 4.

prove that $\Lambda_{0}$ is perturbed to $\Lambda_{\varepsilon}$, which
contains a Cantor set of KAM tori.

  5. 5.

apply Theorem 11 to obtain the existence of diffusing orbits.

Section 4.1 sets out a parallel shooting method, which is then used for the
steps 1 and 3 in sections 4.2 and 4.4, respectively. Sections 4.3, 4.5 and 4.6
deal with steps 2, 4 and 5, respectively.

###  4.1 Parallel shooting for symmetric orbits

Here we consider the PCR3BP and work with the flow $\Phi_{t}$ in
$\mathbb{R}^{4}$. Let us consider a $C^{1}$ function
$p:\mathbb{R}\rightarrow\mathbb{R}^{4}$. For now we leave $p$ unspecified.
Later on, depending on the choice of $p$ we will use the below method to
establish either bounds on the points along a Lyapunov orbit, or for a
homoclinic orbit to a Lyapunov orbit.

Let us define the following function

$F:\mathbb{R}\times\mathbb{R}\times\underset{n}{\underbrace{\mathbb{R}^{4}\times\ldots\times\mathbb{R}^{4}}}\rightarrow\underset{n}{\underbrace{\mathbb{R}^{4}\times\ldots\times\mathbb{R}^{4}}}\times\mathbb{R}\times\mathbb{R},$

as

$\displaystyle F\left(s,\tau,x_{1},\ldots,x_{n}\right):=$ (22)
$\displaystyle\left(\Phi_{\tau}\left(p(s)\right)-x_{1},\Phi_{\tau}\left(x_{1}\right)-x_{2},\ldots,\Phi_{\tau}\left(x_{n-1}\right)-x_{n},\pi_{Y}\Phi_{\tau}\left(x_{n}\right),\pi_{P_{X}}\Phi_{\tau}\left(x_{n}\right)\right).$

We see that if we find a point
$\mathbf{x}^{\ast}=\left(s,\tau,x_{1},\ldots,x_{n}\right)$ such that

$F\left(\mathbf{x}^{\ast}\right)=0,$ (23)

then by taking $x_{0}=p(s)$ and $x_{n+1}=\Phi_{\tau}\left(x_{n}\right)$ we
obtain

$\Phi_{\tau}\left(x_{i}\right)=x_{i+1}\qquad\text{for }i=0,\ldots,n.$

Moreover, since

$\pi_{Y}x_{n+1}=\pi_{Y}\Phi_{\tau}\left(x_{n}\right)=0,\qquad\qquad\pi_{P_{X}}x_{n+1}=\pi_{P_{X}}\Phi_{\tau}\left(x_{n}\right)=0,$

we see that $x_{n+1}$ is $\mathcal{R}$-symmetric, i.e.
$x_{n+1}=\mathcal{R}\left(x_{n+1}\right).$ If we now define
$x_{n+k}=\mathcal{R}\left(x_{n+2-k}\right)$, for $k=2,\ldots,n+2$ then by (20)
we will obtain an $\mathcal{R}$-symmetric orbit $x(t)$, which starts from
$x\left(0\right)=x_{0}$, and for which

$\Phi_{\tau}\left(x_{i}\right)=x_{i+1}\qquad\text{for }i=0,\ldots,2n+1.$

Solving of (23) can be done by means of the Krawczyk method. The first step is
to obtain an approximate solution of (23) by iterating (using non-rigorous
numerics)

$\mathbf{x}_{i+1}=\mathbf{x}_{i}-\left(DF\left(\mathbf{x}_{i}\right)\right)^{-1}F\left(\mathbf{x}_{i}\right).$
(24)

After a few iterates of (24) we obtain a point $\mathbf{x}$ around which we
can construct a cube $\mathbf{X}$ and validate that we have the solution of
(23) inside of $\mathbf{X}$ by using Theorem 1. (We take $C$ as the non-
rigorously computed inverse of the derivative
$C=\left(DF\left(\mathbf{x}\right)\right)^{-1}$.)

###  4.2 Bounds for Lyapunov orbits

When we fix some $\mathcal{X}\in\mathbb{R}$ and choose the function
$p:\mathbb{R}\rightarrow\mathbb{R}^{4}$ as

$p\left(s\right)=\left(\mathcal{X},0,0,s\right),$ (25)

then the methodology form section 4.1 can be used to obtain a sequence of
points $x_{0},\ldots,x_{2n+2}$ along an $\mathcal{R}$-symmetric periodic
orbit. This is because by the definition of $F$

$\Phi_{(n+1)\tau}\left(x_{0}\right)=x_{n+1},$

and since $x_{0}$ and $x_{n+1}$ are self $\mathcal{R}$-symmetric, by (20)

$\Phi_{(n+1)\tau}\left(x_{n+1}\right)=\Phi_{(n+1)\tau}\circ\mathcal{R}\left(x_{n+1}\right)=\mathcal{R}\circ\Phi_{-(n+1)\tau}\left(x_{n+1}\right)=\mathcal{R}\left(x_{0}\right)=x_{0}.$

We thus see that $x_{0}$ is a point on an $\mathcal{R}$-symmetric periodic
orbit with period $T=2\left(n+1\right)\tau$.

An advantage of the method is that we can obtain a bound for a whole family of
Lyapunov orbits. This can be done by considering an interval instead of a
single $\mathcal{X}$, for our validation. The computer assisted computation
allows us to obtain a bound on $F(x)$ and $DF\left(x\right)$ for (22). This
way we obtain a bound for points which lie on Lyapunov orbits for an interval
of values $\mathcal{X}$. In our case we take the interval (21) and use this
method to validate the following result:

$n$ | $X$ | $Y$ | $P_{X}$ | $P_{Y}$
---|---|---|---|---
0 | -0.9499999995 | 0 | 0 | -0.84134724633
1 | -0.95011002908 | 0.010872319337 | -0.012750492306 | -0.84566628682
2 | -0.95027977734 | 0.020942127841 | -0.021798848297 | -0.85701977629
3 | -0.95016249921 | 0.02964511208 | -0.02595601236 | -0.87222441368
4 | -0.94945269037 | 0.036681290981 | -0.026117965229 | -0.88881047965
5 | -0.94799596417 | 0.041912835694 | -0.023738329652 | -0.90554260324
6 | -0.94578584443 | 0.045275924414 | -0.020027448005 | -0.92194698888
7 | -0.94292354313 | 0.046741951127 | -0.015815875271 | -0.93784624844
8 | -0.93958019632 | 0.046310536541 | -0.011642040196 | -0.9531124603
9 | -0.93596962345 | 0.044016227704 | -0.0078569213076 | -0.96756323289
10 | -0.93232909834 | 0.039939546482 | -0.0046935571959 | -0.98092506113
11 | -0.92890403848 | 0.034218254755 | -0.0022984843476 | -0.99282861416
12 | -0.92593321085 | 0.027056429274 | -0.00073280826007 | -1.0028268728
13 | -0.92363228182 | 0.018728871729 | 4.5125264188e-05 | -1.0104387689
14 | -0.92217533629 | 0.0095779013399 | 0.00019712976878 | -1.0152203731
15 | -0.92167641746 | 0 | 0 | -1.0168530766
Table 1: Midpoints of our enclosure of points along the family of Lyapunov
orbits, computed using (23) with $n=14$.

###### Lemma 15

Let $r:=3.633\cdot 10^{-9}$. For every $\mathcal{X}$ from the interval (21)
there exists an

$s(\mathcal{X})\in-0.84134724633+\left[-r,r\right]$

such that

$x_{0}(\mathcal{X})=\left(\mathcal{X},0,0,s(\mathcal{X})\right)$ (26)

is a point, which lies on a Lyapunov orbit, which we denote as
$L_{\mathcal{X}}$. Moreover, we have a sequence of points along
$L_{\mathcal{X}}$, which passes within the $r$ distance (in maximum norm) of
the points from Table 1. Moreover, we have the following bound
$T(\mathcal{X})$ for the period of $L_{\mathcal{X}}$

$T(\mathcal{X})\in[3.0417517493,3.0417517846]$ (27)

###### Remark 16

In Table 1 we write out half of the points along the family of periodic
orbits, since the second half (the remaining fourteen points, to be precise)
follows from the $\mathcal{R}$-symmetry.

###### Remark 17

Our estimate on the distance from the points from Table 1 obtained by our
computer program varies from point to point and is frequently more accurate
than stated in Lemma 15, where we simply write out the single $r$, which is
the upper bound that can be applied to all the points.

###  4.3 Bounds on the unstable manifolds of Lyapunov orbits

In the previous section we have shown how to compute the bound on the family
of Lyapunov orbits $L_{\mathcal{X}}$ containing points $x_{0}(\mathcal{X})$
given by (26), with $\mathcal{X}$ from the interval (21). Here we fix a single
periodic orbit for some $\mathcal{X}$ from (21) and discuss how one can obtain
a computer validated enclosure of its local unstable manifold. Before we give
the method, we introduce some notation.

We consider a Poincaré section $\Sigma=\left\\{Y=0\right\\}$ and define
$\rho:\mathbb{R}^{4}\rightarrow\mathbb{R}$ and
$\mathcal{P}:\mathbb{R}^{4}\rightarrow\mathbb{R}^{4}$ as

$\displaystyle\rho\left(x\right)$
$\displaystyle=\inf\left\\{t>0:\Phi_{t}\left(x\right)\in\Sigma\right\\},$
$\displaystyle\mathcal{P}\left(x\right)$
$\displaystyle=\Phi_{\rho\left(x\right)}\left(x\right).$

In other words, $\rho$ is the time along the flow to the section, and
$\mathcal{P}$ is the map which goes to the section along the flow.

###### Remark 18

The $\rho$ and $\mathcal{P}$ do not need to be globally defined. Whenever we
will use these functions in our computer assisted proofs, the CAPD library
verifies that the considered sets lie within the domains of these maps, and
that they are properly defined throughout the performed validations. (If a set
would not belong to the domain, the program would return an error and
terminate.)

Let us fix $\mathcal{X}$ from (21) and a Lyapunov orbit $L_{\mathcal{X}}$
containing $x_{0}=x_{0}(\mathcal{X})$ given by (26). Denote the period of this
orbit as $T=T(\mathcal{X})$. For convenience, let us write $N=2n+2$, so that
for our points $x_{0},\ldots,x_{N}$ on the Lyapunov orbit, whose bounds we
established in Lemma 15 we have

$\displaystyle\Phi_{T/N}\left(x_{i}\right)$
$\displaystyle=x_{i+1}\qquad\text{for }i=0,\ldots,N-1,$ (28) $\displaystyle
x_{N}$ $\displaystyle=x_{0}.$

Note that from (26) we see that $x_{0}\in\Sigma$.

Let us consider now a sequence of invertible matrices
$A_{i}\in\mathbb{R}^{4\times 4}$ for $i=0,\ldots,N$, with $A_{0}=A_{N}$ and
define the following maps

$f_{i}:\mathbb{R}^{4}\rightarrow\mathbb{R}^{4},\qquad\text{for }i=1,\ldots,N,$

as

$\displaystyle f_{i}\left(v\right)$
$\displaystyle:=A_{i}^{-1}\left(\Phi_{T/N}\left(A_{i-1}v+x_{i-1}\right)-x_{i}\right),\qquad\text{for
}i=1,\ldots,N-1,$ $\displaystyle f_{N}\left(v\right)$
$\displaystyle:=A_{N}^{-1}\left(\mathcal{P}\left(A_{N-1}v+x_{N-1}\right)-x_{N}\right)=A_{0}^{-1}\left(\mathcal{P}\left(A_{N-1}v+x_{N-1}\right)-x_{0}\right).$
(29)

In other words, for $i=1,\ldots,N-1$ we consider time shift maps along the
flow, expressed in local coordinates around $x_{i}$. The last map, $f_{N}$,
maps to the section $\Sigma=\left\\{Y=0\right\\}\subset\mathbb{R}^{4}$. This
means that $f_{N}\circ\ldots\circ f_{1}|_{\Sigma}:\Sigma\rightarrow\Sigma$.
From (28) it follows that

$f_{i}\left(0\right)=0\qquad\text{for }i=1,\ldots,N.$

For $F:\mathbb{R}^{4}\rightarrow\mathbb{R}^{4}$ defined as

$F=f_{N}\circ\ldots\circ f_{1}$ (30)

we see that the origin is a fixed point. Our objective will be to establish
bounds on the unstable manifold of the origin. To be more precise, we shall
establish bounds on the intersection of the unstable manifold of the Lyapunov
orbit with $\Sigma$, which is the unstable manifold of the origin for the map
$F|_{\Sigma}$.

First we introduce some notion.

###### Definition 19

Let $\left\|\cdot\right\|$ be some norm in $\mathbb{R}^{3}.$ Let
$Q:\mathbb{R}^{4}\rightarrow\mathbb{R}$ be the function

$Q\left(v_{1},\ldots,v_{4}\right)=\left|v_{1}\right|-\left\|\left(v_{2},v_{3},v_{4}\right)\right\|.$
(31)

We define the cone centered at a point $v\in\mathbb{R}^{4}$ as

$Q^{+}\left(v\right):=\left\\{w\in\mathbb{R}^{4}:Q\left(w-v\right)\geq
0\right\\}.$

We consider a sequence of cones defined by
$Q_{i}:\mathbb{R}^{4}\rightarrow\mathbb{R}$, for $i=0,\ldots,N$ and assume
that

$Q_{N}=Q_{0}.$

We take the norms $\left\|\cdot\right\|_{i}$ for $Q_{i}$ from (31) as

$\left\|\left(x_{1},x_{2},x_{3}\right)\right\|_{i}=\max\left\\{\left|x_{1}\right|/a_{i,1},\ldots,\left|x_{3}\right|/a_{i,3}\right\\}$

where $a_{i,k}\in\left(0,1\right)$ are fixed coefficients for $i=1,\ldots,N$
and $k=1,2,3$. In other words, we use different norms in (31) to define
different cones. Note that
$\left|y\right|\geq\left\|\left(x_{1},x_{2},x_{3}\right)\right\|_{i}$ is
equivalent to $a_{i,k}\left|y\right|\geq\left|x_{k}\right|$, for $k=1,2,3$, so
the cone $Q_{i}^{+}\left(v\right)$ can be expressed as

$Q_{i}^{+}\left(v\right)=\left\\{v+\left(t,tx_{1},tx_{2},tx_{3}\right):x_{k}\in\left[-a_{i,k},a_{i,k}\right]\text{
for }k=1,2,3\text{ and }t\in\mathbb{R}\right\\}.$ (32)

###### Remark 20

The form (32) is convenient, since then cones defined by $Q_{i}^{+}$ can be
represented in a computer assisted implementation by a set

$V_{i}=\left[1\right]\times\left[-a_{i,1},a_{i,1}\right]\times\ldots\times\left[-a_{i,3},a_{i,3}\right].$

What we mean by this is that

$Q_{i}^{+}\left(v\right)=\left\\{v+tw:w\in V_{i}\text{ and
}t\in\mathbb{R}\right\\}.$ (33)

###### Definition 21

Let $B\subset\mathbb{R}^{4}$. We say that $f_{i}$ satisfies cone conditions in
$B$ iff for every $v\in B$

$f_{i}\left(Q_{i}^{+}\left(v\right)\cap B\right)\subset
Q_{i+1}^{+}\left(f_{i}(v)\right).$ (34)

Below lemma is our main tool for establishing bounds on the unstable manifold
of the origin for the map (30).

###### Lemma 22

[16, Lemma 6.3] Let $B:=\left[-1,1\right]^{4}\subset\mathbb{R}^{4}$. Assume
that $f_{1},\ldots,f_{N}$ satisfy cone conditions in $B$. Let $m>1$ and assume
that for $F=f_{N}\circ\ldots\circ f_{1}$ the matrix $DF\left(0\right)$ has a
single eigenvalue $\lambda$ satisfying
$\left|\operatorname{Re}\lambda\right|>m$ and the absolute values of the real
parts of the remaining eigenvalues below $m$. If also for every $v\in
Q_{i}^{+}\left(0\right)\cap B$,

$\left\|f_{i}\left(v\right)\right\|>m\left\|v\right\|,$ (35)

then the unstable manifold of the origin for the map $F$ is parameterised as a
smooth curve $p^{\mathrm{u}}:\left[-1,1\right]\rightarrow B$ which satisfies

$\displaystyle p^{\mathrm{u}}\left(0\right)$ $\displaystyle=0,$
$\displaystyle\pi_{1}p^{\mathrm{u}}$ $\displaystyle=Id,$ $\displaystyle
p^{\mathrm{u}}\left(\left[-1,1\right]\right)$ $\displaystyle\subset
Q_{0}^{+}\left(p^{\mathrm{u}}\left(u\right)\right),\qquad\text{for every
}u\in\left[-1,1\right],$

and

$\frac{d}{du}p^{\mathrm{u}}\left(u\right)\in
Q_{0}^{+}\left(0\right),\qquad\text{for every }u\in\left[-1,1\right].$

###### Remark 23

In [16, Lemma 6.3] is stated for the setting where $F$ is a full turn along
the periodic orbit. Here we have a sequence of local maps, which when composed
constitute the full turn. The proof in such setting is analogous to [16], by
using the graph transform method. The only needed modification with respect to
[16] is that here the graphs need to be propagated successively by
$f_{1},f_{2},\ldots,f_{N},$ after which they return to the same local
coordinates. The graph transform for each successive $f_{i}$ follows from an
identical construction as the one from [16] (which there is done for the
single $F$). By composing the graph transforms for the successive $f_{i}$ we
obtain a graph transform for $F$.

###### Remark 24

The benefit of considering several maps and shooting is that this requires
shorter integration times, which improves accuracy in the interval arithmetic
computations. This is the only reason why we shoot between 29 points along
Lyapunov orbits, instead of considering a single turn.

###### Remark 25

In our computer assisted proof, we choose $A_{0}$ so that the derivative of
$F=f_{N}\circ\ldots\circ f_{1}$ at the origin is close to diagonal. The
eigenvalues of $F$ are $\lambda,\frac{1}{\lambda},1,0$. (The zero comes from
the fact that $f_{N}$ maps to the section $\Sigma$.) We can validate the bound
on $\lambda$ by using the Gersgorin theorem. We choose the remaining $A_{i}$
so that the derivatives of $f_{i}$ at the origin are close to diagonal.

###### Remark 26

The validation of cone conditions is done as follows. Take an interval set
$V_{i}$ for which we have (33). Then by the mean value theorem the fact that

$\left[Df_{i}\left(B\right)\right]V_{i}\subset Q_{i+1}\left(0\right)$ (36)

implies (34). To check (36) it is enough to compute the interval set
$w=\left[Df_{i}\left(B\right)\right]V_{i}$ and validate that

$\left[\frac{w}{\pi_{1}w}\right]\subset V_{i+1}.$

This means that the cone condition is straightforward to validate from the
interval enclosure of the derivative of the map.

In our computer assisted proof we use the following lemma to validate (35) for
the maximum norm.

###### Lemma 27

Consider
$Q\left(v_{1},\ldots,v_{4}\right)=\left|v_{1}\right|-\left\|\left(v_{2},v_{3},v_{4}\right)\right\|$
with

$\left\|\left(x_{1},x_{2},x_{3}\right)\right\|=\max\left\\{\left|x_{1}\right|/a_{1},\left|x_{1}\right|/a_{2},\left|x_{3}\right|/a_{3}\right\\},$

where $a_{1},a_{2},a_{3}\in\left(0,1\right)$. Take
$f:B\rightarrow\mathbb{R}^{4}$ such that $f\left(0\right)=0$ and let

$\left[Df\left(B\right)\right]=\left(\begin{array}[c]{cc}A_{11}&A_{12}\\\
A_{21}&A_{22}\end{array}\right),$ (37)

where $A_{11}$, $A_{12}$, $A_{21}$ and $A_{22}$ are $1\times 1$, $1\times 3$,
$3\times 1$ and $3\times 3$ interval matrices, respectively. If $A_{11}>c>0$
and $m=c-\left\|A_{12}\right\|_{\max}\max\left(a_{1},a_{2},a_{3}\right)$, then

$\left\|f\left(v\right)\right\|_{\max}\geq
m\left\|v\right\|_{\max}\qquad\text{for every }v\in Q^{+}\left(0\right)\cap
B.$

Proof. The proof is given in A.

A mirror result can be formulated to obtain a bound in the other direction.
This is relevant for the method for obtaining the bound (14) outlined in C. We
place this lemma here since it follows from similar arguments to Lemma 22.

###### Lemma 28

Consider
$Q\left(v_{1},\ldots,v_{4}\right)=\left|v_{1}\right|-\left\|\left(v_{2},v_{3},v_{4}\right)\right\|$
with

$\left\|\left(x_{1},x_{2},x_{3}\right)\right\|=\max\left\\{\left|x_{1}\right|/a_{1},\left|x_{1}\right|/a_{2},\left|x_{3}\right|/a_{3}\right\\},$

where $a_{1},a_{2},a_{3}\in\left(0,1\right)$. Take
$f:B\rightarrow\mathbb{R}^{4}$ and let

$\left[Df\left(B\right)\right]=\left(\begin{array}[c]{cc}A_{11}&A_{12}\\\
A_{21}&A_{22}\end{array}\right),$ (38)

where $A_{11}$, $A_{12}$, $A_{21}$ and $A_{22}$ are $1\times 1$, $1\times 3$,
$3\times 1$ and $3\times 3$ interval matrices, respectively. If

$\displaystyle a$ $\displaystyle:=\max\left(a_{1},a_{2},a_{3}\right),$ (39)
$\displaystyle\bar{m}$
$\displaystyle:=\max\left(\left|A_{11}\right|+a\left\|A_{12}\right\|_{\max},\left\|A_{21}\right\|_{\max}+a\left\|A_{22}\right\|_{\max}\right),$
(40)

then

$\left\|f\left(v\right)\right\|_{\max}\leq\bar{m}\left\|v\right\|_{\max}\qquad\text{for
every }v\in Q^{+}\left(0\right)\cap B.$

Proof. The proof is given in B.

###### Remark 29

In practice, the term $A_{11}$ in (37)–(38), is associated with hyperbolic
expansion, and dominates in (40). Also, in our application the number $a$ from
(39) is small, so we obtain $m\approx\bar{m}\approx\left|A_{11}\right|$;
naturally $m<\bar{m}$.

With the aid of Lemma 22 we have validated the following result.

###### Lemma 30

Let $r=3\cdot 10^{-8}$, $B=[-r,r]^{4}$, and let

$A_{0}=\,\left(\begin{array}[c]{cccc}0.280324&-0.220733&0.280324&0\\\
0&0&0&0.816632\\\ 1&0&-1&-1\\\ -0.343269&1&-0.343269&0\end{array}\right).$
(41)

Consider $L=5\cdot 10^{-6}$ and cones defined as

$Q_{0}^{+}\left(v\right)=\left\\{v+\left(t,tx_{2},tx_{3},tx_{4}\right):x_{k}\in\left[-L,L\right]\text{
for }k=2,3,4\text{ and }t\in\mathbb{R}\right\\}.$

Then the unstable manifold of the origin for the map $F$ is parameterised as a
smooth curve $p^{\mathrm{u}}:\left[-r,r\right]\rightarrow B$ which satisfies

$\displaystyle p^{\mathrm{u}}\left(0\right)$ $\displaystyle=0,$
$\displaystyle\pi_{1}p^{\mathrm{u}}$ $\displaystyle=Id,$ $\displaystyle
p^{\mathrm{u}}\left(\left[-r,r\right]\right)$ $\displaystyle\subset
Q_{0}^{+}\left(p^{\mathrm{u}}\left(u\right)\right),\qquad\text{for every
}u\in\left[-r,r\right],$

and

$\frac{d}{du}p^{\mathrm{u}}\left(u\right)\in
Q_{0}^{+}\left(0\right),\qquad\text{for every }u\in\left[-r,r\right].$

###### Corollary 31

The curve
$w^{\mathrm{u}}\left(u\right):=x_{0}+A_{0}p^{\mathrm{u}}\left(u\right)$ lies
on the unstable manifold of the Lyapunov orbit containing $x_{0}.$ The curve
$w^{\mathrm{u}}\left(u\right)$ is contained in
$x_{0}+A_{0}Q_{0}^{+}\left(0\right)$ and
$\frac{d}{du}w^{\mathrm{u}}\left(u\right)\in A_{0}V_{0}$ where

$V_{0}=\left\\{1\right\\}\times\left[-L,L\right]\times\left[-L,L\right]\times\left[-L,L\right].$

Since to define $F$ in (30) we take $f_{N}=\mathcal{P}$ (see (29)), and
$\mathcal{P}$ maps to $\Sigma$, we know that
$w^{\mathrm{u}}\left(u\right)\subset\Sigma$.

###### Remark 32

The matrices $A_{i}$ used for the coordinate changes as well as the cones
$Q_{i}$ are chosen automatically by our computer program. We do not write them
out here, since the important bounds for the proof are at the point
$x_{0}=x_{0}(\mathcal{X})$, and are given in Lemma 30.

###### Remark 33

To validate assumption (35) of Lemma 22 we have used Lemmas 27 and obtaining
$m=1.2767546773,$ for the bound from (35).

###### Remark 34

We have also obtained the following bound by using Lemma 28 for maps (29)

$\left\|f_{i}\left(v_{i}\right)\right\|_{\max}\leq\bar{m}\left\|v_{i}\right\|_{\max}\qquad\text{for
}v_{i}\in Q_{i}^{+}\left(0\right)\cap B\text{ and }i=1,\ldots,N,$

where $\bar{m}=1.2767636743.$

So far we have discussed how to obtain a bound on the fiber for a fixed point
at the origin of our section-to-section map $F$ defined in (30). Such fixed
point resulted from intersecting the Lyapunov orbit with $\\{Y=0\\}$. In the
future discussion we will need to consider the problem in the extended phase
space, since the PER3BP is not autonomous. In the extended phase space the
fixed point at the origin for $F$ becomes an invariant curve
$\left\\{0\right\\}\times\mathbb{T}$. The unstable fiber for a given point
$\left(0,\lambda\right)$, with $\lambda\in\mathbb{T}$, on this curve is
contained in the extended phase space. Since the return times to the section
$\\{Y=0\\}$ differ from point to point, such fiber does not need to be
contained in $\left\\{\theta=\lambda\right\\}$. Below we discuss a method with
which we establish bounds on the unstable fibers of points from
$\left\\{0\right\\}\times\mathbb{T}$ in the extended phase space.

To make the above discussion more precise, we consider the flow of the PCR3BP
in the extended phase space and denote it as $\tilde{\Phi}_{t}$. We consider
the Poincaré section
$\tilde{\Sigma}=\left\\{Y=0\right\\}=\Sigma\times\mathbb{T}$, define
$\tilde{\rho}:\mathbb{R}^{4}\times\mathbb{T}\rightarrow\mathbb{R}$ and define
$\mathcal{\tilde{P}}:\mathbb{R}^{4}\times\mathbb{T}\rightarrow\mathbb{R}^{4}\times\mathbb{T}$
as

$\displaystyle\tilde{\rho}\left(x\right)$
$\displaystyle=\inf\left\\{t>0:\tilde{\Phi}_{t}\left(x\right)\in\tilde{\Sigma}\right\\},$
$\displaystyle\mathcal{\tilde{P}}\left(x\right)$
$\displaystyle=\tilde{\Phi}_{\tilde{\rho}\left(x\right)}\left(x\right).$

Let also $\tilde{A}_{0}$ be a $5\times 5$ matrix defined as

$\tilde{A}_{0}=\left(\begin{array}[c]{cc}A_{0}&0\\\ 0&1\end{array}\right),$
(42)

where $A_{0}$ is from (41). We define
$\tilde{F}:\tilde{\Sigma}\rightarrow\tilde{\Sigma}$ as

$\tilde{F}\left(x,\theta\right)=\tilde{A}_{0}^{-1}\left(\mathcal{\tilde{P}}^{2}\left(\left(x_{0},0\right)+\tilde{A}_{0}\left(x,\theta\right)\right)-\left(x_{0},0\right)\right).$

Then $\left\\{0\right\\}\times\mathbb{T}$ becomes a an invariant curve for the
map $\tilde{F}$. We will now show how to obtain bounds for the unstable fibers
of a point $\left(0,\lambda\right)$ on such curve, in the extended phase
space.

###### Lemma 35

Let $r$ and $L$ be the constants considered in Lemma 30. Let $M\in\mathbb{R}$,
$M>0$. Consider cones in the extended phase space defined as

$\tilde{Q}_{0}^{+}\left(v\right)=\left\\{v+\left(t,tx_{2},tx_{3},tx_{4},t\theta\right):x_{k}\in\left[-L,L\right]\text{
for }k=2,3,4,\text{ }\theta\in\left[-M,M\right]\text{ and
}t\in\mathbb{R}\right\\}.$

If for every $\lambda\in\mathbb{T}$ the unstable eigenvector of
$D\tilde{F}\left(0,\lambda\right)$ is contained in
$\tilde{Q}_{0}^{+}\left(0\right)$ and if $\tilde{F}$ satisfies
$\tilde{Q}_{0}^{+}$ cone conditions on
$\left[-r,r\right]^{4}\times\mathbb{T}$, then for every $\lambda\in\mathbb{T}$
the unstable fiber $W_{\left(0,\lambda\right)}^{u}(\tilde{F})$ is
parameterised by
$\tilde{p}_{\lambda}^{\mathrm{u}}:\left[-r,r\right]\rightarrow\mathbb{R}^{4}\times\mathbb{T}$
satisfying (below $p^{\mathrm{u}}$ is the function established in Lemma 30)

$\pi_{\mathbb{R}^{4}}\tilde{p}_{\lambda}^{\mathrm{u}}\left(u\right)=p^{\mathrm{u}}\left(u\right)\qquad\text{and
\qquad}\tilde{p}_{\lambda}^{\mathrm{u}}\left(u\right)\in\tilde{Q}_{0}^{+}\left(0,\lambda\right)\qquad\mbox{for
}u\in[-r,r].$ (43)

In particular

$\left|\pi_{\theta}\tilde{p}_{\lambda}^{\mathrm{u}}\left(u\right)-\lambda\right|\leq
rM.$ (44)

For the family of Lyapunov orbits with $\mathcal{X}$ from the interval (21) we
can take $M=3.$

Proof. For every $\lambda$ the unstable fiber
$W_{\left(0,\lambda\right)}^{u}(\tilde{F})$ is tangent at
$\left(0,\lambda\right)$ to the eigenvector of
$D\tilde{F}\left(0,\lambda\right).$ This implies that sufficiently close to
$\left(0,\lambda\right)$ this fiber is in
$\tilde{Q}_{0}^{+}\left(0,\lambda\right)$. Every point $x\in
W_{\left(0,\lambda\right)}^{u}(\tilde{F})$ can be expressed as
$x=\tilde{F}^{m}\left(z\right)$ for an arbitrary $m\in\mathbb{N}$ , and for
some appropriate point $z=z(m)\in
W_{\tilde{F}^{-m}(0,\lambda)}^{u}(\tilde{F})$. By taking sufficiently large
$m$ the point $z$ can be chosen as close to
$\left\\{0\right\\}\times\mathbb{T}$ as we want, which means that we can
choose $m$ large enough so that
$z\in\tilde{Q}_{0}^{+}\left(\tilde{F}^{-m}(0,\lambda)\right).$ Since
$z\in\tilde{Q}_{0}^{+}\left(\tilde{F}^{-m}(0,\lambda)\right)$, by the fact
that $\tilde{F}$ satisfies cone condition we obtain that
$x=\tilde{F}^{m}\left(z\right)\in\tilde{Q}_{0}^{+}\left(0,\lambda\right).$
This means that for
$W_{\left(0,\lambda\right)}^{u}(\tilde{F})\cap(\left[-r,r\right]^{4}\times\mathbb{T})\subset\tilde{Q}_{0}^{+}\left(0,\lambda\right)$.
Since the PCR3BP is autonomous,
$\pi_{\mathbb{R}^{4}}W_{\left(0,\lambda\right)}^{u}(\tilde{F})$ does not
depend on $\lambda$, and is parameterised by $p^{\mathrm{u}}\left(u\right)$
from Lemma 30. We have thus established (43).

The condition (44) follows from the fact that
$\tilde{p}_{\lambda}^{\mathrm{u}}\left(u\right)\in\tilde{Q}_{0}^{+}\left(0,\lambda\right)$
by computing

$\left|\pi_{\mathbb{T}}\tilde{p}_{\lambda}^{\mathrm{u}}\left(u\right)-\lambda\right|\leq
M\left|\pi_{x_{1}}\tilde{p}_{\lambda}^{\mathrm{u}}\left(u\right)\right|=M\left|\pi_{x_{1}}p^{\mathrm{u}}\left(u\right)\right|=M\left|u\right|\leq
rM.$

With computer assistance, we have validated that for every $\mathcal{X}$ from
the interval (21), the unstable eigenvector of
$D\tilde{F}\left(0,\lambda\right)$ is contained in

$\\{1\\}\times[-3\cdot 10^{-6},3\cdot
10^{-6}]^{3}\times[2.7763408157,2.7918430312]\subset\tilde{Q}_{0}^{+}(0).$

We have also validated the cone conditions using the method described in
Remark 26.

###### Remark 36

By choosing $\tilde{A}_{0}$ more carefully, instead of just adding $1$ in the
lower diagonal term in (42), we could reduce the constant $M$ in Lemma 35. We
have opted for (42) for the sake of simplicity, since the established $M$ is
good enough for the proof of the main result.

###  4.4 Intersection of stable/unstable manifolds of Lyapunov orbits

The way we establish intersection of the stable and unstable manifolds of
Lyapunov orbits is similar to the method from section 4.2. We use the parallel
shooting from section 4.1 combined with the bounds on the unstable manifold
established in section 4.3. We fix $\mathcal{X}$ from (21) and consider
$p:\mathbb{R}\rightarrow\mathbb{R}^{4}$ for the shooting operator (22) to be

$p\left(s\right)=x_{0}(\mathcal{X})+A_{0}p^{\mathrm{u}}\left(s\right),$ (45)

where $p^{\mathrm{u}}$ is the function from which parameterises the
intersection of the unstable manifold of $L_{\mathcal{X}}$ with
$\left\\{Y=0\right\\}$, whose bounds we have obtained in Lemma 30. (The choice
of $A_{0}$ and the bounds on $p^{\mathrm{u}}\left(s\right),$ together with its
derivative are written out in Lemma 30).

$n$ | $X$ | $Y$ | $P_{X}$ | $P_{Y}$
---|---|---|---|---
0 | -0.95000000242 | -0 | -1.0427994645e-08 | -0.84134724275
1 | -0.94997599415 | 0.032508255048 | -0.026407880234 | -0.87841641186
2 | -0.94367569216 | 0.046559121408 | -0.016859552423 | -0.93400981538
3 | -0.93185570859 | 0.039269494574 | -0.0043284590814 | -0.98260192685
4 | -0.9227740512 | 0.014134265524 | 0.00018569728993 | -1.0132585163
5 | -0.92342538156 | -0.017741502784 | -8.9504023812e-05 | -1.0111198858
6 | -0.9332873018 | -0.041191398464 | 0.0054636769818 | -0.97749478375
7 | -0.94483412083 | -0.046015557539 | 0.018549243384 | -0.92767343656
8 | -0.95017205833 | -0.029473079334 | 0.025907589742 | -0.87187283933
9 | -0.95002477482 | 0.0043928225728 | -0.0053602054156 | -0.84205223945
10 | -0.94968890529 | 0.035271238551 | -0.026402429149 | -0.88508159957
11 | -0.94246936536 | 0.046807227358 | -0.015244331177 | -0.94027843983
12 | -0.93053821073 | 0.037102593271 | -0.003488949681 | -0.9875100719
13 | -0.92242610924 | 0.010464221796 | -0.00019611370513 | -1.0149474376
14 | -0.92456420774 | -0.021092391132 | -0.00083603746 | -1.0083700678
15 | -0.93548053074 | -0.042380742641 | 0.0047663014948 | -0.97130987256
16 | -0.94730307615 | -0.043772642856 | 0.016236744567 | -0.91841934035
17 | -0.95320906454 | -0.022318932226 | 0.012564889629 | -0.85729285711
18 | -0.96070049024 | 0.017216683341 | -0.061196504134 | -0.84499124342
19 | -0.98081741509 | 0.044653406262 | -0.11087155568 | -0.94184521723
20 | -1.0057774379 | 0.042684693077 | -0.12465612254 | -1.0610257866
21 | -1.0281206713 | 0 | 0 | -1.2112777154
Table 2: Midpoints of our enclosure o homoclinic orbits for our family of
Lyapunov orbits.

By the $\mathcal{R}$-symmetry of the PCR3BP we know that
$\mathcal{R}\left(p\left(s\right)\right)$ is a point on the stable manifold of
$L_{\mathcal{X}}$. If for $F$ defined by (22) we validate that for some
$s\in\left(0,1\right)$ and $h\in\mathbb{R}$ we have

$F\left(s,h,x_{1},\ldots,x_{n}\right)=0,$

then taking $x_{0}:=p\left(s\right)$ and $x_{n+1}:=\Phi_{h}\left(x_{n}\right)$
we obtain a sequence of points

$x_{0},\ldots,x_{n},x_{n+1},\mathcal{R}\left(x_{n}\right),\ldots,\mathcal{R}\left(x_{0}\right),$
(46)

along an intersection of the stable and unstable manifolds of
$L_{\mathcal{X}}$. With this method we have managed to obtain the following
result.

###### Lemma 37

For every $\mathcal{X}$ from (21) the stable and unstable manifolds of
$L_{\mathcal{X}}$ intersect. Moreover, the intersection is along a
$\mathcal{R}$-symmetric homoclinic orbit, which contains a sequence of points
along it that lies $r=1.96\cdot 10^{-7}$ close to the orbit written out in
Table 2; see also Figure 2.

###### Remark 38

In Table 2 we write out half of the points along the homoclinic, since the
second half (the remaining twenty one points, to be precise) follows from the
$\mathcal{R}$-symmetry.

###### Remark 39

Our estimate on the distance from the points from Table 2 obtained by our
computer program varies from point to point and is frequently more accurate
than stated in Lemma 37, where we simply write out the single $r$, which is
the upper bound that can be applied to all the points.

###### Remark 40

From the method we also obtain a bound on the integration time $h$ between the
consecutive points along the homoclinic. We have obtained that

$h\in[0.34246881126,0.34246888642].$ (47)

Figure 2: A homoclinic orbit to $L_{\mathcal{X}}.$

We now show that the established intersection is transversal.

###### Lemma 41

For every $\mathcal{X}$ from the interval (21) the intersection of the stable
and unstable manifolds of $L_{\mathcal{X}}$ is transversal when considered in
the three dimensional constant energy level.

Proof. We consider
$\Sigma_{\\{X<-1\\}}=\left\\{Y=0,X<-1\right\\}\subset\mathbb{R}^{4}$ and we
will study the intersections of the stable/unstable manifolds of
$L_{\mathcal{X}}$ on this section. (See Figure 2.)

Let us denote by $W_{L_{\mathcal{X}}}^{u}$ and $W_{L_{\mathcal{X}}}^{s}$ the
unstable and stable manifolds of $L_{\mathcal{X}},$ respectively. These are
two dimensional tubes, contained in the three dimensional constant energy
level $\left\\{H=h^{\ast}\right\\}$ for
$h^{\ast}=H\left(L_{\mathcal{X}}\right)$. We have established that
$W_{L_{\mathcal{X}}}^{u}$ and $W_{L_{\mathcal{X}}}^{s}$ intersect along a
homoclinic orbit, which passes through the points (46). The point $x_{n+1}$
belongs to $\Sigma_{\\{X<-1\\}}$. (From (22) we know that $\pi_{Y}x_{n+1}=0$
and from Table 2 we see that $\pi_{X}x_{n+1}<-1$, so
$x_{n+1}\in\Sigma_{\\{X<-1\\}}$.) The vector field at $x_{n+1}$ has a non zero
$Y$-component. This means that the tangent spaces to the stable and unstable
manifolds at $x_{n+1}$ span the coordinate $Y$.

The manifolds $W_{L_{\mathcal{X}}}^{u}$ and $W_{L_{\mathcal{X}}}^{s}$
intersect with $\Sigma_{\\{X<-1\\}}$ at $x_{n+1}$ along one dimensional
curves. (Note that some of the points from the unstable/stable manifolds can
collide with Jupiter. Those that reach $\Sigma_{\\{X<-1\\}}$ close to
$x_{n+1}$ intersect the section along one dimensional curves.) The section
$\Sigma_{\\{X<-1\\}}$ is three dimensional, but $L_{\mathcal{X}},$
$W_{L_{\mathcal{X}}}^{u}$ and $W_{L_{\mathcal{X}}}^{s}$ are contained in
$\left\\{H=h^{\ast}\right\\}$. The set
$\Sigma_{\\{X<-1\\}}\cap\left\\{H=h^{\ast}\right\\}$ is two dimensional and
can be parameterised444On $\Sigma_{\left\\{X<-1\right\\}}$ the coordinate
$P_{Y}$ can be computed from $X,P_{X}$ since $H(X,Y=0,P_{X},P_{Y})=h^{*}$. by
coordinates $(X,P_{X})$. The $W_{L_{\mathcal{X}}}^{u}\cap\Sigma_{\\{X<-1\\}}$
and $W_{L_{\mathcal{X}}}^{s}\cap\Sigma_{\\{X<-1\\}}$ are therefore one
dimensional curves contained in a two dimensional space, parameterised by
$(X,P_{X})$, and if we show that

$\pi_{X,P_{X}}\left(W_{L_{\mathcal{X}}}^{u}\cap\Sigma_{\\{X<-1\\}}\right)\text{\quad
intersect transversally
with\quad}\pi_{X,P_{X}}\left(W_{L_{\mathcal{X}}}^{s}\cap\Sigma_{\\{X<-1\\}}\right)\text{\quad
at }x_{n+1},$ (48)

then we will obtain transversal intersections of $W_{L_{\mathcal{X}}}^{u}$
with $W_{L_{\mathcal{X}}}^{s}$ at $x_{n+1}$ in $\left\\{H=h^{\ast}\right\\}.$
(In more detail: the vector field at $x_{n+1}$ is tangent to
$W_{L_{\mathcal{X}}}^{u}$ and $W_{L_{\mathcal{X}}}^{s}$ at $x_{n+1}$ and has a
non zero $Y$-component; from (48) we will have that the tangent spaces to
$W_{L_{\mathcal{X}}}^{u}$ and $W_{L_{\mathcal{X}}}^{s}$ at $x_{n+1}$ span
$X,P_{X}$. In all we span a three dimensional vector space, hence the
intersection is transversal in $\left\\{H=h^{\ast}\right\\}.$)

Let $\tau:\mathbb{R}^{4}\rightarrow\mathbb{R}$,
$\mathcal{P}:\mathbb{R}^{4}\rightarrow\Sigma_{\\{X<-1\\}}$ be defined as

$\displaystyle\tau\left(x\right)$
$\displaystyle:=\inf\left\\{t>0:\Phi_{t}\left(x\right)\in\Sigma_{\left\\{X<-1\right\\}}\right\\},$
$\displaystyle\mathcal{P}\left(x\right)$
$\displaystyle:=\Phi_{\tau\left(x\right)}\left(x\right).$

The curve $W_{L_{\mathcal{X}}}^{u}\cap\Sigma_{\\{X<-1\\}}$ can be obtained by
computing $\mathcal{P}\circ p\left(s\right)$. (See (45) for the definition of
$p\left(s\right).$) By the $\mathcal{R}$-symmetry of the PCR3BP the
$W_{L_{\mathcal{X}}}^{s}\cap\Sigma_{\\{X<-1\\}}$ is equal to
$\mathcal{R}\circ\mathcal{P}\circ p\left(s\right)$. Let
$s^{\ast}\in\mathbb{R}$ be such that $\mathcal{P}\circ
p\left(s^{\ast}\right)=x_{n+1}$. If we establish that

$\pi_{X}\frac{d}{ds}\mathcal{P}\circ
p\left(s\right)|_{s=s^{\ast}}>0\qquad\text{and\qquad}\pi_{P_{X}}\frac{d}{ds}\mathcal{P}\circ
p\left(s\right)|_{s=s^{\ast}}>0,$ (49)

then

$\displaystyle\pi_{X}\frac{d}{ds}\mathcal{R}\circ\mathcal{P}\circ
p\left(s\right)|_{s=s^{\ast}}$
$\displaystyle=\pi_{X}\frac{d}{ds}\mathcal{P}\circ
p\left(s\right)|_{s=s^{\ast}}>0,$ (50)
$\displaystyle\pi_{P_{X}}\frac{d}{ds}\mathcal{R}\circ\mathcal{P}\circ
p\left(s\right)|_{s=s^{\ast}}$
$\displaystyle=-\pi_{P_{X}}\frac{d}{ds}\mathcal{P}\circ
p\left(s\right)|_{s=s^{\ast}}<0,$ (51)

and from (49–51) we will obtain (48).

In section 4.3 we have established that (see Lemma 30 and Corollary 31)

$\frac{d}{ds}p\left(s\right)\in\tilde{A}_{0}V_{0},$ (52)

where $V_{0}$ is a set

$V_{0}=\left\\{1\right\\}\times\left[-L,L\right]\times\left[-L,L\right]\times\left[-L,L\right],$

with $L=5\cdot 10^{-6}$. The set $V_{0}$ represents a cone, as described in
Remark 20. We can propagate the bound (52) to the point $x_{n+1}$ using cone
propagation method described in Remark 26. We have thus validated that

$\pi_{X,P_{X}}\frac{d}{ds}\mathcal{P}\circ
p\left(s\right)|_{s=s^{\ast}}\in\left\\{uV\,|\,u>0\text{ and
}V=\\{1\\}\times[4.06081,4.06404]\right\\}.$

This establish (49) and finishes our proof.

###  4.5 Persistence of the family of Lyapunov orbits

Recall that a Lyapunov orbit starting from $x_{0}(\mathcal{X})$ (see (26)),
which has a period $T(\mathcal{X})$, is given as

$L_{\mathcal{X}}=\left\\{\Phi_{t}\left(x_{0}(\mathcal{X})\right):t\in\left[0,T(\mathcal{X})\right]\right\\}.$

Let us denote the normally hyperbolic invariant manifold consisting of the
family of Lyapunov orbits as

$\Lambda_{L}=\\{L_{\mathcal{X}}:\mathcal{X}\mbox{ are from the interval
(\ref{eq:X-interval})}\\}\subset\mathbb{R}^{4}.$

We shall also use the following notation for the manifold in the extended
phase space:

$\tilde{\Lambda}_{L}=\Lambda_{L}\times\mathbb{T}.$ (53)

To prove persistence of $\tilde{\Lambda}_{L}$ we shall use the following
theorem:

###### Theorem 42

[10] Assume that

$\frac{d}{d\mathcal{X}}T(\mathcal{X})\neq 0$ (54)

and also

$\frac{d}{d\mathcal{X}}H(x_{0}(\mathcal{X}))\neq 0.$ (55)

Then for sufficiently small perturbation $\varepsilon$ from the PCR3BP to the
PER3BP, the manifold $\tilde{\Lambda}_{L}$ is perturbed into a
$O\left(\varepsilon\right)$ close normally hyperbolic manifold
$\tilde{\Lambda}_{L}^{\varepsilon}$, with boundary, which is invariant under
the flow induced by (19). Moreover, there exists a Cantor set of invariant
tori in $\tilde{\Lambda}_{L}^{\varepsilon}$.

We now discuss how we validate (54–55). Let us fix a single $\mathcal{X}$ from
(26). From section 4.2 we know that there exists a $\tau=\tau(\mathcal{X})$,
and an even555In our case we take $N=30$, see Table 1. $N\in\mathbb{N}$, such
that for $\tau=\tau(\mathcal{X})$ and for

$\displaystyle x_{0}$
$\displaystyle=x_{0}\left(\mathcal{X}\right)=\left(\mathcal{X},0,0,s\left(\mathcal{X}\right)\right),$
$\displaystyle x_{i+1}$
$\displaystyle=\Phi_{\tau}\left(x_{i}\right)\qquad\text{for }i=0,\ldots,N-1$

we have

$x_{N}=x_{0},$

and the period of the Lyapunov orbit starting from $x_{0}$ is
$T(\mathcal{X})=N\tau\left(\mathcal{X}\right)$.

Above was established by fixing $\mathcal{X}$ and solving for $\tau$ and $s$
the following equation

$\pi_{Y,P_{X}}\Phi_{\tau N/2}\left(\mathcal{X},0,0,s\right)=0.$ (56)

(See (22) and (25); equation (22) provides us with the solution of (56) using
parallel shooting.)

We can now define a function $g:\mathbb{R}^{3}\rightarrow\mathbb{R}^{2}$ as

$g\left(\mathcal{X},s,\tau\right)=\pi_{Y,P_{X}}\Phi_{N\tau/2}\left(\mathcal{X},0,0,s\right)$

and observe that

$g\left(\mathcal{X},s(\mathcal{X}),\tau(\mathcal{X})\right)=0.$ (57)

This means that we can compute the derivatives of $\frac{ds}{d\mathcal{X}}$
and $\frac{d\tau}{d\mathcal{X}}$ from the implicit function theorem; i.e. by
differentiating (57) with respect to $\mathcal{X}$ we obtain

$\frac{\partial g}{\partial\mathcal{X}}+\frac{\partial
g}{\partial\left(s,\tau\right)}\left(\begin{array}[c]{c}\frac{ds}{d\mathcal{X}}\\\
\frac{d\tau}{d\mathcal{X}}\end{array}\right)=0,$

and provided that $\frac{\partial g}{\partial\left(s,\tau\right)}$ is
invertible we see that

$\left(\begin{array}[c]{c}\frac{ds}{d\mathcal{X}}\\\
\frac{d\tau}{d\mathcal{X}}\end{array}\right)=-\left(\frac{\partial
g}{\partial\left(s,\tau\right)}\right)^{-1}\frac{\partial
g}{\partial\mathcal{X}}.$ (58)

The partials $\frac{\partial g}{\partial\mathcal{X}},$ $\frac{\partial
g}{\partial s}$ and $\frac{\partial g}{\partial\tau}$ are $2\times 1$
matrices, which can be computed as follows. Let $e_{i}\in\mathbb{R}^{4}$, for
$i=1,\ldots,4$, be a vector with $1$ on $i$-th coordinate and zeros on the
remaining coordinates. Let $\mathcal{F}:\mathbb{R}^{4}\to\mathbb{R}^{4}$ stand
for the vector field of the PCR3BP. Then

$\displaystyle\frac{\partial g}{\partial\mathcal{X}}\left(x_{0}\right)$
$\displaystyle=\pi_{Y,P_{X}}D\Phi_{N\tau/2}\left(x_{0}\right)e_{1}=\pi_{Y,P_{X}}D\Phi_{\tau}\left(x_{N/2-1}\right)\ldots
D\Phi_{\tau}\left(x_{0}\right)e_{1},$ (59) $\displaystyle\frac{\partial
g}{\partial s}\left(x_{0}\right)$
$\displaystyle=\pi_{Y,P_{X}}D\Phi_{N\tau/2}\left(x_{0}\right)e_{4}=\pi_{Y,P_{X}}D\Phi_{\tau}\left(x_{N/2-1}\right)\ldots
D\Phi_{\tau}\left(x_{0}\right)e_{4},$ (60) $\displaystyle\frac{\partial
g}{\partial\tau}\left(x_{0}\right)$
$\displaystyle=\frac{d}{d\tau}\pi_{X,P_{X}}\Phi_{N\tau/2}(x_{0})=\frac{N}{2}\pi_{Y,P_{X}}\mathcal{F}\left(\Phi_{N\tau/2}\left(x_{0}\right)\right)=\frac{N}{2}\pi_{Y,P_{X}}\mathcal{F}\left(x_{N/2}\right).$
(61)

Note that from the above $\frac{\partial g}{\partial s}$ and $\frac{\partial
g}{\partial\tau}$ we obtain $\frac{\partial
g}{\partial\left(s,\tau\right)}=\left(\begin{array}[c]{cc}\frac{\partial
g}{\partial s}&\frac{\partial g}{\partial\tau}\end{array}\right)$, which we
can use in (58) to compute $\frac{ds}{d\mathcal{X}}$ and
$\frac{d\tau}{d\mathcal{X}}$.

Once $\frac{ds}{d\mathcal{X}}$ is established, we can easily compute

$\frac{dH}{d\mathcal{X}}(x_{0}(\mathcal{X}))=\frac{\partial H}{\partial
X}(x_{0}(\mathcal{X}))+\frac{\partial H}{\partial
P_{Y}}(x_{0}(\mathcal{X}))\frac{ds}{d\mathcal{X}}(\mathcal{X}).$ (62)

We have used (58–62) to validate, with computer assistance, that we have the
following:

###### Lemma 43

For every $\mathcal{X}$ from (21) we have

$\displaystyle\frac{d}{d\mathcal{X}}x_{0}\left(\mathcal{X}\right)$
$\displaystyle\in\\{1\\}\times\\{0\\}\times\\{0\\}\times\left[-4.530367,-4.530349\right],$
$\displaystyle\frac{d}{d\mathcal{X}}\tau(\mathcal{X})$
$\displaystyle\in\left[-0.5523403,-0.5522624\right],$
$\displaystyle\frac{d}{d\mathcal{X}}T(\mathcal{X})$
$\displaystyle\in\left[-16.57021,-16.56787\right],$
$\displaystyle\frac{d}{d\mathcal{X}}H\left(x_{0}\left(\mathcal{X}\right)\right)$
$\displaystyle\in\left[-0.359187,-0.359185\right].$

###### Corollary 44

By Theorem 42 and Lemma 43 we obtain the persistence of the normally
hyperbolic manifold $\tilde{\Lambda}_{L}$ (see (53)), which consists of our
family of Lyapunov orbits in the extended phase space.

###  4.6 Proof of the main theorem

Recall that $\Phi_{t}^{\varepsilon}$ is the flow of the PER3BP in the extended
phase space. Let
$\rho^{\varepsilon}:\mathbb{R}^{4}\times\mathbb{T}\rightarrow\mathbb{R}$ be
the time to the section $\left\\{Y=0\right\\}$

$\rho^{\varepsilon}\left(x\right)=\inf\left\\{t>0:\Phi_{t}^{\varepsilon}\left(x\right)\in\left\\{Y=0\right\\}\right\\},$

and let
$\mathcal{P}_{\varepsilon}:\mathbb{R}^{4}\times\mathbb{T}\rightarrow\left\\{Y=0\right\\}$
be defined as

$\mathcal{P}_{\varepsilon}\left(x\right)=\Phi_{\rho\left(x\right)}^{\varepsilon}\left(x\right).$
(63)

The section $\left\\{Y=0\right\\}$ in the extended phase space is isomorphic
with $\mathbb{R}^{3}\times\mathbb{T}$, which fits the setting from section
2.4.

We consider a single point $x^{\ast}:=x_{0}\left(-0.95\right)$
($x_{0}(\mathcal{X})$ is defined in (26); see also (21) regarding the choice
of $-0.95$) and consider the matrix $\tilde{A}_{0}$ from (42). We define

$\tilde{\Sigma}=\left\\{\tilde{A}_{0}^{-1}\left(x-x^{\ast},\theta\right):x\in\left\\{Y=0\right\\},\theta\in\mathbb{T}\right\\}.$

In other words, $\tilde{\Sigma}$ is the section $\left\\{Y=0\right\\}$
considered in the extended phase space, in the local coordinates given by the
affine change given by $x^{\ast}$ and $\tilde{A}_{0}$. We decide to work in
these local coordinates since then we can directly use the estimates on the
local unstable manifolds, which were established in section 4.3. This is the
reason why we choose $\tilde{\Sigma}$ as above.

To apply Theorem 11 we will choose our family of maps (11)

$f_{\varepsilon}:\tilde{\Sigma}\rightarrow\tilde{\Sigma}$

to be defined as

$f_{\varepsilon}\left(x,\theta\right)=\tilde{A}_{0}^{-1}\left(\mathcal{P}_{\varepsilon}\circ\mathcal{P}_{\varepsilon}\left(\left(x^{\ast},0\right)+\tilde{A}_{0}\left(x,\theta\right)\right)-\left(x^{\ast},0\right)\right).$
(64)

In other words, we consider the return map to the section $\\{Y=0\\}$
expressed in our local coordinates.

###### Remark 45

The reason for taking a composition of two $\mathcal{P}_{\varepsilon}$ in (64)
is that for $\varepsilon=0$ a Lyapunov orbit intersects $\\{Y=0\\}$ in two
points; or to be more precise, for $\varepsilon=0$ a Lyapunov orbit becomes an
invariant torus in the extended phase space, which intersects $\\{Y=0\\}$
along two disjoint curves. (See Figure 3.) When working with
$\mathcal{P}_{\varepsilon}\circ\mathcal{P}_{\varepsilon}$, each of these
curves becomes an invariant circle for $f_{\varepsilon=0}$.

Prior to the perturbation, for $\varepsilon=0$, we define our normally
hyperbolic invariant manifold $\Lambda_{0}\subset\tilde{\Sigma}$ as (see (53)
for the definition of $\tilde{\Lambda}_{L}$)

$\Lambda_{0}:=\tilde{\Lambda}_{L}\cap\tilde{\Sigma}.$

The manifold $\Lambda_{0}\ $is foliated by invariant circles for the map
$f_{0}$.

Figure 3: The manifold $\Lambda_{0}$ projected onto coordinates $(X,Y,\theta)$
on the left, and in coordinates $(I,\theta)$ on the right.

Note that for $\varepsilon=0$ the energy is preserved, so for

$I(x):=H_{0}\left(\left(x^{\ast},0\right)+\tilde{A}_{0}x\right)$ (65)

we see that

$I\left(f_{0}\left(x\right)\right)=I(x).$ (66)

We treat $I$ as our conserved variable when $\varepsilon=0$ and write

$\pi_{I}x=I(x).$

Thus, condition (66) plays the role of (12).

We can write $f_{\varepsilon}$ as

$f_{\varepsilon}\left(x\right)=f_{0}\left(x\right)+\varepsilon
g\left(\varepsilon,x\right)$

where

$g\left(\varepsilon,x\right)=\frac{1}{\varepsilon}\left(f_{\varepsilon}\left(x\right)-f_{0}\left(x\right)\right).$
(67)

To compute the change of the energy after an iterate of
$f_{\varepsilon}\left(x\right)$ we compute

$\pi_{I}f_{\varepsilon}\left(x\right)=\pi_{I}f_{0}\left(x\right)+\varepsilon\pi_{I}g\left(\varepsilon,x\right)=\pi_{I}x+\varepsilon\pi_{I}g\left(\varepsilon,x\right),$

where $\pi_{I}g\left(\varepsilon,x\right)$ is

$\pi_{I}g\left(\varepsilon,x\right)=\frac{1}{\varepsilon}\left[I\left(f_{\varepsilon}\left(x\right)\right)-I\left(f_{0}\left(x\right)\right)\right]=\frac{1}{\varepsilon}\left[I\left(f_{\varepsilon}\left(x\right)\right)-I\left(x\right)\right].$

It follows from the above that

$\pi_{I}g\left(0,x\right)=\frac{d}{d\varepsilon}I(f_{\varepsilon}(x))|_{\varepsilon=0}=DI\left(f_{0}\left(x\right)\right)\frac{\partial
f_{0}}{\partial\varepsilon}\left(x\right)=\nabla
I\left(f_{0}\left(x\right)\right)\cdot\frac{\partial
f_{0}}{\partial\varepsilon}\left(x\right).$ (68)

(In the above equation $\cdot$ is the scalar product.)

We are now ready for the proof of our main result.

Proof of Theorem 12. We start by describing our system for $\varepsilon=0$.
While discussing the system for $\varepsilon=0$ we will recall some results
for the PCR3BP established in the previous sections. We need to keep in mind
that these were considered in coordinates
$\left(X,Y,P_{X},P_{Y}\right)\in\mathbb{R}^{4}$; without the extended time
coordinate $\theta$.

The manifold $\Lambda_{0}$ is invariant under $f_{0}$ and in coordinates
$X,Y,P_{X},P_{Y},\theta$ can be written as

$\Lambda_{0}=\left\\{\left(x_{0}\left(\mathcal{X}\right),\theta\right):\text{
}\mathcal{X}\text{ is from (\ref{eq:X-interval}),
}\theta\in\mathbb{T}\right\\}\subset\mathbb{R}^{4}\times\mathbb{T}.$

For $\varepsilon=0$ the inner dynamics produced by $f_{0}$ on the manifold is
given as

$\left(x_{0}\left(\mathcal{X}\right),\theta\right)\mapsto\left(x_{0}\left(\mathcal{X}\right),\theta+T\left(\mathcal{X}\right)\right),$
(69)

where $T\left(\mathcal{X}\right)$ is the period of the Lyapunov orbit
$L_{\mathcal{X}}$. Thus $\Lambda_{0}$ is an invariant cylinder, foliated by
invariant curves.

Recall that for a given single Lyapunov orbit $L_{\mathcal{X}}$, we have
established in Lemma 30 the bounds on a curve $p^{\mathrm{u}}(u)$, with
$u\in\left(-r,r\right)$, which lies along the intersection of the two
dimensional local unstable manifold of $L_{\mathcal{X}}$ (for the PCR3BP in
$\mathbb{R}^{4}$) with the section $\\{Y=0\\}$. Let us emphasize the
dependence of $p^{\mathrm{u}}(u)$ on $\mathcal{X}$ by writing
$p_{\mathcal{X}}^{\mathrm{u}}\left(u\right)$. The unstable manifold of
$\Lambda_{0}$ for $f_{0}$, considered in $\\{Y=0\\}$ (in the extended phase
space) is three dimensional, and in the coordinates $X,P_{X},P_{Y},\theta$,
can locally be written as

$W_{\Lambda_{0}}^{u}=\left\\{\left(p_{\mathcal{X}}^{\mathrm{u}}\left(u\right),\theta\right):\text{
}\mathcal{X}\text{ is from (\ref{eq:X-interval}),
}u\in\left(-r,r\right),\theta\in\mathbb{T}\right\\}.$

By considering the $\mathcal{R}$-symmetry of the PCR3BP, in the extended phase
space, and restricted to $\\{Y=0\\}$, i.e.

$\mathcal{\tilde{R}}\left(X,P_{X},P_{Y},\theta\right):=\left(X,-P_{X},P_{Y},\theta\right),$

we obtain the local stable manifold

$W_{\Lambda_{0}}^{s}=\mathcal{\tilde{R}}\left(W_{\Lambda_{0}}^{u}\right).$

We will now show that for $\varepsilon=0$ we have a well defined (and global)
scattering map

$\sigma:\Lambda_{0}\rightarrow\Lambda_{0}.$ (70)

For this we first need to establish a homoclinic channel. By Lemmas 37, 41 we
know that the two dimensional stable and unstable manifolds of
$L_{\mathcal{X}}$ in the PCR3BP in $\mathbb{R}^{4}$ intersect transversally
(when considered on a three dimensional fixed energy set
$\left\\{H_{0}=I\right\\}$, where $I=H_{0}\left(L_{\mathcal{X}}\right)$) along
an $\mathcal{R}$-symmetric homoclinic orbit, which contains a point which we
shall denote here as $\gamma\left(I\right)\in\\{Y=0\\}$. The two dimensional
stable and unstable manifolds of a given Lyapunov orbit $L_{\mathcal{X}}$,
when intersected with $\\{Y=0\\}$, become one dimensional curves in
$\\{Y=0\\}$ which intersect at $\gamma\left(I\right)$. Let us denote these
curves as $w_{I}^{\mathrm{u}}\left(u\right)$ and
$w_{I}^{\mathrm{s}}\left(s\right)$, and work under a convention that
$w_{I}^{\mathrm{u}}\left(0\right)=w_{I}^{\mathrm{s}}\left(0\right)=\gamma\left(I\right)$.
(This can always be ensured by re-parameterising the curves.) We have added
the subscript $I$ to emphasize the dependence of the curves on the choice of
the energy level: on different energy levels we have a different Lyapunov
orbits, that lead to different curves. We have shown during the proof of Lemma
41 that the tangent vectors to these curves span the $X,P_{X}$ plane, i.e.

$\mathrm{span}\left(\pi_{X,P_{X}}\frac{d}{du}w_{I}^{\mathrm{u}}(u)|_{u=0},\pi_{X,P_{X}}\frac{d}{ds}w_{I}^{\mathrm{s}}(s)|_{s=0}\right)=\mathbb{R}^{2}.$
(71)

We shall take

$\Gamma:=\left\\{\left(\gamma\left(I\right),\theta\right):\text{
}I=H_{0}\left(L_{\mathcal{X}}\right)\text{, }\mathcal{X}\text{ is from
(\ref{eq:X-interval}) and }\theta\in\mathbb{T}\right\\}\subset\tilde{\Sigma}$
(72)

and prove that it is a well defined homoclinic channel by showing (5–8).

It will be convenient for us to check the transversality conditions (5–8) in
coordinates $X,P_{X},I,\theta$. In these coordinates we can parameterise
$\Lambda_{0}$ by $I,\theta$ (see Figure 3)

$\Lambda_{0}=\left\\{\left(\pi_{X,P_{X}}x_{0}\left(\mathcal{X}\right),I,\theta\right):I\in\mathbb{R},\text{
}\theta\in\mathbb{T},\,H_{0}\left(L_{\mathcal{X}}\right)=I\right\\}.$

Close the intersection of the stable and unstable manifold at
$\left(\gamma\left(I\right),\theta\right)$, we can parameterise the manifold
$W_{\Lambda_{0}}^{u}$ by $I,\theta,u$ as follows

$W_{\Lambda_{0}}^{u}=\left\\{\left(\pi_{X,P_{X}}w_{I}^{\mathrm{u}}\left(u\right),I,\theta\right):I\in\mathbb{R},\text{
}\theta\in\mathbb{T},\,u\in\mathbb{R}\right\\}.$

We can similarly parameterise the manifold $W_{\Lambda_{0}}^{u}$ by
$I,\theta,s$

$W_{\Lambda_{0}}^{s}=\left\\{\left(\pi_{X,P_{X}}w_{I}^{\mathrm{s}}\left(s\right),I,\theta\right):I\in\mathbb{R},\text{
}\theta\in\mathbb{T},\,s\in\mathbb{R}\right\\},$

and parameterise $\Gamma$ by $I,\theta$ as in (72). We see that at a point
$\left(\gamma\left(I\right),\theta\right)\in\Gamma$

$\displaystyle
T_{\left(\gamma\left(I\right),\theta\right)}W_{\Lambda_{0}}^{u}$
$\displaystyle=\mathrm{span}\left\\{\left(\pi_{X,P_{X}}\frac{d}{du}w_{I}^{\mathrm{u}}\left(u\right)|_{u=0},0,0\right),\left(\pi_{X,P_{X}}\frac{d}{dI}w_{I}^{\mathrm{u}}\left(0\right),1,0\right),\left(0,0,0,1\right)\right\\},$
(73) $\displaystyle
T_{\left(\gamma\left(I\right),\theta\right)}W_{\Lambda_{0}}^{s}$
$\displaystyle=\mathrm{span}\left\\{\left(\pi_{X,P_{X}}\frac{d}{ds}w_{I}^{\mathrm{s}}\left(s\right)|_{s=0},0,0\right),\left(\pi_{X,P_{X}}\frac{d}{dI}w_{I}^{\mathrm{s}}\left(0\right),1,0\right),\left(0,0,0,1\right)\right\\},$
(74) $\displaystyle T_{\left(\gamma\left(I\right),\theta\right)}\Gamma$
$\displaystyle=\mathrm{span}\left\\{\left(\pi_{X,P_{X}}\frac{d}{dI}\gamma\left(I\right),1,0\right),\left(0,0,0,1\right)\right\\}.$
(75)

We know that $\gamma\left(I\right)$ results from the intersection of
$w_{I}^{\mathrm{u}}\left(u\right)$ and $w_{I}^{\mathrm{s}}\left(s\right)$ at
$u=s=0$

$\gamma\left(I\right)=w_{I}^{\mathrm{u}}\left(0\right)=w_{I}^{\mathrm{s}}\left(0\right),$

so

$\pi_{X,P_{X}}\frac{d}{dI}\gamma\left(I\right)=\pi_{X,P_{X}}\frac{d}{dI}w_{I}^{\mathrm{u}}\left(0\right)=\pi_{X,P_{X}}\frac{d}{dI}w_{I}^{\mathrm{s}}\left(0\right).$
(76)

We now see that from (71), (73–76) we have (5–6).

We now turn to proving (7–8). For a fixed $I$ let us take
$x=x_{0}\left(\mathcal{X}\right)$ where $\mathcal{X}$ is such that
$H_{0}\left(L_{\mathcal{X}}\right)=I$. Let us also consider a fixed
$\lambda\in\mathbb{T}$. We see that the unstable and stable fibres of the
point $(x,\lambda)\in\Lambda_{0}$ are parameterised by $u$ and $s$,
respectively, as

$\displaystyle W_{\left(x,\lambda\right)}^{u}$
$\displaystyle=\left\\{\left(\pi_{X,P_{X}}w_{I}^{\mathrm{u}}\left(u\right),I,\theta_{I,\lambda}^{\mathrm{u}}\left(u\right)\right):u\in\mathbb{R}\right\\},$
$\displaystyle W_{\left(x,\lambda\right)}^{s}$
$\displaystyle=\left\\{\left(\pi_{X,P_{X}}w_{I}^{\mathrm{s}}\left(s\right),I,\theta_{I,\lambda}^{\mathrm{s}}\left(s\right)\right):s\in\mathbb{R}\right\\},$

where
$\theta_{I,\lambda}^{\mathrm{u}},\theta_{I,\lambda}^{\mathrm{s}}:\mathbb{R}\rightarrow\mathbb{T}$
are some functions, which parameterise the fibers along the angle coordinate.
Hence for every $I,\theta$

$\displaystyle
T_{\left(\gamma\left(I\right),\theta\right)}W_{\left(x,\theta\right)}^{u}$
$\displaystyle=\mathrm{span}\left\\{\left(\pi_{X,P_{X}}\frac{d}{du}w_{I}^{\mathrm{u}}\left(u\right)|_{u=0},0,\frac{d}{du}\theta_{I,\theta}^{\mathrm{u}}\left(u\right)|_{u=0}\right)\right\\},$
(77) $\displaystyle
T_{\left(\gamma\left(I\right),\theta\right)}W_{\left(x,\theta\right)}^{s}$
$\displaystyle=\mathrm{span}\left\\{\left(\pi_{X,P_{X}}\frac{d}{ds}w_{I}^{\mathrm{s}}\left(s\right)|_{s=0},0,\frac{d}{ds}\theta_{I,\theta}^{\mathrm{s}}\left(s\right)|_{s=0}\right)\right\\}.$
(78)

From (74), (75), (76) and (78) we obtain (7). Similarly, from (73), (75), (76)
and (77) we have (8). We have thus shown that $\Gamma$ is a well defined
homoclinic channel.

We now discuss the scattering map (70) associated with $\Gamma$. For fibers
$W_{\left(x_{-},\theta_{-}\right)}^{u}$ and
$W_{\left(x_{+},\theta_{+}\right)}^{s}$ to intersect we must have
$x_{+}=x_{-}$, since if this was not the case then the points $x_{-},x_{+}$
would lie on different energy levels, and their fibres would not meet. This
means that

$\sigma\left(x,\lambda\right)=\left(x,\pi_{\theta}\sigma\left(x,\lambda\right)\right),$

We now show how to obtain estimates for
$\pi_{\theta}\sigma\left(x,\lambda\right)$.

For every $\left(x,\lambda\right)\in\Lambda_{0}$ we have the homoclinic orbit
established in section 4.4 (see in particular Table 2) with the initial point
$x_{0}=x_{0}\left(x,\lambda\right)$ lying on the unstable fiber established in
Lemma 35. From (44) it follows that

$\left|\pi_{\theta}x_{0}\left(x,\lambda\right)-\lambda\right|\leq Mr=3\cdot
3\cdot 10^{-8}.$

As the point $x_{0}\left(x,\lambda\right)$ is iterated by $f_{0}$ the angle
$\theta$ changes. From Remark 40 we know that after the full excursion along
the homoclinic from Table 2 we return to the neighbourhood of $\Lambda_{0}$ at
an angle $\pi_{\theta}x_{0}\left(x,\lambda\right)+42\cdot h$, where $h=h(x)$
is from the interval (47). Such homoclinic excursion takes five iterates of
$f_{0}$ (see Table 2) so

$\pi_{\theta}f_{0}^{5}\left(x_{0}\left(x,\lambda\right)\right)=\pi_{\theta}x_{0}\left(x,\lambda\right)+42\cdot
h.$

We know that $f_{0}^{5}\left(x_{0}\left(x,\lambda\right)\right)$ lies in the
stable fiber of $f_{0}^{5}\left(\sigma\left(x,\lambda\right)\right)$. We also
know that

$f_{0}^{5}\left(\sigma\left(x,\lambda\right)\right)=f_{0}^{5}\left(x,\pi_{\theta}\sigma\left(x,\lambda\right)\right)=\left(x,\pi_{\theta}\sigma\left(x,\lambda\right)+5T(x)\right),$

where $T\left(x\right)$ is the period of the Lyapunov orbit. From Lemma 35 and
the $\mathcal{R}$-symmetry of the system since
$f_{0}^{5}\left(x_{0}\left(x,\lambda\right)\right)\in
W_{f_{0}^{5}\left(\sigma\left(x,\lambda\right)\right)}^{s}$ we have

$\left|\pi_{\theta}\left(f_{0}^{5}\left(x_{0}\left(x,\lambda\right)\right)-f_{0}^{5}\left(\sigma\left(x,\lambda\right)\right)\right)\right|\leq
Mr.$

This allows us to obtain the following estimate for the scattering map

$\displaystyle\pi_{\theta}\sigma\left(x,\lambda\right)$
$\displaystyle=\pi_{\theta}\sigma\left(x,\lambda\right)+5T(x)-5T(x)$
$\displaystyle=\pi_{\theta}f_{0}^{5}\left(\sigma\left(x,\lambda\right)\right)-5T(x)$
$\displaystyle=\pi_{\theta}f_{0}^{5}\left(x_{0}\left(x,\lambda\right)\right)+\pi_{\theta}\left(f_{0}^{5}\left(\sigma\left(x,\lambda\right)\right)-f_{0}^{5}\left(x_{0}\left(x,\lambda\right)\right)\right)-5T(x)$
$\displaystyle\in\pi_{\theta}x_{0}\left(x,\lambda\right)+42\cdot
h+\left[-Mr,Mr\right]-5T(x)$
$\displaystyle=\lambda+\left(\pi_{\theta}x_{0}\left(x,\lambda\right)-\lambda\right)+42\cdot
h+\left[-Mr,Mr\right]-5T(x)$
$\displaystyle\in\lambda+\left[-Mr,Mr\right]+42\cdot
h+\left[-Mr,Mr\right]-5T(x)$
$\displaystyle\in\lambda+[-0.82506903038,-0.82506533656].$ (79)

We have thus obtained our bound for the scattering map of the unperturbed
system.

We have finished our discussion about the unperturbed system. Now we turn to
showing that for sufficiently small $\varepsilon>0$ the manifold $\Lambda_{0}$
persists.

Recall the notation $\tilde{\Lambda}_{L}$ from (53), which stands for our
family of Lyapunov orbits in the extended phase space. Recall also that
$\Lambda_{0}=\tilde{\Lambda}_{L}\cap\\{Y=0\\}$. By Corollary 44 we know that
for sufficiently small $\varepsilon$ the normally hyperbolic invariant
manifold $\tilde{\Lambda}_{L}$ is perturbed to a nearby manifold
$\tilde{\Lambda}_{L}^{\varepsilon}$. We thus obtain
$\Lambda_{\varepsilon}:=\tilde{\Lambda}_{L}^{\varepsilon}\cap\\{Y=0\\}$ as the
perturbation of $\Lambda_{0}$. Since by Theorem 42 we know that
$\tilde{\Lambda}_{L}^{\varepsilon}$ is a manifold with a boundary that
consists of two, two dimensional invariant KAM tori. These two tori
intersected with $\\{Y=0\\}$ produce two curves, which are invariant under
$f_{\varepsilon}$. They become boundaries of $\Lambda_{\varepsilon}$. Thus
$\Lambda_{\varepsilon}$ is a normally hyperbolic invariant manifold, with
boundary, for $f_{\varepsilon}$. Similarly, all the two dimensional KAM tori
in $\tilde{\Lambda}^{\varepsilon}_{L}$ become one dimensional invariant tori
in $\Lambda_{\varepsilon}$.

We now turn to validating the assumptions of Theorem 11 to obtain our result.
This will be done in the local coordinates given by the affine change
involving $\tilde{A}_{0}$ and $x^{\ast}$, in which our $f_{\varepsilon}$ is
expressed. Recall that in Lemma 30 and Corollary 31 we have established bounds
on the intersection of the local unstable manifold of a Lyapunov orbit with
$\\{Y=0\\}$. This bound is valid in the following neighbourhood of
$\Lambda_{0}$

$U=\left[-r,r\right]^{4}\times\mathbb{T},$ (80)

where $r=3\cdot 10^{-8}$ is the constant from Lemma 30. We have performed a
computer assisted validation that the constant $L_{g}$ from (14) is

$L_{g}=90943.$ (81)

This value was computed using the method described in detail in C. The value
obtained in (81) is a large overestimate. By performing more careful checks,
for instance by subdividing $U$ into small fragments, the bound can be
significantly improved. Due to the small size of the neighbourhood (80) in
which we consider the local unstable manifold, we see that the constant $C$
from (15) can is very small:

$C=r=3\cdot 10^{-8}.$ (82)

Thus the large value of $L_{g}$ is not a problem for us, since $L_{g}$ enters
condition (17) multiplied by $C$. We use the bound $\bar{m}$ from Remark 34 to
compute

$\lambda=\left(\bar{m}\right)^{-30}=\frac{1}{1525.16}.$

The power $30$ comes from the fact that to complete a full turn round a
Lyapunov orbit involves $30$ local maps (29); see Table 1. We consider the
following strips666The positioning of the strips was motivated by computing
$\sum_{i=0}^{4}\pi_{I}g\left(0,f_{0}^{i}\left(x\right)\right)$, which
represents the change of energy after the homoclinic excursion along the
homoclinic. These five terms in the sum play the leading role in (83). In
Figure 4 we provide a plot of computer assisted bounds on
$\sum_{i=0}^{4}\pi_{I}g\left(0,f_{0}^{i}\left(x\right)\right)$, for different
choices of $\pi_{\theta}z$, and place our strips $S^{+},S^{-}$ for reference
in the figure.

$\displaystyle S^{-}$
$\displaystyle=\Lambda_{0}\cap\left\\{\theta\in\left[0.65-0.125,0.65+0.125\right]\right\\},$
$\displaystyle S^{+}$
$\displaystyle=\Lambda_{0}\cap\left\\{\theta\in\left[\pi+0.65-0.125,\pi+0.65+0.125\right]\right\\}.$

For each point $z\in S^{+}$ we compute $m=m\left(z\right)$ such that condition
(16) is fulfilled. Depending on the choice of $z$ the resulting $m$ can
differ. To compute it we used the bound on $\sigma$ from (79) and the fact
that the inner dynamics is given by (69) with the bound on $T(\mathcal{X})$
from (27). We check the assumptions of Theorem 9 by subdividing $S^{+}$ into
$25$ fragments along the $\theta$ coordinate, and validated the assumptions
for each fragment independently. For the first three fragments for which
$\theta$ was closest to $\pi+0.65-0.125$ it turned out that a good choice is
$m(z)=21$; for three fragments with $\theta$ close to $\pi+0.65+0.125$ we used
$m(z)=25$; and for the remaining we took $m(z)=23$. For $z\in S^{+}$ the point
$x\in W_{z}^{u}$ from condition 2. from Theorem 9, is taken as the first point
from homoclinic orbit from Lemma 37. For the alignment of $x$ along $\theta$
we use the estimate (44) from Lemma 35. We then validate that

$\sum_{i=0}^{m\left(z\right)-1}\pi_{I}g\left(0,f_{0}^{i}\left(x\right)\right)-\frac{1+\lambda}{1-\lambda}L_{g}C>0.$
(83)

The inequality (83) is validated with the aid of computer assisted estimates.

Figure 4: The computer assisted bound for the first five terms in the sum in
(83), depending on $\pi_{\theta}z$.

In the similar fashion we validate the assumptions of Theorem 10.

Conditions 1. and 2. from Theorem 11 are simple to validate since for every
point $z=\left(x_{0}\left(\mathcal{X}\right),\lambda\right)\in\Lambda_{0}$ we
have

$\pi_{\theta}f\left(z\right)=\lambda+T\left(\mathcal{X}\right),$

where $T(\mathcal{X})$ is the period of the Lyapunov orbit, whose bound is
known to us in (27).

###### Remark 46

Due to the fact that we work with the particular family of Lyapunov orbits,
which correspond to the energy of the comet Oterma, it turned out that the $m$
used in conditions 1. and 2. from Theorems 9 and 10 was a large number. This
resulted in the need of large number of iterates of $f_{0}$ when computing
(83), which meant long integration time. Also, we needed relatively wide
strips $S^{+},S^{-}$. This meant that we needed many subdivisions of the
strips to perform our validation. The long integration time and the large
number of sets increased the computational time of our proof. By making a more
careful choice of the energy level, one could focus on Lyapunov orbits for
which $m$ would be smaller. We have chosen not to do so to demonstrate that
the method is applicable for the choice of energy dictated by a concrete
physical object.

##  Appendix A Proof of Lemma 27

Proof. Take $v\in Q^{+}\left(0\right)\cap B$. Since $a_{i}\in\left(0,1\right)$
we see that
$\left|v_{1}\right|\geq\left|v_{i}\right|/a_{i-1}>\left|v_{i}\right|$, for
$i=2,3,4$. Since $f\left(0\right)=0$ from the mean value theorem we obtain

$\displaystyle\left\|f\left(v\right)\right\|_{\max}$
$\displaystyle=\left\|f\left(v\right)-f\left(0\right)\right\|_{\max}$
$\displaystyle\in\left\|\left[Df\left(B\right)\right]v\right\|_{\max}$
$\displaystyle=\left\|\left(A_{11}v_{1}+A_{12}\left(v_{2},v_{3},v_{4}\right),A_{21}v_{1}+A_{22}\left(v_{2},v_{3},v_{4}\right)\right)\right\|_{\max}$
$\displaystyle\geq\left|A_{11}v_{1}+A_{12}\left(v_{2},v_{3},v_{4}\right)\right|$
$\displaystyle\geq\left(c-\left\|A_{12}\right\|_{\max}\max\left(a_{1},a_{2},a_{3}\right)\right)\left|v_{1}\right|$
$\displaystyle=\left(c-\left\|A_{12}\right\|_{\max}\max\left(a_{1},a_{2},a_{3}\right)\right)\left\|v\right\|_{\max},$

as required.

##  Appendix B Proof of Lemma 28

Proof. Let $v=\left(v_{1},v_{2}\right)$ where $v_{1}=\pi_{1}v$ and
$v_{2}=\pi_{2,3,4}v$. Since $v\in Q^{+}\left(0\right)$ we see that

$\left|v_{1}\right|\geq a^{-1}\left\|v_{2}\right\|_{\max}.$

From the mean value theorem we obtain

$\displaystyle\left\|f\left(v\right)\right\|_{\max}$
$\displaystyle\in\left\|\left[Df\left(B\right)\right]v\right\|_{\max}$
$\displaystyle=\left\|\left(A_{11}v_{1}+A_{12}v_{2},A_{21}v_{1}+A_{22}v_{2}\right)\right\|_{\max}$
$\displaystyle\leq\max\left(\left|A_{11}\right|\left|v_{1}\right|+\left\|A_{12}\right\|_{\max}\left\|v_{2}\right\|_{\max},\left\|A_{21}\right\|_{\max}\left|v_{1}\right|+\left\|A_{22}\right\|_{\max}\left\|v_{2}\right\|_{\max}\right)$
$\displaystyle\leq\max\left(\left|A_{11}\right|+a\left\|A_{12}\right\|_{\max},\left\|A_{21}\right\|_{\max}+a\left\|A_{22}\right\|_{\max}\right)\left|v_{1}\right|$
$\displaystyle=\max\left(\left|A_{11}\right|+a\left\|A_{12}\right\|_{\max},\left\|A_{21}\right\|_{\max}+a\left\|A_{22}\right\|_{\max}\right)\left\|v\right\|_{\max}$

as required.

##  Appendix C Lipschitz bounds for the perturbation term

Here we give a method, with which we can check (14) in the case when $I(x)$ is
defined as (65). Below we start with Lemma 47, that can be applied to achieve
this. We have found though that in our particular case of the PER3BP, due to
long integration times, a direct application of Lemma 47 leads to
overestimates, which were too large for our needs. We therefore follow with
Lemma 48, which can be used for an inductive computation of $L_{g}$ by
expressing $f_{\varepsilon}$ as a composition of maps. This allowed us to
avoid long integration times and improved the estimate.

###### Lemma 47

Let $D^{2}I(x)$ stand for the Hessian of $I$ at $x$. Consider
$x_{0},x_{1}\in\mathbb{R}^{3}\times\mathbb{S}$ and let
$\mathbf{x}\subset\mathbb{R}^{3}\times\mathbb{S}$ be a convex set which
contains $x_{0}$ and $x_{1}$. Let $\mathbf{v}\subset\mathbb{R}^{4}$ be the
following interval enclosure

$\mathbf{v}^{\top}=\left[\left(\frac{\partial
f_{0}}{\partial\varepsilon}\left(\mathbf{x}\right)\right)^{\top}D^{2}I\left(f_{0}\left(\mathbf{x}\right)\right)Df_{0}\left(\mathbf{x}\right)+DI\left(f_{0}\left(\mathbf{x}\right)\right)\frac{\partial^{2}f_{0}}{\partial
x\partial\varepsilon}\left(\mathbf{x}\right)\right].$

If $\left\|\mathbf{v}\right\|\leq L_{g}$ then

$\left|\pi_{I}g\left(0,x_{1}\right)-\pi_{I}g\left(0,x_{0}\right)\right|\leq
L_{g}\left\|x_{1}-x_{0}\right\|.$

Proof. Let

$x_{s}:=x_{0}+s\left(x_{1}-x_{0}\right).$

From (68) we obtain

$\displaystyle\pi_{I}g\left(0,x_{1}\right)-\pi_{I}g\left(0,x_{0}\right)$
$\displaystyle=\int_{0}^{1}\frac{d}{ds}\nabla
I\left(f_{0}\left(x_{s}\right)\right)\cdot\frac{\partial
f_{0}}{\partial\varepsilon}\left(x_{s}\right)ds$
$\displaystyle=\int_{0}^{1}D^{2}I\left(f_{0}\left(x_{s}\right)\right)Df_{0}\left(x_{s}\right)\left(x_{1}-x_{0}\right)\cdot\frac{\partial
f_{0}}{\partial\varepsilon}\left(x_{s}\right)+\nabla
I\left(f_{0}\left(x_{s}\right)\right)\cdot\frac{\partial^{2}f_{0}}{\partial
x\partial\varepsilon}\left(x_{s}\right)\left(x_{1}-x_{0}\right)ds$
$\displaystyle=\int_{0}^{1}\frac{\partial
f_{0}}{\partial\varepsilon}\left(x_{s}\right)\cdot
D^{2}I\left(f_{0}\left(x_{s}\right)\right)Df_{0}\left(x_{s}\right)\left(x_{1}-x_{0}\right)+\nabla
I\left(f_{0}\left(x_{s}\right)\right)\cdot\frac{\partial^{2}f_{0}}{\partial
x\partial\varepsilon}\left(x_{s}\right)\left(x_{1}-x_{0}\right)ds$
$\displaystyle=\left(\int_{0}^{1}\left(\frac{\partial
f_{0}}{\partial\varepsilon}\left(x_{s}\right)\right)^{\top}D^{2}I\left(f_{0}\left(x_{s}\right)\right)Df_{0}\left(x_{s}\right)+DI\left(f_{0}\left(x_{s}\right)\right)\frac{\partial^{2}f_{0}}{\partial
x\partial\varepsilon}\left(x_{s}\right)ds\right)\left(x_{1}-x_{0}\right)$
$\displaystyle\in\mathbf{v}\cdot\left(x_{1}-x_{0}\right),$

so

$\left|\pi_{I}g\left(0,x_{1}\right)-\pi_{I}g\left(0,x_{0}\right)\right|\leq\left\|\mathbf{v}\right\|\left\|x_{1}-x_{0}\right\|\leq
L_{g}\left\|x_{1}-x_{0}\right\|,$

as required.

We now consider the case where $f$ is a composition of two functions
$f=f_{\varepsilon}^{1}\circ f_{\varepsilon}^{2}$. Our objective will be to
compute Lipschitz bound in terms of $x$ for

$g\left(\varepsilon,x\right):=\frac{1}{\varepsilon}\left(I\circ
f_{\varepsilon}^{2}\circ
f_{\varepsilon}^{1}\left(\varepsilon,x\right)-I(x)\right).$ (84)

The following lemma gives the bound for $\pi_{I}g\left(0,x\right)$ from bounds
for $f_{\varepsilon}^{1}$ and $f_{\varepsilon}^{2}.$

###### Lemma 48

Assume that

$\displaystyle f_{\varepsilon}$ $\displaystyle=f_{\varepsilon}^{2}\circ
f_{\varepsilon}^{1},$ $\displaystyle g\left(\varepsilon,x\right)$
$\displaystyle=\frac{1}{\varepsilon}\left(f_{\varepsilon}\left(x\right)-f_{0}\left(x\right)\right),$
$\displaystyle g_{1}\left(\varepsilon,x\right)$
$\displaystyle=\frac{1}{\varepsilon}\left(f_{\varepsilon}^{1}\left(x\right)-f_{0}^{1}\left(x\right)\right),$
$\displaystyle g_{2}\left(\varepsilon,x\right)$
$\displaystyle=\frac{1}{\varepsilon}\left(f_{\varepsilon}^{2}\left(x\right)-f_{0}^{2}\left(x\right)\right),$

and that

$\displaystyle\left\|\pi_{I}\left(g_{1}\left(0,x_{1}\right)-g_{1}\left(0,x_{0}\right)\right)\right\|$
$\displaystyle\leq L_{g}^{1}\left\|x_{1}-x_{0}\right\|,$
$\displaystyle\left\|\pi_{I}\left(g_{2}\left(0,x_{1}\right)-g_{2}\left(0,x_{0}\right)\right)\right\|$
$\displaystyle\leq L_{g}^{2}\left\|x_{1}-x_{0}\right\|,$
$\displaystyle\left\|f_{0}^{1}\left(x_{1}\right)-f_{0}^{1}\left(x_{0}\right)\right\|$
$\displaystyle\leq L_{f}^{1}\left\|x_{1}-x_{0}\right\|.$

Then

$\left\|\pi_{I}\left(g\left(0,x_{1}\right)-g\left(0,x_{0}\right)\right)\right\|\leq\left(L_{g}^{2}L_{f}^{1}+L_{g}^{1}\right)\left\|x_{1}-x_{0}\right\|.$

Proof. Consider fixed $x_{1},x_{0}$. Then

$\displaystyle\left\|\pi_{I}g\left(\varepsilon,x_{1}\right)-\pi_{I}g\left(\varepsilon,x_{0}\right)\right\|$
$\displaystyle=\left\|\frac{1}{\varepsilon}\left(I\circ
f_{\varepsilon}^{2}\circ
f_{\varepsilon}^{1}\left(x_{1}\right)-I\left(x_{1}\right)-\left(I\circ
f_{\varepsilon}^{2}\circ
f_{\varepsilon}^{1}\left(x_{0}\right)-I\left(x_{0}\right)\right)\right)\right\|$
$\displaystyle\leq\left\|\frac{1}{\varepsilon}\left(I\circ
f_{\varepsilon}^{2}\circ f_{\varepsilon}^{1}\left(x_{1}\right)-I\circ
f_{\varepsilon}^{1}\left(x_{1}\right)-\left(I\circ f_{\varepsilon}^{2}\circ
f_{\varepsilon}^{1}\left(x_{0}\right)-I\circ
f_{\varepsilon}^{1}\left(x_{0}\right)\right)\right)\right\|$
$\displaystyle\quad+\left\|\frac{1}{\varepsilon}\left(I\circ
f_{\varepsilon}^{1}\left(x_{1}\right)-I\left(x_{1}\right)-\left(I\circ
f_{\varepsilon}^{1}\left(x_{0}\right)-I\left(x_{0}\right)\right)\right)\right\|$
$\displaystyle=\left\|\pi_{I}g_{2}\left(\varepsilon,f_{\varepsilon}^{1}\left(x_{1}\right)\right)-\pi_{I}g_{2}\left(\varepsilon,f_{\varepsilon}^{1}\left(x_{0}\right)\right)\right\|+\left\|\pi_{I}g_{1}\left(\varepsilon,x_{1}\right)-\pi_{I}g_{1}\left(\varepsilon,x_{0}\right)\right\|$
$\displaystyle\leq
L_{g}^{2}\left\|f_{\varepsilon}^{1}\left(x_{1}\right)+f_{\varepsilon}^{1}\left(x_{0}\right)\right\|+L_{g}^{1}\left\|x_{1}+x_{0}\right\|+o\left(\varepsilon\right)$
$\displaystyle\leq\left(L_{g}^{2}L_{f}^{1}+L_{g}^{1}\right)\left\|x_{1}+x_{0}\right\|+o\left(\varepsilon\right)$

and the result follows by passing with $\varepsilon$ to zero.

The bound (14) computed in the proof (see section 4.6) is for the Poincaré map
$f_{\varepsilon}$ defined in (64). This can be computed by considering

$f_{\varepsilon}=f_{\varepsilon,N}\circ\ldots\circ f_{\varepsilon,1},$ (85)

where $f_{\varepsilon,i}$ are local maps of the form

$\displaystyle f_{\varepsilon,i}\left(x\right)$
$\displaystyle=A_{i}^{-1}\left(\Phi_{t}^{\varepsilon}\left(x_{i-1}+A_{i-1}x\right)-x_{i}\right)\qquad\text{for
}i=1,\ldots,N-1,$ $\displaystyle f_{\varepsilon,N}\left(x\right)$
$\displaystyle=A_{0}^{-1}\left(\mathcal{P}^{\varepsilon}\left(x_{N-1}+A_{N-1}x\right)-x_{0}\right).$

In the inductive application of Lemma 48 we can use the bound $\bar{m}$ from
Remark 34 and take $L_{f}^{1}=\bar{m}^{k}$ at the $k$-th inductive step.

The energy after an iterate of a local map is

$I\left(f_{\varepsilon,i}\left(x\right)\right)=H_{0}\left(\Phi_{t}^{\varepsilon}\left(x_{i-1}+A_{i-1}x\right)\right).$

Since

$\displaystyle\frac{d}{dx}\frac{d}{d\varepsilon}H_{0}\left(\Phi_{t}^{\varepsilon}\left(x_{i-1}+A_{i-1}x\right)\right)$
$\displaystyle=\frac{d}{dx}\left(\nabla
H_{0}\left(\Phi_{t}^{\varepsilon}\left(x_{i-1}+A_{i-1}x\right)\right)\cdot\frac{\partial\Phi_{t}^{\varepsilon}}{\partial\varepsilon}\left(x_{i-1}+A_{i-1}x\right)\right)$
$\displaystyle=\left(\frac{\partial\Phi_{t}^{\varepsilon}}{\partial\varepsilon}\left(x_{i-1}+A_{i-1}x\right)\right)^{\top}D^{2}H_{0}\left(\Phi_{t}^{\varepsilon}\left(x_{i-1}+A_{i-1}x\right)\right)\frac{\partial\Phi_{t}^{\varepsilon}}{\partial
x}\left(x_{i-1}+A_{i-1}x\right)A_{i-1}$
$\displaystyle\quad+DH_{0}\left(\Phi_{t}^{\varepsilon}\left(x_{i-1}+A_{i-1}x\right)\right)\frac{\partial^{2}\Phi_{t}^{\varepsilon}}{\partial
x\partial\varepsilon}\left(x_{i-1}+A_{i-1}x\right)A_{i-1}$

the $\mathbf{v}^{\top}$ used for the computation of $L_{g}$ for a local map
$f_{\varepsilon,i}$ by means of Lemma 47 is computed as

$\mathbf{v}^{\top}=\left[\left(\frac{\partial\Phi_{t}^{\varepsilon}}{\partial\varepsilon}\left(\mathbf{x}\right)\right)^{\top}D^{2}H_{0}\left(\Phi_{t}^{\varepsilon}\left(\mathbf{x}\right)\right)\frac{\partial\Phi_{t}^{\varepsilon}}{\partial
x}\left(\mathbf{x}\right)A_{i-1}+DH_{0}\left(\Phi_{t}^{\varepsilon}\left(\mathbf{x}\right)\right)\frac{\partial^{2}\Phi_{t}^{\varepsilon}}{\partial
x\partial\varepsilon}\left(\mathbf{x}\right)A_{i-1}\right].$ (86)

(For the last local map $f_{\varepsilon,N}$ we take
$\mathcal{P}^{\varepsilon}$ instead of $\Phi_{t}^{\varepsilon}$ in (86).)

## References

  * [1] V. I. Arnol’d, Instability of dynamical systems with many degrees of freedom, Dokl. Akad. Nauk SSSR 156 (1964) 9–12.
  * [2] A. Delshams, R. de la Llave, T. M. Seara, Geometric properties of the scattering map of a normally hyperbolic invariant manifold, Adv. Math. 217 (3) (2008) 1096–1153.
URL https://doi-org.ezproxy.fau.edu/10.1016/j.aim.2007.08.014

  * [3] A. Delshams, R. Ramírez-Ros, Poincaré-Mel’nikov-Arnol’d method for analytic planar maps, Nonlinearity 9 (1) (1996) 1–26.  doi:10.1088/0951-7715/9/1/001.
URL https://doi.org/10.1088/0951-7715/9/1/001

  * [4] M. Gidea, R. de la Llave, T. M-Seara, A general mechanism of diffusion in Hamiltonian systems: qualitative results, Comm. Pure Appl. Math. 73 (1) (2020) 150–209.  doi:10.1002/cpa.21856.
URL https://doi.org/10.1002/cpa.21856

  * [5] M. Gidea, R. de la Llave, T. M. Seara, A general mechanism of instability in Hamiltonian systems: skipping along a normally hyperbolic invariant manifold, Discrete Contin. Dyn. Syst. 40 (12) (2020) 6795–6813.  doi:10.3934/dcds.2020166.
URL https://doi.org/10.3934/dcds.2020166

  * [6] M. J. Capiński, J. Gonzalez, J.-P. Marco, J. D. Mireles James, Computer assisted proof of drift orbits along normally hyperbolic manifolds, https://arxiv.org/abs/2102.06436.
URL https://arxiv.org/abs/2102.06436

  * [7] A. Delshams, V. Kaloshin, A. de la Rosa, T. M. Seara, Global instability in the restricted planar elliptic three body problem, Comm. Math. Phys. 366 (3) (2019) 1173–1228.  doi:10.1007/s00220-018-3248-z.
URL https://doi.org/10.1007/s00220-018-3248-z

  * [8] M. J. Capiński, M. Gidea, Arnold diffusion, quantitative estimates and stochastic behavior in the three-body problem (2020).
URL https://arxiv.org/abs/1812.03665

  * [9] M. J. Capiński, P. Zgliczyński, Transition tori in the planar restricted elliptic three-body problem, Nonlinearity 24 (5) (2011) 1395–1432.  doi:10.1088/0951-7715/24/5/002.
URL https://doi.org/10.1088/0951-7715/24/5/002

  * [10] M. J. Capiński, M. Gidea, R. de la Llave, Arnold diffusion in the planar elliptic restricted three-body problem: mechanism and numerical verification, Nonlinearity 30 (1) (2017) 329–360.  doi:10.1088/1361-6544/30/1/329.
URL https://doi.org/10.1088/1361-6544/30/1/329

  * [11] T. Kapela, M. Mrozek, D. Wilczak, P. Zgliczyński, Capd::dynsys: a flexible c++ toolbox for rigorous numerical analysis of dynamical systems, Communications in Nonlinear Science and Numerical Simulation (2020) 105578doi:https://doi.org/10.1016/j.cnsns.2020.105578.
URL https://www.sciencedirect.com/science/article/pii/S1007570420304081

  * [12] G. Alefeld, Inclusion methods for systems of nonlinear equations—the interval Newton method and modifications, in: Topics in validated computations (Oldenburg, 1993), Vol. 5 of Stud. Comput. Math., North-Holland, Amsterdam, 1994, pp. 7–26.
  * [13] M. W. Hirsch, C. C. Pugh, M. Shub, Invariant manifolds, Bull. Amer. Math. Soc. 76 (1970) 1015–1019.  doi:10.1090/S0002-9904-1970-12537-X.
URL https://doi.org/10.1090/S0002-9904-1970-12537-X

  * [14] V. Szebehely, Theory of Orbits: The Restricted Problem of Three Bodies, Academic Press, 1967.
  * [15] W. S. Koon, M. W. Lo, J. E. Marsden, S. D. Ross, Heteroclinic connections between periodic orbits and resonance transitions in celestial mechanics, Chaos 10 (2) (2000) 427–469.  doi:10.1063/1.166509.
URL https://doi.org/10.1063/1.166509

  * [16] M. J. Capiński, Computer assisted existence proofs of Lyapunov orbits at $L_{2}$ and transversal intersections of invariant manifolds in the Jupiter-Sun PCR3BP, SIAM J. Appl. Dyn. Syst. 11 (4) (2012) 1723–1753.  doi:10.1137/110847366.
URL https://doi.org/10.1137/110847366