text:
\begin{document}

\title{A positive solution to Hilbert's 10th problem}

\author{Michael Pfender\footnote{michael.pfender@alumni.tu-berlin.de}}
\date{April 2014\\last revised \today}
\maketitle

\abstract{
Polynome codes and code evaluation; arithmetical theory frames;
$\mu$-recursive race for decision; decision correctness;
decision termination in Zermelo-Fraenkel set theory ZFC+ with
axiom of choice and consistency provability; decision correctness
in theory $\T = \PR$ of Primitive Recursion;
comparison with the negative result of Matiyasevich; positive solution
for each single diophantine polynomial in
\emph{p.\,r.\ non-infinite-descent theory} $\piR=\PR+(\pi).$
}

\tableofcontents

\section*{Introduction}

Within theory $\ZFC^+ = \ZFC+\Con_{\ZFC}$ of Zermelo-Fraenkel
\textbf{set} theory with axiom of choice $\AC,$ strengthened by
formula $\Con_{\ZFC}$ which is to express $\ZFC$'s \emph{internal,}
\emph{g\"odelised} consistency, we solve Hilbert's 10th
problem positively: we organise decision of
diophantine polynome codes---decision on overall
\emph{non-nullity}---as an enumerative $\mu$-recursive race for
a (first) zero (\emph{counterexample}), against race for a first
internal $\ZFC$-non-nullity \emph{proof} for a given such polynomial
code, given as the (nested) list of coefficients. Comparison with
Matiyasevich's negative solution of Hilbert's 10th problem gives
inconsistency of theory $\ZFC+\Con_{\ZFC}$ whence self-inconsistency
$\ZFC \derives \neg \Con_{\ZFC}.$

In a final section we
plug our positive solution of the problem into the constructive
framework of \emph{p.\,r.\ non-infinite descent theory} $\piR=\PR+(\pi)$
out of \emph{Arithmetical Foundations} in the References.

This is to give a decision algorithm for each single
diophantine equation (in a uniform way), as asked in the original
Hilbert's 10th problem.

\section{Hilbert's 10th Problem}

We attempt a positive solution to Hilbert's 10th problem. In its original
form it reads:

\bigskip
10. DETERMINATION OF THE SOLVABILITY OF A DIOPHANTINE EQUATION
Given a diophantine equation with any number of unknown quantities
and with rational integer numerical coefficients: \emph{To devise a process
according to which it can be determined by a finite number of operations
whether the equation is solvable in rational integers.}

[translation quoted from \NAME{Matiyasevich} 1993.]

\bigskip
Formally, this text allows for a separate decision algorithm
(``process'') for each diophantine polynomial. But it is clear that
a decision-\emph{family} must be \emph{uniform} in a suitable sense.

\emph{Correctness} of our alleged $\mu$-recursive decision algorithm
$\nabla_{\ZFC}: \mr{DIO} \parto \two=\set{0,1}$ builds, within $\ZFC^+,$
on diophantine soundness inferred by $\Con_{\ZFC}$ over $\ZFC.$
Termination follows from (countable) Choice. This already within
$\ZFC.$ Together this gives the wanted decision $\nabla=\nabla_{\ZFC}$
within $\ZFC^+,$ of all polynome codes in $\mr{DIO}\subset\N.$

Comparison with Matiyasevich's negative Theorem,
\emph{unsolving} Hilbert's 10th Problem,
theorem in particular of (classically quantified
Arithmetical Theory) $\ZFC^+,$ gives a contradiction within
$\ZFC^+,$ hence {self-inconsistency} of $\ZFC,$
and from that in particular $\omega$-inconsistency.

In a final section we show correctness and irrefutable termination
of \emph{localised} decision $\nabla[D]$---for each single diophantine
polynomial $D=D(\vec{x})$---within the constructive framework of
p.\,r.\ \emph{finite-descent-theory} $\piR=\piR+\Con_{\piR}$ out
of op.\,cit.

\section{Polynome coding and code evaluation}

Diophantine polynomials $D = D(\vec{\xi}): \Z^* \to \Z$ (``in $\DIO$'')
are \LaTeX/\textbf{ASCII} coded into
  $$\mr{DIO} \defeq \Z^{\an{*}} \iso \union_{m\geq1} \Z[\xi_1,\ldots,\xi_m]
                 = \union_{m \geq 1} \Z[\xi_1][\xi_2] \ldots [\xi_m]$$
as nested coefficient \emph{lists} $\Z^{\an{*}} \subset \N.$

\smallskip
\noindent
$[\,$The symbols $\xi_i$ are the \emph{indeterminates.}$\,]$

\subsubsection*{Example:}
\begin{align*}
D = D(\xi_1,\xi_2)
& = (2\cdot\xi_1^{\,0}+3\cdot \xi_1^{\,1}- 4\cdot\xi_1^{\,3})\cdot\xi_2^{\,0} \\
& \quad
    +(0\cdot\xi_1^{\,0}+3\cdot\xi_1^{\,1}-7\cdot\xi_1^{\,2})\cdot\xi_2^{\,1}
                                           + (1-4\cdot\xi_1)\cdot\xi_2^{\,2}
\end{align*}
is coded 1-1 as (nested) {coefficient list}
\begin{align*}
\ccode{D}
& = \an{\an{2;3;0;4};\an{0;3;-7};\an{0};\an{1;-4}}: \\
& \one \to \mr{DIO} \bydefeq \Z^{\an{*}} \subset \N: \\
& \text{\emph{defined element, point} of $\mr{DIO}$}
\end{align*}

\subsubsection*{PR evaluation of $\mbf{DIO}$ codes:}

Evaluation $\mr{ev} = \mr{ev}(d,\vec{x}): \mr{DIO} \times \Z^*$
is PR {defined}
\begin{align*}
& \mr{ev}(d,\an{\vec{x};x_{m+1}\,}) = \mr{ev}(d,\an{x_1;\ldots;x_m;x_{m+1}\,})\\
& \defeq \mr{ev}(\horner(d,x_{m+1}),\an{\,\vec{x}\,}): \\
& \mr{DIO} \times \Z^* \supset \Z[\vec\xi,\xi_{m+1}] \times \Z^{m+1}
    \ovs{\iso} (\Z[\vec\xi] [\xi_{m+1}] \times (\Z^m \times \Z) \\
& \xto{\iso} (\Z[\vec\xi][\xi_{m+1}] \times \Z) \times \Z^m
              \ovs{\horner \times \id} \Z[\vec\xi] \times \Z^m \xto{\ev} \Z,
\end{align*}
recursively by iterative application of Horner's schema
to the hitherto trailing argument, until all of the arguments
(constants or variables) are substituted into
their corresponding indeterminates $\xi_j.$

\smallskip
Result then is the integer $\mr{ev}(d,\vec{x}),$
constant or integer variable.

\smallskip
For the \textbf{example} above, $D = D(\xi_1,\xi_2),$ with argument string
$\an{x_1;x_2} :\,= \an{23;64} \in \Z^*,$ we get
\begin{align*}
& \mr{ev}(d,\an{x_1;x_2})
    = \mr{ev}(\an{\an{2;3;0;4};\an{0;3;-7};\an{0};\an{1;-4}},\an{23;64}) \\

& = \horner(\,
    ((((-4\cdot64+1)\cdot\xi_1+0))
                       \cdot\xi_1 +(-7\cdot64+3)\cdot64)\cdot\xi_1 \\
& \qquad\qquad\qquad\qquad\qquad\qquad\quad
    + ((4\cdot64+0)\cdot64+3)\cdot64+2\,,23\,) \\
& = ((((-4\cdot64+1)\cdot23+0))\cdot23 +(-7\cdot64+3)\cdot64)\cdot23 \\
& \qquad\qquad\qquad\qquad\qquad\qquad
    + ((4\cdot64+0)\cdot64+3)\cdot64+2
\end{align*}
\textbf{First step:} apply Horner's schema to coefficient list
$d \in \mr{DIO}$ und (trailing) Argument $x_2:$ {indeterminate}
$\xi_1$ is {coded} by {list nesting} and is seen
as a \emph{constant}, as an element of intermediate ring $\Z[\xi_1]:$
  $$\Z[\xi_1,\xi_2] \bydefeq \Z[\xi_1][\xi_2] \bydefeq (\Z[\xi_1])\,[\xi_2].$$
\textbf{Last---here second---step:} evaluation of
$\Z[\xi_1]$ polynomial in remaining indeterminate $\xi_1$ on
remaining argument $x_1,$ by a last application of Horner's schema.

\section{Arithmetical frame theories}

We consider here as frame theories---for our decision algorithm --
\textbf{on one hand} classically quantified arithmetical theories
$\T = \bfQ+\AC$ with (countable) axiom of choice, as in particular
Zermelo-Fraenkel set theory $\T = \ZFC = \ZF+\AC.$ Frame then is the strengthening
  $$\T^+ = \T+\Con_\T = \ZFC+\Con_{\ZFC}$$
of $\T$ by its own consistency-\emph{formula}
\begin{eqnarray*}
\Con_\T &=& \neg\,(\exists\,k \in \N)\,\Pro_\T(k,\code{\false}) \\
&=& (\forall k)\,\neg\,\Pro_\T(k,\code{\false})\ (\text{G\"odel}),
\end{eqnarray*}
see \NAME{Smorynski} 1977 and op.\,cit.

Strengthening by this consistency formula will provide for
\emph{correctness} of our \emph{decision process} (Hilbert).

\smallskip
\textbf{On the other hand} we take as frame the
Free-Variables (categorical) theory $\T = \PR = \PRa$ of
\emph{Primitive Recursion with predicate abstraction into subsets}
$$(\chi = \chi(a): A \to \two)
  \boldsymbol\mapsto \set{A:\chi}=\set{a\in A:\chi(a)}$$
out of op.\,cit.\,,\
$\T = \bfS$
in Smorynski's notation, as well as \emph{descent theory}
$\piR = \piR^+ = \piR+\Con_{\piR}:$ that theory is self-consistent,
$\piR \derives \Con_{\piR,}$ main result of op.\,cit.

\section{A $\mu$-recursive race for decision}

We \textbf{define} an enumerative \emph{race}---for
$d \in \mr{DIO}$ thought \emph{passive, fixed,} and $k \in \N$
\emph{running}---for satisfaction of
\begin{align*}
& \ph_0(d,k) = [\,\ev(d,\ct_*k)=0\,]\ \text{against} \\
& \ph_1(d,k) = \Pro_{\T}(k,\code{(\vec{x})\ev(d,\vec{x})\neq 0}):
  \mr{DIO} \times \N \to \two = \set{0,1}, \\
& \ct_* = \ct_*\,k: \N \ovs{\iso} \Zlist
  \ \text{Cantor-type \emph{count},}\
    \vec{x} \in \Z^*\ \text{free under code.}
\end{align*}

This race towards \emph{termination} is defined as
a---formally partial---$\mu$-recursive mapping as follows
within the theory $\hatT$ of \emph{partial PR maps,} i.\,e.\ of
(partially defined) \emph{$\mu$-recursive maps,} cf.\ again op.\,cit.:
  $$t = t(d) = \mu\set{k\,|\,\ph_0(d,k)\,\lor\,\ph_1(d,k)}:
                                          \mr{DIO} \parto \N. \quad (*)$$

\bigskip
\textbf{Decision candidate} then is
\begin{align*}
\nabla d =
& \begin{cases}
  0\ \myif\ \ph_0(d,t(d)) \\
  1\ \myif\ \ph_1(d,t(d)) \\
  \end{cases} \\
& \quad \\
=
& \begin{cases}
  0\ \myif\ \ev(d,\ct_*(t(d)))=0 \\
        \qquad (\emph{zero found}) \\
  1\ \myif\ \Pro_{\T}(t(d),\code{\ev(d,\vec{x})\neq 0}) \\
        \qquad (\text{\emph{internal proof} found
                        for \emph{global non nullity}})
  \end{cases} \\
& : \mr{DIO} \overset{(\id,t)} {\parto} \mr{DIO} \times \N \to \two.
\end{align*}

\textbf{Question:} Is $\nabla$ \emph{well-defined} as a partial map?
In which frame?

\subsubsection*{Well-definedness of the decision within
$\T^+ = \ZFC^+=\ZFC+\Con_{\ZFC} = \T+\Con_\T:$}
\begin{align*}
\T^+\ \derives\
& \ph_0(d,k)\,\land\,\ph_1(d,k') \\
&     \qquad (\text{\emph{cases-overlap}
                        {$\mathit{Assumption}$}}) \\
& \implies \ev(d,\ct_*k)=0 \\
& \quad\,\land\,\Pro_\T(k',\code{(\vec{x})\,\ev(d,\vec{x})\neq 0}) \\
& \implies \Pro_\T(j(k,k'),\code{\false}) \\

& \implies \neg\,\Con_\T \implies \false, \\
& j=j(k,k'): \N^2\to \N\ \text{suitable.}
\end{align*}

\textbf{Consequence:}
  $$\T^+\ \derives\ \neg\,[\ph_0(d,k)\,\land\,\ph_1(d,k')\,]:
                                     \mr{DIO} \times \N^2 \to \two,$$
$\nabla = \nabla_{\T}(d): \mr{DIO} \parto \N$ is \emph{well-defined}
as a (\emph{formally partial}) $\mu$-recursive map, within
$\T^+ = \T+\Con_\T.$

\subsubsection*{Well-definedness of decision within \emph{descent}
theory $\piR:$}

We consider now \emph{descent theory} $\piR$ out of op.\,cit.\
strengthening $\PR$ by axiom $(\pi)$ of
\emph{non-infinite endo driven descending complexity with complexity
values in polynomial semiring $\N[\omega],$} and its logical properties,
in particular \emph{soundness} giving $\piR\derives \Con_{\piR}.$

Decision $\nabla=\nabla_{\piR}(d):\mr{DIO}\parto \two$
is in fact well-defined as a partial PR map, within theory $\piR,$
since---in parallel to the above case $\T = \ZFC:$
\begin{align*}
\piR\ \derives\
& \ph_0(d,k)\,\land\,\ph_1(d,k') \\
&     \qquad (\text{\emph{cases-overlap}
                        {$\mathit{Assumption}$}}) \\
& \implies \ev(d,\ct_*k)=0 \\
& \quad\,\land\,\Pro_{\piR}(k',\code{(\vec{x})\,\ev(d,\vec{x})\neq 0}) \\
& \implies \Pro_{\piR}(j(k,k'),\code{\false}) \\

& \implies \text{``$\neg\,\Con_{\piR}$''} \implies \false, \\
& j=j(k,k'): \N^2\to \N\ \text{suitable.}
\end{align*}
The latter since $\piR\derives \Con_{\piR}.$

\subsubsection*{Well-definedness of DIO-decision within $\PR$ itself}
Decision $\nabla=\nabla_{\PR}(d):\mr{DIO}\parto \two$
is well-defined as a partial PR map, within theory $\hatPRa$ of
partial PR maps since
\begin{align*}
\hatPRa\ \derives\
& \ph_0(d,k)\,\land\,\ph_1^{\DIO}(d,k') \\
&     \qquad (\text{\emph{cases-overlap}
                        {$\mathit{Assumption}$}}) \\
& \Iff \ev(d,\ct_*k)=0 \\
& \quad \land\,\Pro_{\DIO}(k',\code{(\vec{x})\,\ev(d,\vec{x})\neq 0}) \\
& \implies \Pro_{\DIO}(j(k,k'),\code{\false}) \\

& \implies \false, \\
& j=j(k,k'): \N^2\to \N\ \text{suitable.}
\end{align*}
The latter by \emph{diophantine soundness} of $\T=\PR,$ see
\NAME{Smorynski} 1977, \textsc{Theorem} \textbf{4.1.4}.

\section{Decision Correctness}

\textbf{Decision Correctness, result-0-case:}
\begin{align*}
\T \derives\
& [\,\ph_0(d,t(d))
    \implies \mr{ev}(d,\ct_*\,\circ\,t(d))=0\,] \\
& \subseteq\,\true_{\mr{DIO}}:
    \mr{DIO} \overset{(\id,t)} {\parto} \mr{DIO} \times \N \to \two:
\end{align*}
\textbf{If} race-for-decision $\nabla$ \emph{terminates} on
DIO-code $d,$ with \textbf{result} $0,$ \textbf{then}
(evaluation of) $d$ has (at least) one zero, namely
$$\ct_*\,\circ\,t(d) \in \N.$$

\subsubsection*{Correctness, result-1-case:}
\begin{align*}
\T \derives\
& \ph_1(d,k) \implies \Pro_{\DIO}(k,\code{\ev(d,\vec{x})\neq 0}) \\
& \implies \ev(d,\vec{x})\neq 0:
  (\mr{DIO} \times \N) \times \Z^* \to \two, \\
& (d \in \mr{DIO},\ k \in \N,\ \vec{x} \in \Z^*\ \text{all free}), \\
& \quad
    \text{or, with quantifier decoration:} \\
\T \derives\
& (\forall\,d \in \mr{DIO})(\forall\,k \in \N)
    (\forall\, \vec{x} \in \Z^*) \\
& [\,\ph_1^{\T} (d,k) \implies \Pro_{\DIO}(k,\code{\ev(d,\vec{x})\neq 0}) \\
& \implies \ev(d,\vec{x})\neq 0\,].
\end{align*}
\textbf{If} race-for-decision $\nabla$ \emph{terminates} on
DIO-code $d,$ with \textbf{result} $1,$ \textbf{then}
(evaluation of) $d$ has no zeroes.

This because of \emph{Diophantine Soundness} of $\T,$ see
\NAME{Smorynski} 1977, \textsc{Theorem} \textbf{4.1.4} again.

\subsubsection*{Correctness in result-1-case,
under termination condition:}
Substitution of $t(d)$ for $k$ in the above gives
\begin{align*}
\T^+,\piR,\PR\ \derives\
& [\,\ph_1^{\DIO} (d,t) \implies \ev(d,\vec{x})\neq 0\,]
  \subseteq\,\true_{\mr{DIO} \times \Z^*},\\
& d \in \mr{DIO},\ \vec{x} \in \Z^*\ \text{both free}:
\end{align*}
\emph{Correctness} of $\nabla(d)$ where defined, in \emph{both}
defined cases: in case of reaching
\textbf{result} 0, as well as in case of reaching
\textbf{result} 1.

\smallskip
[\,For {partial maps} $f,\,g: A \parto B,$
$f\,\subseteq\,g$ designates \emph{inclusion} of the
\emph{graphs} of $f$ and $g.$]

\section{Termination}

We show first

\smallskip
\textbf{Pointwise non-derivability of non-termination:}

For no diophantine \emph{point} $d_0: \one \to \mr{DIO}$
$\T$ derives non-termination of $t$ at $d_0.$

\textbf{Proof:}
\inference{}
{ {$\mr{Assumption}$} \\
& $\T \derives\ (\vec{x})\ev(d_0,\vec{x}) \neq 0$ \qquad $(\bullet)$ \\
& \qquad $\land\,(k)\,\neg\,\Pro_\T
                  (k,\code{(\vec{x})\,\mr{ev}(d_0,\vec{x}) \neq 0})$
}
{ $\T \derives\ \Pro_\T
     (\num\,\ulj,\code{(\vec{x})\,\mr{ev}(d_0,\vec{x}) \neq 0})$ \\
& $\land\,(k)\,\neg\,\Pro_\T
                (k,\code{(\vec{x})\,D(\vec{x}) \neq_Z 0})$
}
a contradiction:
appropriate $\ulj$ is available from $(\bullet)$ via
derivation-to-$\Proof$-internalisation (\emph{g\"odelisation}).

\bigskip
$[\,$For the time being we consider $\T$ as frame, not (yet)
$\T^+=\T+\Con_\T.\,]$

\bigskip
For $\T = \bfQ$ quantified, with (countable) \emph{axiom of choice}
$\ACC,$ in particular $\bfQ = \PA+\ACC$ Peano Arithmetic with choice,
we define the \emph{undecided part} of $\mr{DIO}$ as
\begin{eqnarray*}
\Psi &=& \Psi^{\bfQ} \\
&=& \set{d \in \mr{DIO}:\forall\,k\ \mr{ev}(d,\ct_*\,k) \neq 0 \\
&& \,\land\,\forall\,k
  \ \neg\,\Pro_{\bfQ}(k,\code{(\vec{x})\,\mr{ev}(d,\vec{x}) \neq 0)})} \\
&\subset& \mr{DIO} = \Zlist \subset \N.
\end{eqnarray*}
With this definition we get
\begin{align*}
\bfQ \derives\
& \Psi \neq \emptyset
    \implies \choice_{\Psi}: \one \to \Psi \subset \N\ \emph{total} \\
& (\text{choice available by $\ACC:$ non-empty sets have
        \emph{defined points}}) \\
& \implies \mu\set{d:t(d)\ \text{non-terminating}}:\one \to \Psi\ \emph{total.}
\end{align*}
This means: the {assumption} of (formal) \emph{existence} of a
$d \in \mr{DIO}$ for which decision race $t: \mr{DIO} \parto \N$ does \emph{not} terminate, leads to a (\emph{defined}) point
  $$d_0:\one \to \mr{DIO}$$
for which $t$ derivably does not terminate.

But this is \textbf{excluded} by pointwise non-derivability above
of non-termination, within frame $\bfQ$ assumed consistent.

So we have shown
\begin{align*}
& \bfQ,\PA+\ACC \derives\ \Psi = \emptyset,\ \text{i.\,e.} \\
& \bfQ \derives\ (\forall d\in \mr{DIO})[\exists k\,\ev(d,\mr{ct}_*k)=0 \\
& \qquad\qquad
  \lor \exists k\,\Pro_{\DIO}(k,\code{(\vec{x})\,\ev(d,\vec{x}}) \neq 0)],
\end{align*}
whence

\smallskip
\textbf{Termination Theorem:} $\bfQ,\ZFC,\PA+\ACC$ derive race $t$ to
terminate on all diophantine codes $d,$ on all $d \in \mr{DIO} = \Zlist.$

\section{Correct termination of decision $\nabla$}

\textbf{In particular} ($\bfQ^+ = \bfQ+\ACC$ stronger than $\bfQ$):
\begin{align*}
& \bfQ^+\ \textbf{derives} \\
& \quad
    \text{overall termination of $\mu$-recursive} \\
& \quad
    \text{termination race}\
                         t=t^{\bfQ}(d):\mr{DIO} \to \N: \\
& \bfQ^+\ \derives\
    [\,(\forall\,d \in \mr{DIO})\ t(d) \in \N
                   \ \text{{\emph{defined}}}\,]
\end{align*}

\textbf{Hence,} by Decision Correctness within $\bfQ^+:$

\smallskip
$\bfQ^+$ {$\mathbf{derives}$} \\
{overall \emph{correct} termination} of $\mu$-recursive \emph{decision}\\
$\nabla: \mr{DIO} \to \two,$ \textbf{main result} here:
\begin{align*}
& \nabla(d) \\
& = \begin{cases}
      0\ \myif\ \mr{ev}(d,t(d)) = 0 \\
      \qquad
        \ [ \ \implies d\ \emph{has}
                         \ \text{a zero}\ \vec{z} \in \Z^*\ ] \\
      1\ \myif\ \Pro_{\DIO}
              (t,\code{(\forall\,\vec{x})
                      \ \mr{ev}(d,\vec{x}) \neq 0}) \\
      \qquad
        \ [ \ \implies d\ \text{has}\ \emph{no}\ \text{zero}\ ]
    \end{cases}: \mr{DIO} \to \two.
\end{align*}

\section{Comparison with Matiyasevich's\\negative result}

\emph{Main result} above says in terms of the theory $\mathbf{TM}$
of TURING machines, by the established part of CHURCH's thesis:

\smallskip
\emph{For concrete diophantine polynomials} $D = D(\vec{x}): \Z^m \to \Z:$

For quantified arithmetical choice theories $\bfQ+\ACC$ like
$\ZFC$ and already $\PA+\ACC,$

$\bfQ^+ = \bfQ+\Con_{\bfQ}$ {$\mathbf{derives}$}:

\smallskip
\emph{TURING machine ${{\mrTM}}_{\nabla_{\bfQ}}$
corresponding---CHURCH---to
totally defined $\mu$-recursive decision map
  $$\nabla_{\bfQ}: \mr{DIO} \to \set{0,1},$$
\emph{when written} \emph{coefficient list} $\ccode{D}$
of a diophantine polynomial $D$ on its (initial) TAPE,
eventually \textbf{reaches} \emph{HALT state,}
leaves result $0$ (as its \emph{final} TAPE)
$\mathbf{iff}$ $D$ \emph{has} a zero\\
$\vec{z}:$ $D(\vec{z}) = 0,$
and \textbf{result} $1$ $\mathbf{iff}$ $D$ is overall \emph{non-null}: \\
$(\forall\,\vec{x} \in \Z^*)\,[\,D(\vec{x}) \neq 0\,].$}

\medskip
This contradicts \textbf{Matiyasevich's} THEOREM \emph{unsolving}
Hilbert's 10th problem, within theory $\bfQ^+$ which strengthens
his framework of Peano Arithmetic $\PA+\ACC$ with countable axiom
of choice. Whence

\subsubsection*{Conclusion:}
\begin{itemize}
\item
$\ZFC^+ = \ZFC+\ConZFC$ is contradictory, so

\item
$\ZFC\derives\neg\,\ConZFC:$ $\ZFC$ \emph{is internally inconsistent,}

\item
same for theory $\PA+\ACC:$

\emph{Peano-Arithmetic with axiom of countable choice is internally
inconsistent}

\item
\textbf{Question:} is already Peano Arithmetic $\PA$ by itself
internally inconsistent? It would be if axiom $\ACC$ of countable
choice were derivable within $\PA$ or independent from $\PA,$
as is axiom of choice $\AC$ from \textbf{set} theory. This would
mean that formal existential quantification is incompatible with
free-variables Primitive Recursive Arithmetic $\PR.$
\end{itemize}

\subsubsection*{Discussion}
\begin{itemize}
\item
After his talk at Humboldt University Berlin, I have mailed to
Matiyasevich the question, if his \emph{unsolving} of Hilbert's
10th problem is really constructive: it depends heavily on
formal existential quantification. No reply: may be he considers this
question when present paper will be brought to his attention.

\item
I have submitted the 200? version of present work, claiming
self-inconsistency $\PA\derives \neg\,\Con_{\PA},$ to the
\emph{Journal of Symbolic Logic.} The (anonymous) referee:

\emph{... this is certainly false. ...} Robert 'Rob' \NAME{Goldblatt} ed.:
\emph{under these circumstances etc.}
\end{itemize}

\textbf{What is such editorial policy good for?}

\section{Hilbert 10 constructively}

In this section we show that the \emph{local} version $\nabla[D]: 1 \parto 2$
of the $\mu$-recursive \emph{decision algorithm}
$\nabla = \nabla_{\DIO}(d): \mathit{DIO} \parto \two$
\emph{irrefutably} \emph{decides} \emph{each (single)} diophantine
equation---\emph{correctly}---when placed in p.\,r.
\emph{non-infinite-descent theory} $\piR=\PR+(\pi)$ of
op.\,cit.\  in the References.

This will give a positive solution to Hilbert's 10th problem
in that constructive framework, at least when
stated in its original form quoted in first section above.

\medskip
Formally, this \textbf{problem} allows for solution by a separate decision algorithm  (``process'') for each diophantine polynomial.

By \emph{localisation} at a given polynomial, we extract such a
decision-\emph{family} from the forgoing sections, and formalise it
within $\piR.$

We index that family (externally) by the
\emph{diophantine constants} $\delta:\one \to \mr{DIO} \subset \N,$
among which the diophantine polynomials
$$D = D(\vec{x}) = D(x_1, \ldots, x_{\bs{m}}): \Z^{\bs{m}} \to \Z$$
are represented by their coefficient list codes
$\ccode{D}: \one \to \mr{DIO}.$

\medskip
\textbf{Definition:} For PR predicates $\ph_0, \ph_1: A \times \N \to \two$
we define the \emph{race winner predicate}
 $$\mu_{\lor}[\ph_0,\ph_1]: A \to \two$$
between $\ph_0$ and $\ph_1$ slightly assymmetrically by
\begin{eqnarray*}
&& \mu_{\lor} [\ph_0,\ph_1] = \mu_{\lor} [\ph_0,\ph_1](a)\\
&& \defeq
(\mathit{dc} \circ (\ph_0,\ph_1))
  \parcirc (A \times \mu[\ph_0 \,\lor\, \ph_1]) \parcirc \Delta_A: \\
&& A \to A \times A \parto A \times \N \to \two \times \two
       \overset{\mathit{dc}} {\longrightarrow} \two, \ \text{with} \\
&& \mathit{dc} = \mathit{dc}(u,v): \two \times \two \parto \two
   \ \text{defined by} \\
&& \mathit{dc}(u,v) \defeq
\begin{cases}
0 \ \text{if}\ u = 1, \\
1 \ \text{if}\ u = 0\ \land\ v = 1, \\
\text{\emph{definably undefined} if}\ u = v = 0.
\end{cases}
\end{eqnarray*}
This (partial) race winner predicate
$\mu_{\lor}[\ph_0,\ph_1](a): A \parto \two$
is characterised---within $\bs{S}=\PR$ as well as in $\bs{S}=\piR$---by
\begin{align*}
  \bs{S} \derives\,
           & [\,\ph_0 (a,n)
                \,\land\, \underset{i<n} {\land}\,\neg\,\ph_1(a,n)
                 \implies \mu_{\lor}[\,\ph_0,\ph_1\,](a)=0\,] \\
           & \land\, [\,\ph_1(a,n)
             \,\land\,\underset{i \leq n}{\land}\,\neg\,\ph_0(a,n)
                \implies \mu_{\lor}[\,\ph_0,\ph_1\,](a)=1\,].
\end{align*}
We allow us to write for this intuitively---in classical terms of a (partial)
case-distinction:
$$
  \mu_{\lor}[\,\ph_0,\ph_1\,](a) =
  \begin{cases}
    0 \ \text{if}\  \mu\ph_0(a) < \infty
                       \,\land\,\mu\ph_0(a) \leq \mu\ph_1(a), \\
    1 \ \text{if}\  \mu\ph_1(a) < \infty
                             \,\land\,\mu\ph_1(a) < \mu\ph_0(a).
  \end{cases}
$$
Our decision family
  $$\nabla[\delta]: 1 \parto \two,\ \delta: \one \to \mr{DIO} \subset \N$$
now is defined in the present $\mu$-recursive
frame as this type of race winning, of PR search for a zero
(in the evaluation) of $\delta$ against PR search for a (first)
internal non-nullity \emph{proof} for (the evaluation) of $\delta,$
namely by
\begin{align*}
\nabla[\delta]
& \defeq \mu_{\lor} [\ph_0[\delta],\ph_1[\delta]]: 1 \parto \two,
  \ \text{with} \\
\ph_0[\delta](k)
& \defeq [\,\ev(\delta,\ct_*(k)) = 0\,]: \N \to \two, \\
\ph_1[\delta](k)
& \defeq \Pro_{\bs{S}}(k,\code{(\vec{x})\ev(\delta,\vec{x}) \neq 0}.
\end{align*}
Here
  $$\ev = \ev(d,x): \N \times \N \supset \mr{DIO} \times \Z^* \to \Z$$
is evaluation with the characteristic \textbf{evaluation property}
\begin{align*}
& \ev(\ccode{D},(x_1,\ldots,x_{\bs{m}}))
  = D(x_1, \ldots, x_{\bs{m}}): Z^{\bs{m}} \to \Z,
\end{align*}
realised by (iterated) \NAME{Horner}'s schema (each application reduces
the number of remaining variables by 1), or by ``brute force'' evaluation of monomials.

\subsection{Decision Correctness}

\textbf{Soundness Recall:} Main result of op.\,cit.\
in the References is (logical) \emph{soundness} of theory $\piR:$
\begin{itemize}
\item
For a (p.\,r.\,) predicate  $\chi = \chi(a): A \to \two$ we have
$$\piR\,\derives\ \Pro_{\piR}(k,\code{\chi})
  \implies \chi(a): \N \times A \to \two,$$
$a \in A$ free, meaning here \emph{for all} $a\in A,$ and $k \in \N$ free,
meaning here \emph{exists} $k\in\N.$ This entails

\item
\emph{$\PR$ soundness} of $\piR:$ For a p.\,r.\ predicate
$\chi = \chi(a): A \to \two,$
$$\piR\,\derives\ \Pro_{\PR}(k,\code{\chi})
  \implies \chi(a): \N \times A \to \two,$$
as well as in particular

\item
\emph{Diophantine soundness} of $\piR:$
for a diophantine polynomial $D=D(\vec{x}):\Z^*\to \two$
$$\piR\,\derives\ \Pro_{\piR}(k,\code{(\vec{x})D(\vec{x})\neq 0})
  \implies D(\vec{x})\neq 0,$$
$k\in \N,\ \vec{x}\in \Z^*$ free.

\item
Already $\PR^+ = \PR+\Con_{\PR}$ is diophantine sound.
This needs an extra Proof.

We consider here frame $\bs{S}=\piR,$
$$\piR^+ = \piR+\Con_{\piR} = \piR,$$
the latter by op.\,cit.\ equivalent to soundness of theory $\piR.$
\end{itemize}

Namely from PR Soundness we get the

\bigskip
\textbf{Local Correctness-Lemma} for $\nabla[\delta]$ in $\piR:$
The partial $\PR$-map $\nabla[\delta]:\one \parto \two$
has the following correctness properties:

\smallskip
$\piR \derives\ :$
\begin{itemize}
\item $\delta$ does not fall in
\emph{both} of the two defined-cases stated for $\nabla[\delta],$

\item
$\nabla[\delta] = 0
 \implies \ev(\delta,\ct_{*}\circ \mu\ph_0 [\delta])=0:$
$\delta$ is implied to have available a zero in its \emph{evaluation,}

\item
$\nabla[\delta] = 1  \implies \ev(\delta,\vec{x}) \neq_Z 0,$
$\vec{x}$ free in $\Z^*$:
$\delta$ is implied to be evaluated globally non-null, in particular:

\item By diophantine evaluation for
$D = D(x_1, \dots x_{\bs{m}}): \Z^* \to \Z$ diophantine:
\begin{itemize}
\item
$\nabla[D] := \nabla[\ccode{D}] = 0
 \implies D(\ct_*(\mu\ph_0 [\ccode{D}]))=0:$

$D$ is implied to have a zero, as well as

\item $\nabla[D]=1 \implies [\,D(\vec{x}) \neq 0\,],$
here again $\vec{x}$ free over $\Z^*:$

$D$ is implied to be globally non-null\ \,\textbf{\qed}
\end{itemize}
\end{itemize}

\subsection{Decision Termination}

The final question to treat for this---canonical---family
$$\nabla = \nabla_{\mbf{DIO}}[\delta]:\one \parto \two,
  \ \delta: \one \to \mr{DIO} \subset \N$$
of \emph{local}---$\mu$-recursive---decision algorithms, is
\emph{termination,} for each $\delta,$ in particular for
$\delta = \ccode{D},$ $D = D(\vec{x})$ diophantine.

\smallskip
\textbf{Assume} $\nabla[d_0]$ \emph{not} to terminate for a particular
\emph{constant} $d_0: \one \to \mr{DIO},$ in particular $d_0$ of
form $D_0 = D_0(\vec{x}).$

Since we argue here purely \emph{syntactically}---within the \emph{theory}
$\widehat{\bs{S}} \bs{\supset} \bs{S} \bs{=} \PR+(\mr{abstr})$ of \emph{partial}
p.\,r.\ maps---no modelling in mind except some primitive recursive
\emph{Meta}mathematics (these in turn g\"odelised within $\bs{S}$)---we
discuss the stronger assumption

$\nabla[d_0]$ $\T$-\emph{derivably}
does \emph{not} terminate for a given diophantine constant
$d_0: \one \to \mathit{DIO},$ $\T$ an extension of $\bs{S}.$

\medskip
This \textbf{assumption} reads:

$$\T \derives\,(k)\psi[d_0](k):$$
here $k$ is free over $\N,$ and the PR predicate
$\psi[d_0](k): \N \to \two$ is defined by
\begin{align*}
& \psi[d_0](k) = \psi_0[d_0](k) \land \psi_1[d_0](k)\ \text{with} \\
& \psi_0[d_0](k) = [\,\ev(d_0,\ct_*(k)) \neq 0 \,], \text{and} \\
& \psi_1[d_0](k) = \neg\,\Pro_{\T}(k,\code{\ev(d_0,\vec{x})\neq 0}).
\end{align*}
So the assumption (``of the contrary'') reads:
\begin{align*}
\T \derives\,
& [\,\ev(d_0,\ct_{*}(k) ) \neq 0 \,] \\
& \land \neg\,\Pro_{\T}(k,\code{(\vec{x})\ev(d_0,\vec{x}) \neq 0}).
\end{align*}
Here $k\in\N$ is the only free variable in the \emph{accessible} level,
$\vec{x}$ is free over $\Z^*,$ but \emph{encapsulated} within g\"odelisation,
\emph{not visible} on the object language level.

\smallskip
The derivably-non-termination assumption
$$\T \derives\, \psi[d_0](k),\ k\ \free,$$
would entail in particular (first conjunct $\psi_0[d_0]$):
$$\T \derives\, \ev(d_0,\ct_*(k)) \neq 0:\N\to\two.$$

\emph{Internalising} (\emph{formalising}) this metamathematical statement,
we (would) get by Proof-Internalisation---\cf \NAME{Smorynski} 1977---a
\emph{constant} $p_0: \one \to \mathit{Proof}_\T \subset \N$
\emph{guilty} for this last statement:
$$\T \derives\ \Pro_\T(p_0,\code{\ev(d_0,\vec{x}) \neq 0});$$
this would give, by definition of $\nabla[d_0]:$
$$\T \derives \nabla[d_0] = 1,$$
a contradiction to our assumption that $d_0$ be derivably
\emph{not decided} by $\nabla_{\DIO},$ \ie to $\T \derives \psi[d_0].$

\bigskip
\textbf{Conclusion:}

\begin{itemize}
\item $\piR=\piR+\Con_{\piR}$ derives the
alleged decision algorithm (family)
$\nabla = \nabla_{\DIO}[D]:\one \parto \two$ to be \emph{correct}
for each diophantine polynomial (if defined).

\item no diophantine polynomial $D=D(\vec{x})$
can come with a $\T$-proof (i.\,p.\ a $\piR$-proof) showing $\nabla[D]$
to be \emph{undefined,} \emph{not} to terminate,
in other words:

\item \emph{correct termination} of the
$\mu$-recursive \emph{decision family} $\nabla = \nabla_{\DIO}[D]$
at each diophantine polynomial is $\piR$-\emph{irrefutable,}
in the sense that \textbf{otherwise}---refutation---
\begin{align*}
& \piR \derives \Pro_{\piR}(q,\code{\mr{false}}),
  \ q: \one \to \N\ \text{a suitable PR point,}
\end{align*}
inconsistency of (self-consistent) theory $\piR$ would be the consequence.

\end{itemize}

\subsection*{Outlook}

Irrefutable correct termination of \emph{uniform} decision algorithm
$$\nabla_{\DIO}=\nabla_{\DIO}(d):\mr{DIO}\parto\two,\ d\in\mr{DIO}\ \free$$
is treated within the general framework of

\emph{Arithmetical Decision}
to come.


text:
This work focuses on modelling the motion of a spacecraft in a three-body configuration. This requires developing novel perturbation methods and comparing their performance to two well-established models of motion: the two-body problem (2BP) and the circular restricted three-body problem (CR3BP). As such, \textcolor{External}{both} are here presented in detail, together with all the necessary information regarding reference frames and associated dynamical structures.

The 2BP is employed for a quick estimation of trajectories and manoeuvres of a space object. This model is used when \textcolor{Internal}{the motion of said object is simplified as having only one main celestial body governing it. Trajectories computed in the 2BP are typically represented in a fixed, \textit{inertial reference frame}.} In contrast, the CR3BP is the highest fidelity model employed in this study, acting as the baseline with respect to which all trajectory errors are computed. \textcolor{Internal}{It is utilised to determine the motion of an object subject to the gravitational influence of two other great celestial bodies, which move in circles. Trajectories computed in the CR3BP are regularly represented in a rotating, \textit{synodical reference frame}.}

\textcolor{Internal}{When computing trajectories for any given space object, the main celestial body governing it is called the \textit{primary} or \textit{central body}. If the influence of an additional great celestial body is considered, both are named the \textit{primaries}---alternatively, the bigger body is termed again \textit{primary} and the smaller \textit{secondary}. If the space object whose motion is being modelled has a negligible mass when compared to the primaries, it is referred to as the \textit{massless particle}.}

\textcolor{Internal}{It is common occurrence to normalise the quantities involved in the computation of space trajectories, especially when assuming a three-body system. In this case, the unit of mass is taken to be $m_1 + m_2$ (subscripts $1$ and $2$ referring to the primary and secondary, respectively); the unit of length is selected to be the distance between the centres of these bodies; the unit of time is chosen such that the orbital period of the primary and the secondary about their common centre of mass is $2\pi$. Then, the universal constant of gravitation becomes $G$ = 1; it follows that the common mean motion of the primaries, $n$, is also unity.}

\section{Reference Frames}
\label{sec:rf}

Several different \textcolor{Internal}{coordinate} frames of reference are going to be used throughout this document, depending on the employed motion model. \textcolor{Internal}{The main ones are the inertial, the synodic and the Local Vertical, Local Horizontal (LVLH) reference frames.} This section describes their characteristics, together with the transformation matrices that can be used \textcolor{Internal}{to transform the state vector from one to another}.

\textcolor{Both}{Figure \ref{fig:srfvscar} illustrates the typical representation of the inertial and the synodic reference frames. The example motion is the one of asteroid 2011 MD, computed in the Sun-Earth system, in the CR3BP. Figure \ref{fig:srfvscar} a) depicts the inertial frame with the asteroid trajectory in red, together with the Earth's two-body motion with the Sun in blue; Figure \ref{fig:srfvscar} b) shows the same motion in the synodic reference frame, where the Earth appears fixed.}

\begin{figure}[h]
	\centering
	\begin{minipage}[t]{0.49\linewidth}
		\includegraphics[width=\textwidth]{Images/car.pdf}
		\caption*{\textcolor{External}{a) Inertial Cartesian frame: the red and blue lines correspond to the movement of the asteroid and the Earth, respectively}}
	\end{minipage}
	\hfill
	\begin{minipage}[t]{0.49\linewidth}
		\includegraphics[width=\textwidth]{Images/srf.pdf}
		\caption*{\textcolor{External}{b) Synodic reference frame: the red line depicts the asteroid's motion, while the Earth appears fixed}}
	\end{minipage}
	\caption{\textcolor{Both}{Motion of asteroid 2011 MD in the Sun-Earth system.} Earth is scaled to three times the Hill radius, for visibility\label{fig:srfvscar}}
\end{figure}

\subsection{Inertial Reference Frame}

An inertial reference frame is one in which Newton's first law of motion is valid. It follows that the centre of mass of the system is always at rest or in uniform motion: from this point, the X-axis points towards a pre-established reference direction (e.g. Vernal equinox), the \textcolor{Internal}{Z-axis is perpendicular to a reference plane (e.g. ecliptic) and the Y-axis completes the right-hand coordinate system}.

\textcolor{Internal}{This is the frame of reference commonly used to plot the 2BP for the motion of a massless particle around a primary. It is typically centred on the body that exerts the gravitational attraction---in the Sun-Earth system's case, the Sun. However, it is important to mention that the primary is not fixed in space, and the 2BP equations describe the motion around a fixed centre of attraction. As such, the adoption of the term \textit{Sun-centred} inertial reference frame from literature implies the assumption that the motion of the Sun is negligible.}

\subsection{Synodic Reference Frame}

The synodic or rotating frame is a non-inertial frame of reference typically employed for \textcolor{Internal}{restricted} three-body systems, in which the massless particle is deemed to have an insignificant mass when compared to the primaries. This frame is characterised by having the X-axis lying along the vector that connects the centres of mass of the two primaries. Therefore, their positions are kept constant with respect to this frame, which rotates at the same angular speed as the mean motion of the secondary around the primary.

\textcolor{Internal}{The quantities involved in this reference frame are regularly presented in non\-/dimensional units: the distance between the primaries is equal to one, the frame rotates with unit angular velocity, the normalised gravitational parameter is determined as $\mu = \frac{\mu_2}{\mu_1 + \mu_2}$ and the two primaries are located on the X-axis at the points ($-\mu$, 0) and ($1 - \mu$, 0), with the barycentre as the origin.}

Figure \ref{fig:rf} illustrates the rotating relationship between the synodic frame \textcolor{Internal}{($O_{R xyz}$)} centred on the barycentre and the Sun-centred inertial reference frame \textcolor{Internal}{($O_{I xyz}$)}, again using the Sun-Earth system as example, with the spacecraft as the massless particle. The angle $\gamma$ is positive in the counter-clockwise direction and is equivalent to the dimensionless time.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.6\linewidth]{Images/RF2.pdf}
	\caption{\label{fig:rf}Relationship between coordinate systems for the rotating and inertial frames}
\end{figure}

\textcolor{Internal}{In order to transform the position and velocity of a body from the inertial to the synodic frame, the first step is to move the frame's origin to the barycentre. This is shown in Eqs. \eqref{eq:trans1} and \eqref{eq:trans2}, where the superscript $(B)$ represents the transitional barycentric origin for a quantity in the inertial frame:}
\begin{align}\label{eq:trans1}
\begin{bmatrix} x_I^{(B)} & y_I^{(B)} & z_I^{(B)} \end{bmatrix}^T &= \begin{bmatrix} x_I & y_I & z_I \end{bmatrix}^T - \begin{bmatrix} \mu & 0 & 0 \end{bmatrix}^T\\
\begin{bmatrix} \dot{x}_I^{(B)} & \dot{y}_I^{(B)} & \dot{z}_I^{(B)} \end{bmatrix}^T &= \begin{bmatrix} \dot{x}_I & \dot{y}_I & \dot{z}_I \end{bmatrix}^T - \begin{bmatrix} 0 & \mu & 0 \end{bmatrix}^T\label{eq:trans2}
\end{align}

\textcolor{Internal}{The rotation of the inertial frame to the synodic frame} is done using the following \textcolor{Internal}{matrix}:
\begin{align}\label{eq:rot}
\text{\textbf{R}}^R_I = \begin{bmatrix}
\cos\gamma & \sin\gamma & 0 \\
-\sin\gamma & \cos\gamma & 0 \\
0 & 0 & 1\\
\end{bmatrix}
\end{align}

Finally, by employing Eqs. \eqref{eq:trans1} to \eqref{eq:rot}, the transformation of the position and velocity of a body from the inertial to the synodic frame becomes:
\begin{align}
\begin{bmatrix} x_R \\ y_R \\ z_R \\ \end{bmatrix} =
\text{\textbf{R}}^R_I\begin{bmatrix} x_I^{(B)} \\ y_I^{(B)} \\ z_I^{(B)} \\ \end{bmatrix}, \hspace{1cm}
\begin{bmatrix} \dot{x}_R \\ \dot{y}_R \\ \dot{z}_R \\ \end{bmatrix} =
\text{\textbf{R}}^R_I\begin{bmatrix} \dot{x}_I^{(B)} + y_I^{(B)} \\ \dot{y}_I^{(B)} - x_I^{(B)} \\ \dot{z}_I^{(B)} \\ \end{bmatrix}
\end{align}

\subsection{Local Vertical, Local Horizontal Reference Frame}
\label{sec:LVLH}

\textcolor{Internal}{The LVLH reference frame is a non-inertial frame of reference whose origin is placed at the centre of mass of the spacecraft. The axes are represented as $\vec{e}_r$, $\vec{e}_{\theta}$ and $\vec{e}_h$: respectively for the radial, in-track and cross-track directions. This reference frame can be visualised in Figure \ref{fig:lvlh}, together with the spacecraft's orbit in the inertial reference frame.}

\begin{figure}[hbt!]
	\centering
	\includegraphics[width=0.56\linewidth]{Images/LVLH.pdf}
	\caption{\label{fig:lvlh}Relationship between LVLH and inertial reference frames}
\end{figure}

\textcolor{Internal}{The rotation of the LVLH frame ($O_{r \theta h}$) to the orbital plane frame ($O_{eph}$) is done using Eq. \eqref{eq:battin_eq}. In turn, Eq. \eqref{eq:battin_eq} is used to convert vectors from the orbital plane ($O_{eph}$) to the inertial reference frame ($O_{I xyz}$), complying with the definitions in Battin \cite{battin}.}
\begin{align}\label{eq:battin_eq2}
\text{\textbf{R}}_{r \theta h}^{eph} =
\begin{bmatrix}
	\cos\nu & -\sin\nu & 0 \\
	\sin\nu & \cos\nu & 0 \\
	0 & 0 & 1 \\
\end{bmatrix}
\end{align}
\begin{align}\label{eq:battin_eq}
\text{\textbf{R}}_{eph}^{I} = \begin{bmatrix}
\cos\Omega & -\sin\Omega & 0 \\
\sin\Omega & \cos\Omega & 0 \\
0 & 0 & 1 \\
\end{bmatrix}
\begin{bmatrix}
1 & 0 & 0 \\
0 & \cos i & -\sin i \\
0 & \sin i & \cos i \\
\end{bmatrix}
\begin{bmatrix}
\cos\omega & -\sin\omega & 0 \\
\sin\omega & \cos\omega & 0 \\
0 & 0 & 1 \\
\end{bmatrix}
\end{align}

\section{The Two-Body Problem}
\label{sub:2bp}

The mathematical problem of designing a trajectory in the Solar System, an environment composed of countless bodies, is named the \textit{N-body problem}, in which $N$ \textcolor{External}{stands for the number of bodies interacting with object whose motion is being computed, including the latter}. Using Newton's Law of Gravitation, its equations of motion can be represented by Eq. \eqref{eq:nbody}:
\begin{align}
\ddot{\bm{r}_i} &= G \sum_{\substack{j=1 \\ j \neq i}}^{N} \frac{m_j}{r^3_{ij}}(\bm{r}_j - \bm{r}_i)\label{eq:nbody}
\end{align}
\textcolor{External}{in which the subscript $i$ singles out the object whose motion is being computed, $m$ represents each bodies' mass, $\bm{r}$ is the position vector, with $r = \norm{\bm{r}}$ and $r_{ij} = \norm{\bm{r_i} - \bm{r_j}}$.}

This problem has $6N$ variables, but only 10 constraint equations are obtained using conservation laws---6 equations for the conservation of momentum of the center of mass, 3 equations for the conservation of the total angular momentum and one equation for total energy conservation. Each constraint equation defines a conserved quantity, so-called \textit{integral of motion}. Thus, there are $6N - 10$ quantities to be determined, making the general $N$-body problem analytically not integrable without performing a range of approximations. This is one reason for models of motion of the Solar System, called full ephemerides models, to require statistical models and historical data in order to reach approximate solutions \cite{ephemerides}.

The 2BP is the case when $N$ is 2; this particular scenario describes the motion of two bodies that interact only with each other. \textcolor{Internal}{In most cases, only the motion of one of the objects relative to the other is of interest. This is the situation presented throughout this thesis, where the singled-out body is treated as a massless particle and the central one is deemed fixed.} The equations of motion are the particular case of Eq. \eqref{eq:nbody} for $N = 2$:
\begin{align}
\ddot{\bm{r}} &= -\mu_C \cdot \frac{\bm{r}}{r^3}\label{eq:2bp}
\end{align}
where $\bm{r}$ is the position vector relative to the central body, $\mu_C = G m_C$ is the gravitational parameter of the central body, equal to the product of the gravitational constant and the object's mass (subscript $C$ used to contrast with the case when two primaries exist, which bears no extra notation).

The 2BP is the basis of the patched-conics method: an approach used in trajectory design that divides space into areas in which an object is only gravitationally perturbed by one central body at a time, within each so-called \textit{sphere of influence} or Hill's sphere \cite{fundamentals}.

\section{The Circular-Restricted Three-Body Problem}

The three-body problem refers to the case where $N = 3$: three bodies move under their mutual gravitational attraction. As previously stated, this has no closed-form solutions. \textcolor{Internal}{A common} simplification can be applied: the CR3BP. In this model, the primaries move in circular, co-planar orbits around their common barycentre.

\subsection{Equations of Motion}

The full derivation of the equations of motion of the CR3BP is omitted here for simplicity's sake, since it can be found in many classical works \cite{szebehely1969theory}. Following the nomenclature defined by Koon et. al. \cite{koon2008}, the final model is depicted by Eq. \eqref{eq:cr3bp} in non-dimensional units in the synodic reference frame:
\begin{align}
\ddot{x} - 2\dot{y} &= -\frac{\partial \bar{U}}{\partial x}\nonumber\\
\ddot{y} + 2\dot{x} &= -\frac{\partial \bar{U}}{\partial y}\nonumber\\
\ddot{z} &= -\frac{\partial \bar{U}}{\partial z}
\label{eq:cr3bp}
\end{align}
where $\bar{U}$ is the effective potential function of the system, computed by Eq. \eqref{eq:pseudop}:
\begin{align}
\bar{U} &= -\frac{1}{2}(x^2 + y^2) - \frac{1 - \mu}{r_1} - \frac{\mu}{r_2}
\label{eq:pseudop}
\end{align}
\textcolor{Internal}{in which $r_1 = \norm{\bm{r_1}}$ and $r_2 = \norm{\bm{r_2}}$ are the distances from the massless particle to the primary and the secondary, respectively} (as depicted in Figure \ref{fig:rf}). These equations describe the state of the body in the synodic reference frame, generally with normalised components.

Given that the CR3BP is not analytically solvable, its solutions are obtained by employing numerical integrators. In contrast to the two-body case, the orbital elements of the trajectory will not be constant, instead changing over time. Some interesting properties and particular solutions of the motion under this model can be derived and will be detailed in the following sections.

\section{Numerical Targeting}
\label{sec:num_targ}

\textcolor{Internal}{As stated in Section \ref{sub:2bp}, there cannot be generic closed-form analytical solutions} to dynamical systems of higher complexity than the 2BP. So, when a specific orbital behaviour is sought and no analytical equations can be utilised, a two-point boundary value problem has to be solved by employing numerical targeting techniques.

The tools here presented are not unique to any motion model, but are essential in the computation of the invariant manifolds and periodic orbits of the CR3BP and any trajectory targeting presented throughout this work.

\subsection{State Transition Matrix}

\textcolor{Internal}{Given a set of initial conditions, the numerical integration process of a motion model will propagate the state in time. Yet, the final condition obtained by this means may not necessarily be the one that achieves the trajectory design goals (e.g. targeting of a specific point in space, rendezvous with an object, computation of a periodic orbit). Thus, in order to reach the desired target, the initial conditions of the propagated motion may require modification. To correctly and efficiently adjust a trajectory, the computation of the state transition matrix (STM) is necessary.}

\textcolor{Internal}{In order to illustrate how a trajectory can be adjusted, Figure \ref{fig:stm} depicts a two motions: one of them, highlighted in purple, is a naturally propagated reference trajectory $\bm{s}_R(t)$, from an initial condition $\bm{s}_R(t_0)$ to a final state $\bm{s}_R(t_1)$; the other, coloured red, is the perturbed trajectory $\bm{s}_P(t)$, computed when a perturbation $\delta \bm{s}(t_0)$ is introduced in the reference motion. After integrating this perturbed state in time, the resulting trajectory is given by:}
\begin{align}
\textcolor{Both}{\bm{s}_P(t) = \bm{s}_R(t) + \delta \bm{s}(t)}\label{eq:stm0}
\end{align}

\textcolor{Internal}{Integrating both the reference and the perturbed initial states in the CR3BP results in two distinct paths---at time $t_1$, the states along these two paths are not equal. In order to estimate their difference $\delta \bm{s}(t_1)$, the motion is linearised using a first-order Taylor series expansion about the reference trajectory.}

The perturbed arc is expressed by Eq. \eqref{eq:stm1}.
\begin{align}
\textcolor{Both}{\dot{\bm{s}}_P = \dot{\bm{s}}_R + \delta \dot{\bm{s}}  = f(\bm{s}_R + \delta \bm{s})}\label{eq:stm1}
\end{align}
where the function $f$ represents the non-linear differential equations of the motion model. The Taylor series expansion on Eq. \eqref{eq:stm1} yields the first-order terms:
\begin{align}
\dot{\bm{s}}_P &\approx f(\bm{s}_R, t) + \delta \dot{\bm{s}}\\
\delta \dot{\bm{s}} &\approx \frac{\partial f}{\partial \bm{s}} \bigg\rvert_{\bm{s}_R} \delta \bm{s}\label{eq:stm2}
\end{align}

Eq. \eqref{eq:stm2} has solutions of the form:
\begin{align}
\delta \bm{s}(t_1) &= \bm{\Phi}(t_1, t_0)\delta \bm{s}(t_0)
\end{align}
where $\bm{\Phi}(t_1, t_0)$ denotes the STM. Thus, the STM provides a linear mapping from time $t_0$ to a time $t_1$, establishing a relationship between initial and final deviations that can be used to adjust trajectories in order to match a final outcome. The propagation of the STM is described by the following equation:
\begin{align}
\textcolor{Both}{\dot{\bm{\Phi}}} &= \text{\textbf{D}}f(\bm{s}) \bm{\Phi} \label{eq:1}
\end{align}
in which $\text{\textbf{D}}f(\bm{s})$ is the Jacobian of the STM.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\linewidth]{Images/SingleShooting.pdf}
	\caption{\label{fig:stm}Reference path and perturbed solution at times $t_0$ and $t_1$}
\end{figure}

\subsection{Differential Correction}

Given an initial and final set of trajectory parameters, the goal of a differential corrector is to use the STM to determine which modifications of the initial conditions are required in order to get a closer approach to the desired state at the end of the integration period. Thus, this method can be used to determine the necessary orbital manoeuvres to \textcolor{Internal}{be performed} on an initial state \textcolor{Both}{$\bm{s}_0$}, so that an adjusted transfer will end up at the desired final state \textcolor{Both}{$\bm{s}_1$} \cite{koon2008}. \textcolor{Internal}{Given a} reference trajectory governed by the non-linear dynamical equations $f$:
\begin{align}
\textcolor{Both}{\dot{\bm{s}}} &= \textcolor{Both}{f(\bm{s}(t))} \\
\textcolor{Both}{\delta \bm{s}(t_1)} &= \textcolor{External}{\bm{\Phi}(t_1, t_0)\delta \bm{s}(t_0)}
\end{align}
in which $\bm{\Phi}$ is the state transition matrix (STM), whose propagation is described by Eq. \eqref{eq:1}.

The Jacobian \textcolor{Both}{$\text{\textbf{D}}f(\bm{s})$} is composed of the partial derivatives of the velocity and acceleration vectors, in the form of Eq. \eqref{eq:bigmatrix}:
\begin{align}
\text{\textbf{D}}f(\bm{s}) = \begin{bmatrix}
\frac{\partial \dot{x}}{\partial x} & \frac{\partial \dot{x}}{\partial y} & \frac{\partial \dot{x}}{\partial z} & \frac{\partial \dot{x}}{\partial \dot{x}} & \frac{\partial \dot{x}}{\partial \dot{y}} & \frac{\partial \dot{x}}{\partial \dot{z}}\\
\frac{\partial \dot{y}}{\partial x} & \frac{\partial \dot{y}}{\partial y} & \frac{\partial \dot{y}}{\partial z} & \frac{\partial \dot{y}}{\partial \dot{x}} & \frac{\partial \dot{y}}{\partial \dot{y}} & \frac{\partial \dot{y}}{\partial \dot{z}}\\
\frac{\partial \dot{z}}{\partial x} & \frac{\partial \dot{z}}{\partial y} & \frac{\partial \dot{z}}{\partial z} & \frac{\partial \dot{z}}{\partial \dot{x}} & \frac{\partial \dot{z}}{\partial \dot{y}} & \frac{\partial \dot{z}}{\partial \dot{z}}\\
\frac{\partial \ddot{x}}{\partial x} & \frac{\partial \ddot{x}}{\partial y} & \frac{\partial \ddot{x}}{\partial z} & \frac{\partial \ddot{x}}{\partial \dot{x}} & \frac{\partial \ddot{x}}{\partial \dot{y}} & \frac{\partial \ddot{x}}{\partial \dot{z}}\\
\frac{\partial \ddot{y}}{\partial x} & \frac{\partial \ddot{y}}{\partial y} & \frac{\partial \ddot{y}}{\partial z} & \frac{\partial \ddot{y}}{\partial \dot{x}} & \frac{\partial \ddot{y}}{\partial \dot{y}} & \frac{\partial \ddot{y}}{\partial \dot{z}}\\
\frac{\partial \ddot{z}}{\partial x} & \frac{\partial \ddot{z}}{\partial y} & \frac{\partial \ddot{z}}{\partial z} & \frac{\partial \ddot{z}}{\partial \dot{x}} & \frac{\partial \ddot{z}}{\partial \dot{y}} & \frac{\partial \ddot{z}}{\partial \dot{z}}\\
\end{bmatrix}\label{eq:bigmatrix}
\end{align}

Considering these mathematical relations, the process to obtain a differential corrector can be detailed. First of all, the objective is to find a solution $\bm{X}$ that satisfies the following expression,
\begin{align}\label{eq:GX}
\text{\textbf{G}}(\bm{X}) = \begin{bmatrix} G_1(\bm{X}) \\ ... \\ G_m(\bm{X}) \end{bmatrix} = 0
\end{align}
in which $\bm{X}$ is a vector of free variables and $\text{\textbf{G}}(\bm{X})$ is the vector of $m$ constraints to which the free variables are subject, so that the trajectory meets the required specifications. This equation can be expanded as a first-order approximation using a Taylor series about an initial guess $\bm{X}_0$:
\begin{align}
\text{\textbf{G}}(\bm{X}) = \text{\textbf{G}}(\bm{X}_0) + \frac{\partial \text{\textbf{G}}(\bm{X}_0)}{\partial \bm{X}_0}(\bm{X} - \bm{X}_0)
\end{align}

Given that the constraint vector should be zero per the equality in Eq. \eqref{eq:GX}, Eq. \eqref{eq:lsm} is obtained:
\begin{align}
\bm{X}^{i + 1} = \bm{X}^{i } - \text{\textbf{D}}\text{\textbf{G}}(\bm{X}^{i})^{-1} \text{\textbf{G}}(\bm{X}^{i})\label{eq:lsm}
\end{align}
in which $\text{\textbf{D}} \text{\textbf{G}}(\bm{X}) = \frac{\partial \text{\textbf{G}}(\bm{X})}{\partial \bm{X}}$ and represents the Jacobian matrix of the constraint vector and superscript $i$ represents each iteration of the differential correction algorithm. Thus, an iterative process is required to cycle through the error ($|\text{\textbf{G}}(\bm{X}^{i+1})|$) until it falls below a pre-defined tolerance ($\epsilon$). The different steps of this iterative process are outlined below \cite{sanchez}:

\begin{enumerate}
	\item Define the problem to solve and identify the free variables $\bm{X}$
	\item Determine a first guess of the free variables vector $\bm{X}_0$
	\item Specify constraints to which variables are subject $\text{\textbf{G}}(\bm{X})$
	\item Calculate the Jacobian matrix of the system $\text{\textbf{D}} \text{\textbf{G}}(\bm{X})$
	\item Solve the system of equations $\bm{X}^{i+1} = \bm{X}^{i} - \text{\textbf{D}} \text{\textbf{G}}(\bm{X}^{i})^{-1} \text{\textbf{G}}(\bm{X}^{i})$
	\item Check the error associated to the new solution ($|\text{\textbf{G}}(\bm{X}^{i+1})|$). If it is larger than a convergence tolerance ($\epsilon$), repeat step 5 using $\bm{X}^{i+1}$ as the new $\bm{X}^i$. If the error falls below $\epsilon$, the method has converged.
\end{enumerate}

\subsubsection{Single Shooting Method}
\label{sec:ssm}

A single shooting method is a differential correction technique to connect the initial and final points of a trajectory with one single segment. It is suitable for the computation of simple trajectories (i.e. periodic orbits or direct transfers between two points). The algorithm can employ a fixed or variable time approach, although the mission design in this project concerns only the first case.

One commonly used approach is the targeting of the final position by changing the initial velocity. This is especially important in the computation of Lambert arcs, in order to compute the $\Delta v$ necessary to move the spacecraft to a certain position in space.

In this scenario, the free variables vector $\bm{X}$ consists of the initial velocity components, which are the only design parameters that can be modified to reach the desired final state \textcolor{Both}{$\bm{s}_d$:}
\begin{align}
\bm{X} &= \begin{bmatrix}\dot{x}_0 & \dot{y}_0 & \dot{z}_0 \end{bmatrix}^T\\
\textcolor{Both}{\bm{s}_d} &= \begin{bmatrix} x_d & y_d & z_d \end{bmatrix}^T
\end{align}

The Jacobian of the system can be determined by computing the partial derivatives of the constraint vector $\text{\textbf{G}}(\bm{X})$. These correspond, in this case, to the section of the STM associated with the velocity of the initial state vector and the position of the final state.

\subsubsection{Multiple Shooting Method}
\label{sec:MSM}

A multiple shooting method is a differential correction technique utilised to connect the initial and final points of a trajectory with multiple segments, i.e. using several single shooting methods between intermediate points, the so called patch points. In this way, the errors associated with the long integration of segments are reduced. Plus, the shape of the trajectory can be more easily manipulated by adding constraints to the different patch points.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.8\linewidth]{Images/MultipleShooting.pdf}
	\caption{\label{fig:ms}Multiple shooting scheme}
\end{figure}

In the case in which the integration times are fixed, \textcolor{External}{the vector of free variables} corresponds to the linkage of each patch point's state vector \textcolor{External}{($\bm{s}_j, j = 1, ... , m$)}.
\begin{align}
\bm{X} = \begin{bmatrix}\textcolor{External}{\bm{s}_1} & ... & \textcolor{External}{\bm{s}_m} \end{bmatrix}^T
\end{align}

The free variables vector $\bm{X}$ becomes a column vector of $6m$ rows if only position and velocity are considered, or $9m$ rows if the acceleration is added (in which $m$ is the number of patch points). To ensure position or velocity continuity at the patch points, conditions have to be added to the constraint vector in the following manner:
\begin{align}
\text{\textbf{G}}(\bm{X}) =
\begin{bmatrix}
\textcolor{External}{\bm{s}_1(t_1 + T_1) - \bm{s}_2(t_2)} \\...\\ \textcolor{External}{\bm{s}_{m-1}(t_{m-1} + T_{m-1}) - \bm{s}_m(t_m)}
\end{bmatrix}
\end{align}

Besides the named constraints, any additional ones can be added, including fixing the initial and final states of the motion.

Eq. \eqref{eq:exampleG} shows the first row of the Jacobian of the constraint vector:
\begin{align}\nonumber
\frac{\partial \text{\textbf{G}}_1(\bm{X})}{\partial \bm{X}} &=
\begin{bmatrix}
\frac{\partial \text{\textbf{G}}_1}{\partial X_1} & ... & \frac{\partial \text{\textbf{G}}_1}{\partial X_{m}}
\end{bmatrix} =
\begin{bmatrix}
\textcolor{External}{\frac{\partial \bm{s}_1(t_1 + T_1)}{\partial \textcolor{External}{\bm{s}_1(t_1)}}} & \textcolor{External}{-\frac{\partial \bm{s}_2(t_2)}{\partial \bm{s}_2(t_2)}} & ... & \textcolor{External}{\frac{\partial G_1}{\partial \bm{s}_m(t_m)}}
\end{bmatrix} \\ &=
\begin{bmatrix}
\bm{\Phi}_1 & -\text{\textbf{I}} & ... & 0
\end{bmatrix}\label{eq:exampleG}
\end{align}

The same procedure applies to the remaining rows of $\text{\textbf{G}}(\bm{X})$. As such, every segment will contribute with its own STM and one identity matrix. Finally, the entire Jacobian for all the constraints of the multiple shooting method can be assembled considering all the segments:
\begin{align}\label{eq:bigJ}
\text{\textbf{DG}}(\bm{X}) =
\begin{bmatrix}
\bm{\Phi}_1 & -\text{\textbf{I}} & \text{\textbf{0}} & 0 & 0 & 0 \\
\text{\textbf{0}} & \bm{\Phi}_2 & -\text{\textbf{I}} & 0 & 0 & 0\\
&  &  & \ddots &  & \\
0 & 0 & 0 & 0 & \bm{\Phi}_m & -\text{\textbf{I}}\\
\end{bmatrix}
\end{align}
in which each STM is computed using an adequate Jacobian matrix.

Given the following auxiliary variable:
\begin{align}
r_3 &= \sqrt{1 -2x + x^2 + y^2 + z^2}
\end{align}
and using the previously defined quantities $r$, $r_1$ and $r_2$, the Jacobian matrix of the CR3BP model \textcolor{External}{is the following}:
\begin{align}\label{eq:jac_cr}
\text{\textbf{D}}f(\bm{s}) = \begin{bmatrix}
0 & 0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 1 \\
\textcolor{External}{-\partial_{xx}\bar{U}} & \textcolor{External}{-\partial_{xy}\bar{U}} & \textcolor{External}{-\partial_{xz}\bar{U}} & 0 & 2 & 0\\
\textcolor{External}{-\partial_{yx}\bar{U}} & \textcolor{External}{-\partial_{yy}\bar{U}} & \textcolor{External}{-\partial_{yz}\bar{U}} & -2 & 0 & 0\\
\textcolor{External}{-\partial_{zx}\bar{U}} & \textcolor{External}{-\partial_{zy}\bar{U}} & \textcolor{External}{-\partial_{zz}\bar{U}} & 0 & 0 & 0\\
\end{bmatrix}
\end{align}
where the values of \textcolor{External}{$\partial_{\alpha\beta}\bar{U}$ ($\alpha \in \{x, y, z\} \text{ and }\beta \in \{x, y, z\}$)} are the second partial derivatives of the effective potential presented in Eq. \eqref{eq:pseudop}. These are expanded into the following:
\begin{align}
\partial_{xx}\bar{U} &= -1 - \frac{3\mu(x - 1 + \mu)^2}{r_2^5} + \frac{\mu}{r_2^3} - \frac{3 (1 - \mu)(x + \mu)^2}{r_1^5} + \frac{1 - \mu}{r_1^3}\\
\partial_{xy}\bar{U} &= -3y\bigg(\frac{\mu(x - 1 + \mu)}{r_2^5} + \frac{(1 - \mu)(x + \mu)}{r_1^5}\bigg) \\
\partial_{xz}\bar{U} &= -3z\bigg(\frac{\mu(x - 1 + \mu)}{r_2^5} + \frac{(1 - \mu)(x + \mu)}{r_1^5}\bigg) \\
\partial_{yx}\bar{U} &= \partial_{xy}\bar{U} \\
\partial_{yy}\bar{U} &= -1 - \frac{3 y^2 \mu}{r_2^5} + \frac{\mu}{r_2^3} - \frac{3 y^2 (1 - \mu)}{r_1^5} + \frac{1 - \mu}{r_1^3}\\
\partial_{yz}\bar{U} &= -3yz\bigg(\frac{\mu}{r_2^5} + \frac{(1 - \mu)(x + \mu)}{r_1^5}\bigg) \\
\partial_{zx}\bar{U} &= \partial_{xz}\bar{U} \\
\partial_{zy}\bar{U} &= \partial_{yz}\bar{U} \\
\partial_{zz}\bar{U} &= -\frac{3 z^2 \mu}{r_2^5} + \frac{\mu}{r_2^3} - \frac{3 z^2 (1 - \mu)}{r_1^5} + \frac{1 - \mu}{r_1^3}
\end{align}

\section{Particular Solutions}
\label{sub:ps}

\textcolor{Internal}{Although} the non-linear equations of motion in the CR3BP are non-integrable and autonomous (not time-dependant), some particular solutions can be derived: specifically, equilibrium points, periodic orbits and quasi-periodic orbits \cite{parker}.

\subsection{Equilibrium Points}

The libration points \textcolor{Internal}{of the CR3BP} are locations in space where an object of negligible mass, affected by the gravitational interactions between the primaries, can theoretically maintain a constant position in the synodic reference frame. This characteristic makes them very attractive for a great number of missions, e.g. to keep telescopes or other observation-type spacecraft, since the fuel consumption required to perform station-keeping is very low \cite{koon2008}.

The five libration points are the equilibrium solutions of the equations of motion of the CR3BP; they can be obtained by setting the velocity and acceleration of Eq. \eqref{eq:cr3bp} to zero:
\begin{align}
0 &= x - \frac{(1 - \mu)(x + \mu)}{r_1^3} - \frac{\mu(x + \mu - 1)}{r_2^3}\nonumber\\
0 &= y\bigg(1 - \frac{(1 - \mu)}{r_1^3} - \frac{\mu}{r_2^3}\bigg)\nonumber\\
0 &= -z\bigg(\frac{(1 - \mu)}{r_1^3} - \frac{\mu}{r_2^3}\bigg)\label{eq:lpoints}
\end{align}

By analysing Eq. \eqref{eq:lpoints}, some insights can be taken. These equations imply that $z = 0$, meaning the equilibrium solutions are on the XY plane of motion. Furthermore, from setting $y = 0$, three solutions can be achieved: these represent the collinear points $L_1$, $L_2$ and $L_3$. By defining $r_1 = r_2 = 1$, the remaining solutions are obtained: points $L_4$ and $L_5$, which form an equilateral triangle with the primaries. Their position can be visualised in Figure \ref{fig:lpoints}, depicted in the synodic reference frame.

\begin{figure}[!hbt]
	\centering
	\includegraphics[width=0.6\linewidth]{Images/LibPoints.pdf}
	\caption{\label{fig:lpoints}\textcolor{Internal}{Distribution of the libration points in the CR3BP}}
\end{figure}

\subsection{Jacobi Constant}

The equations of the CR3BP have an energy integral of motion. Following the notation in Koon \cite{koon2008}, its formula is depicted by:
\begin{align}
E = \textcolor{Internal}{\frac{1}{2}(\dot{x}^2 + \dot{y}^2 + \dot{z}^2)} + \bar{U}\label{eq:Jacobi}
\end{align}
with $\bar{U}$ computed by Eq. \eqref{eq:pseudop}.

The quantity $-2E$ is the one generally found in literature: it is named the Jacobi constant. This is also typically represented by the variable $C$, by employing the following reformulation of Eq. \eqref{eq:Jacobi}:
\begin{align}
C = -2\bar{U} - \norm{\bm{v}}^2\label{eq:C}
\end{align}
\textcolor{Internal}{where $v$ is the norm of the spacecraft's velocity.}

The Jacobi constant can be very helpful in the characterisation of the system's dynamics, namely to determine the regions in space that are accessible to a theoretical spacecraft. Given that ${v^2 \geq 0}$ is a mathematical constraint of Eq. \eqref{eq:C}, the set of variables which obey this condition forms the accessible region, where the spacecraft can move. By setting $v = 0$, and therefore $-2\bar{U} - C > 0$, the obtained solution is a region that delimits the space in which the body can possibly move. This is also called the zero-velocity curve.

Figure \ref{fig:ZVC} shows the zero-velocity curves for several increasing energy levels (i.e. decreasing Jacobi constant) in the Earth-Moon system. \textcolor{External}{Figure \ref{fig:ZVC} a)} shows the energy level of $L_1$, where there can be no motion between primaries. \textcolor{External}{In Figure \ref{fig:ZVC} b)} ($C = 3.172$), a trajectory from the Earth to the Moon is possible for the first time, showing that flying to the Moon via $L_1$ is energetically most favourable. The energy corresponding to \textcolor{External}{Figure \ref{fig:ZVC} c)} ($C = 3.162$) allows the test body to reach the $L_2$ point. \textcolor{External}{Figure \ref{fig:ZVC} d)} ($C = 3.012$) shows the energy of the test body increased to the level of $L_3$. If the energy is further increased, the $L_3$, $L_4$ and $L_5$ points can finally be achieved (\textcolor{External}{Figures \ref{fig:ZVC} e) and f)}, with $C = 2.989$ and $C = 2.979$).

\begin{figure}[!htb]
	\centering
	\begin{minipage}[t]{0.4\textwidth}
		\centering
		\includegraphics[width=\textwidth]{Images/L1.png}\caption*{\label{fig:L1}a) C = 3.188}
	\end{minipage}
	\begin{minipage}[t]{0.4\textwidth}
		\centering
		\includegraphics[width=\textwidth]{Images/L2.png}\caption*{\label{fig:L2}b) C = 3.172}
	\end{minipage}
	\begin{minipage}[t]{0.4\textwidth}
		\centering
		\includegraphics[width=\textwidth]{Images/L2bridge.png}\caption*{\label{fig:L2bridge}c) C = 3.162}
	\end{minipage}
	\begin{minipage}[t]{0.4\textwidth}
		\centering
		\includegraphics[width=\textwidth]{Images/L3.png}\caption*{\label{fig:L3}d) C = 3.012}
	\end{minipage}
	\begin{minipage}[t]{0.4\textwidth}
		\centering
		\includegraphics[width=\textwidth]{Images/L5.png}\caption*{\label{fig:L5}e) C = 2.989}
	\end{minipage}
	\begin{minipage}[t]{0.4\textwidth}
		\centering
		\includegraphics[width=\textwidth]{Images/zero.pdf}\caption*{\label{fig:L0}f) C = 2.979}
	\end{minipage}
	\caption{Zero-velocity curves in the Earth-Moon system with decreasing Jacobi constant. Axes in \textcolor{Internal}{normalised units (unit length equal to the Earth-Moon distance)}}\label{fig:ZVC}
\end{figure}

\subsection{Libration Point Orbits}

The CR3BP contains a wide variety of periodic and quasi-periodic orbits. The former retrace their path over time: at the end of a certain time frame, the orbit's position and velocity repeat themselves. In contrast, quasi-periodic orbits have their motion confined to a particular region in space, changing slightly with each period \textcolor{Internal}{but never repeating the same orbit}. In the vicinity of the libration points, these are called the libration point orbits (LPOs).

Three types of \textcolor{Internal}{LPOs} can be highlighted: horizontal Lyapunov orbits, which are in the fundamental plane, vertical Lyapunov orbits, that are horizontally symmetric and shaped like a figure-eight, and halo orbits. Examples of these three orbital categories can be observed in Figure \ref{fig:lpo} \textcolor{Internal}{in the synodic reference frame}, for the same energy level.

\begin{figure}[!htb]
	\centering
	\includegraphics[width=0.7\linewidth]{Images/LPOv2.png}
	\caption{\label{fig:lpo}\textcolor{Internal}{Three types of LPOs in the Sun-Earth system, for C = 3.0007.}}
\end{figure}

\textcolor{External}{The computation of periodic orbits, as well as any trajectory with a set of desired characteristics in the CR3BP, may require the definition of a two-point boundary value problem, solved with the application of a differential corrector. The latter is} a method that, using the STM, determines how to change an initial trajectory in order to target a final set of orbital conditions.

\subsection{Invariant Manifolds}

\textcolor{Both}{Manifold theory is essential to the understanding of the dynamical environment of the CR3BP. The invariant manifold structures connected to LPOs are particularly interesting for mission design, since they can be travelled by a spacecraft without any fuel consumption.}

Invariant manifolds are dynamical structures composed of countless orbits. They exist for a range of energies and form a series of 'tubes' that connect different regions around the primaries, as it can be seen on Figure \ref{fig:inv_man}. To compute them, it is necessary to acquire information about the local stability characteristics of each point along the LPO. This is obtained by computing the eigenvalues and eigenvectors of the monodromy matrix \cite{koon2008}, defined as the STM when propagated for precisely one period of the orbit.

\begin{figure}[h]
	\centering
	\includegraphics[width=\linewidth]{Images/stableplusun.pdf}
	\caption{\label{fig:inv_man}Invariant manifold orbits connected to a halo orbit around the $L_1$ point of the Sun-Earth system, \textcolor{Internal}{C = 3.0007. On the left: stable manifold. On the right: unstable manifold}}
\end{figure}

\textcolor{Internal}{Therefore, the computation of the invariant manifold structures of an LPO starts with the determination of its monodromy matrix. The latter possesses a certain number of distinct eigenvalues $\lambda$, with corresponding eigenvectors. The eigenvalues that have non-zero real parts are associated with the hyperbolic invariant manifolds, and the stability of the latter is determined by their real sign. If $Re(\lambda) < 0$, the manifold is stable; if $Re(\lambda) > 0$, the manifold is unstable.}

\textcolor{Internal}{After the determination of the associated eigenvectors and eigenvalues, the LPO is discretised into numerous initial conditions, $\bm{s}_{P_0}$. A small perturbation step $\epsilon$ is introduced in the direction of the eigenvectors at each fixed point:}
\begin{align}\label{eq:stable}
\bm{s}_{S_0} = \bm{s}_{P_0} \pm \epsilon \bm{l}_S\\
\bm{s}_{U_0} = \bm{s}_{P_0} \pm \epsilon \bm{l}_U\label{eq:unstable}
\end{align}
where $\bm{l}_S$ and $\bm{l}_U$ are the stable and unstable eigenvectors computed from the monodromy matrix at the fixed point of interest, while $\epsilon$ is typically chosen to be equal to $10^{-6}$ (in normalised units). \textcolor{Internal}{The displacement $\epsilon$ is sufficiently small to avoid violating the linear approximation, but large enough to allow the manifold to depart the LPO after a reasonable time interval. Then, the the initial conditions of Eq. \eqref{eq:stable} are numerically integrated backwards in time to obtain the stable manifold set, while the ones of Eq. \eqref{eq:unstable} are, in contrast, integrated forwards in time to generate the unstable set. The explanations by Koon et. al. \cite{koon2001, koon_heteroclinic} are very helpful in understanding and applying the presented manifold theory.}

\section{Perturbation Methods}

Perturbation methods are a class of mathematical techniques used to generate solutions that describe the motion of a body subject to disturbing forces. \textcolor{Internal}{Depending on the orbit that is being analysed, the latter can be caused by the non-spherical shape of planets, atmospheric drag, solar radiation pressure, the gravitational attraction of other celestial bodies or other disturbing forces.} There are many techniques counted as perturbation methods, which can then be distinguished based on several different features.

The techniques to solve perturbation equations fall into three broad categories: analytical, semi-analytical and numerical solutions \cite{vallado}. The analytical approaches were, historically, developed first due to the lack of computational power before the 20\textsuperscript{th} century. Semi-analytical and numerical methods employ computational routines to solve the differential equations. Still, even with the current computing capabilities, analytical solutions are sometimes preferable for very fast or real-time computations, provided they respect a certain level of accuracy.

As such, perturbation methods can be divided in two main categories: \textit{special} and \textit{general} approaches. The first concerns the cases where the equations are solved using numerical methods. A specific, or \textit{special} answer is produced, only valid for the given data and initial conditions. Although special perturbation methods tend to be very accurate, they often suffer from over-specificity, since they cannot easily be extrapolated to different data. In contrast, general approaches employ anaytical and semi-analytical techniques, which can be used in several distinct case scenarios.

When describing different perturbation methods, some important characteristics should also be outlined. One refers to the difference between \textit{fast} and \textit{slow} variables. Fast variables change greatly during an orbit, even in the absence of perturbations (e.g. true, eccentric and mean anomalies). Slow variables change very little and are not modified at all if there is no perturbing effect (e.g. semi-major axis, eccentricity, inclination, argument of the periapsis and longitude of the ascending node).

The perturbed motion of a spacecraft can be computed using either \textit{mean} or \textit{osculating elements}. Mean elements are averaged over a chosen interval of time or true, eccentric, or mean anomaly, so they are relatively smoothly varying. By doing this, the perturbation is averaged throughout the motion instead of computed at each time-step. Single-averaged elements result from removing the high-frequency short-periodic motions, while double-averaged elements remove also the long periodic variations. On the other hand, osculating elements include all periodic (long and short-periodic) and secular effects. The osculating orbit is equivalent to the two-body orbit that the spacecraft would follow if the perturbing forces were suddenly removed at that instant: they are time-varying, so each trajectory point has a corresponding set of osculating elements. This representation is useful for highly accurate simulations, including real-time pointing and tracking operations.

\subsection{Variation of Parameters}
\label{sub:VOP}

The variation of parameters (VOP) method is a mathematical tool to solve first-order differential equations, which yields a formulation of the equations of motion to use in perturbed dynamical systems. This theory is based on the premise that the solution for the unperturbed system can be used to \textcolor{Internal}{find} the perturbed one. Concretely, the unperturbed system is governed by the two-body problem discussed in Section \ref{sub:2bp}; when a small perturbative force is added, its effect is described by equations of motion that represent the time-varying, osculating orbital elements \cite{vallado}. These equations were originally developed by Euler and Lagrange as a general perturbation method.

For the perturbed problem, the equations of motion of the 2BP change from Eq. \eqref{eq:2bp} to the general form:
\begin{align}
\ddot{\bm{r}} + \mu_C \cdot \frac{\bm{r}}{r^3} = \nabla \mathpzc{R} \label{eq:pert}
\end{align}
where $\mathpzc{R}$ is the disturbing function that fully describes the perturbation. This equation can be further expanded into:
\vspace{-5mm}
\begin{align}
\ddot{x}+ \textcolor{Internal}{\mu_C}\frac{x}{r^{\,3}} = \frac{\partial {\mathpzc R}}{\partial x} = a_x\nonumber\\
\ddot{y}+ \textcolor{Internal}{\mu_C}\frac{y}{r^{\,3}} = \frac{\partial {\mathpzc R}}{\partial y} = a_y\nonumber\\
\ddot{z}+ \textcolor{Internal}{\mu_C}\frac{z}{r^{\,3}} = \frac{\partial {\mathpzc R}}{\partial z} = a_z\label{eq:lpe1}
\end{align}

Different VOP formulations can be distinguished depending on the forces acting on the spacecraft, which are either conservative (e.g. third-body effect) or non-conservative (e.g. pressure, thrust, drag). For the conservative case, Lagrange planetary equations (LPE) \textcolor{Internal}{can be} employed. For the non-conservative case, Gauss' variational equations (GVE) are more adequate to compute the spacecraft's motion. These methods employ different ways of describing the perturbing effect: either using a disturbing function or a disturbing force (or acceleration). The disturbing function is the difference between perturbed and unperturbed potential functions, whereas the disturbing force expresses the specific acceleration being exerted on the spacecraft.

\paragraph{Lagrange Planetary Equations:} this method was derived to describe the osculating motion of a body subject to a conservative perturbation (i.e any gravitational effect).

The formulation of the LPE employs the first equality of Eq. \eqref{eq:lpe1}, using the disturbing function as the way to describe the perturbation. The final set of equations can be found \textcolor{Internal}{in Eq. \ref{eq:LPE} and their derivation can be followed in Battin \cite{battin}.}
\begin{align}
\nonumber\frac{da}{dt} &= \frac{2}{n a} \frac{\partial{\mathpzc{R}}}{\partial M_0}\\\nonumber
\frac{de}{dt} &= \frac{1 - e^2}{n\,a^2 e}\,\frac{\partial{\mathpzc{R}}}{\partial M_0} - \frac{\sqrt{1 - e^2}}{n\,a^2 e}\,\frac{\partial{\mathpzc{R}}}{\partial \omega}\\\nonumber
\frac{di}{dt} &= -\frac{1}{\,n a^2 \sqrt{1 - e^2} \sin i}\,\frac{\partial{\mathpzc{R}}}{\partial \Omega} - \frac{\cos i}{\,n a^2 \sqrt{1 - e^2} \sin i}\,\frac{\partial{\mathpzc{R}}}{\partial \omega}\\\nonumber
\frac{d\Omega}{dt} &= \frac{1}{n\,a^{\,2}\,\sqrt{1-e^2}\sin i}\,\frac{\partial{\mathpzc{R}}}{\partial i}\\\nonumber
\frac{d\omega}{dt} &= \frac{\sqrt{1 - e^2}}{n\,a^2 e}\,\frac{\partial{\mathpzc{R}}}{\partial e} - \frac{\cos i}{\,n a^2 \sqrt{1 - e^2} \sin i}\,\frac{\partial{\mathpzc{R}}}{\partial i}\\
\frac{dM_0}{dt} &= -\frac{2}{n\,a}\,\frac{\partial{\mathpzc{R}}}{\partial a} - \frac{1 - e^2}{\,n a^2 e}\,\frac{\partial{\mathpzc{R}}}{\partial e} \label{eq:LPE}
\end{align}
in which $a$ is the semi-major axis, $e$ is the eccentricity, $i$ is the inclination, $\Omega$ is the longitude of the ascending node, $\omega$ is the argument of the periapsis, $M$ is the mean anomaly and $n$ is the mean motion.

\paragraph{Gauss' Variational Equations:} Gauss' form of the planetary equations is perfectly equivalent to the LPE, but instead uses the second equality of Eq. \eqref{eq:lpe1}. \textcolor{Internal}{The full derivation of the equations can again be found in Battin \cite{battin}.} Since the perturbation is computed based a disturbing force, the GVE can be employed for non-conservative effects (i.e. pressure, drag or thrust). These equations can be found below:
\begin{align}
\nonumber\frac{da}{dt} &= \frac{2 a^2}{L} \Big(a_r e \sin\nu + a_{\theta}\frac{p}{r}\Big)\\\nonumber
\frac{de}{dt} &= \frac{1}{L} \Big(a_r p\sin\nu + a_{\theta}((p + r)\cos\nu + r e)\Big)\\\nonumber
\frac{di}{dt} &= a_h\frac{r \cos \upsilon}{L}\\\nonumber
\frac{d\Omega}{dt} &= a_h\frac{r \sin \upsilon}{L\sin i}\\\nonumber
\frac{d\omega}{dt} &= \frac{1}{L e} \Big(-a_r p\cos\nu + a_{\theta}(p + r)\sin\nu \Big) - a_h\frac{r \sin \upsilon \cos i}{L\sin i}\\
\frac{d\nu}{dt} &=  \frac{L}{r^2} + \frac{1}{L e} \Big(a_r p\cos\nu - a_{\theta}(p + r)\sin\nu \Big)\label{eg:gve}
\end{align}
in which \textcolor{Internal}{$a_{r}$, $a_{\theta}$ and $a_{h}$ are the acceleration components in the LVLH frame,} $L$ is the angular momentum, $p$ is the semilatus rectum, $b$ is the semi-minor axis, $r$ is the orbital position and $\upsilon = \nu + \omega$ is the argument of latitude, with $\nu$ as the true anomaly.

\section{Flow Maps}
\label{sec:FM}

In order to study of the evolution of a perturbed orbit, the computation of the motion's dynamics at each instant may not necessarily be of interest. As an alternative, one of the main techniques to study the behaviour of complex dynamical systems is the flow map: a method that maps points from their initial location at time $t_0$ to their state at time $t$ \cite{koon2008}.

The concept of flow maps is very helpful in exploring the behaviour of numerical propagations. Many techniques that fall under this category are commonly applied in astrodynamics: a classical example is the Poincar\'e map. This method is used to replace the flow of a $N$-dimensional continuous time system by a $(N-1)$-dimensional discrete time one \cite{parker}. When used in orbital propagation, it maps the state of the orbit between two consecutive crossings of the motion with a hyperplane, i.e. a Poincar\'e section. Poincar\'e maps are frequently used to study the stability of periodic orbits and can elegantly show the orbital evolution of different trajectories in the domain of interest (see e.g. \cite{villac}).

There are two basic types of Poincar\'e maps: the first one is the Poincar\'e surface of section \textit{in \textcolor{Internal}{phase} space}, which can be seen in Figure \ref{fig:poincare2}. This consists of a curve on a hyperplane, $\Sigma$, transverse to the flow in $\mathbb{R}^N$, which reflects the evolution of a trajectory point in a certain space region (for a 3-D system, $\Sigma$ is a 2-D plane). In other words, starting from some initial conditions $\bm{x}_0$, the $\Sigma$ plane contains a set of points that indicate the intersection \textcolor{External}{layout} of the orbit.

\begin{figure}[h]
	\centering
	\includegraphics[width=0.5\linewidth]{Images/poincare.pdf}
	\caption{\label{fig:poincare2}Poincar\'e surface of section with crossings $\bm{x}_0$, $\bm{x}_1$ and $\bm{x}_2$}
\end{figure}

A second version of a Poincar\'e map is obtained by a Poincar\'e surface of section \textit{in time}. For this purpose, the time-continuous dynamical system is sampled not with respect to a constraint in phase space but at discrete times $t_m = t_0 + m\cdot T,\text{ }m \in \mathbb{N}$.  In this case, the map is called a stroboscopic sampling of the phase space---the method is then named \textit{stroboscopic map}.

\textcolor{Internal}{The use of stroboscopic maps is common when computing motion in perturbed Keplerian dynamics \cite{gondelach}. The time sampling of orbital elements is often done by employing semi-analytical models \cite{roth, alessi_semi}. These are generally implemented by transforming the chosen equations of motion from being time-dependant to becoming a function of the fast angular variable (e.g. true anomaly) of the orbital elements of the massless particle. Since this fast variable changes by a fixed amount in a time period ($0$ to $2\pi$), orbital conditions can be updated at each time period by integrating over this parameter.}

\textcolor{Internal}{In order to write the equations of motion as a function of the fast variable, the former have to be multiplied by the derivative of time with respect to the fast variable $dt/dv$, where $v$ indicates the fast variable. Finally, the propagation result over one period in the fast angle $v$ returns the stroboscopic map \cite{gondelach}.}

config: arxiv