# Three-dimensional Lagrangian Coherent Structures in the Elliptic-Restricted
Three-body Problem

Tyler, Jack and Wittig, Alexander Corresponding author, PhD Candidate,
Astronautics Research Group, University of Southampton, University Road,
Southampton, SO17 1BJ<EMAIL_ADDRESS>Professor, Astronautics
Research Group, University of Southampton, University Road, Southampton, SO17
1BJ

###### Abstract

In the preliminary design of space missions it can be useful to identify
regions of dynamics that drive the system’s behaviour or separate
qualitatively different dynamics. The Lagrangian Coherent Structure (LCS) has
been widely used in the analysis of dynamical systems, and generalises the
concept of the stable and unstable manifolds to systems with arbitrary time-
dependence. However, the use of three-dimensional LCS in astrodynamics has
thus far been limited. This paper presents the application of a new numerical
method introduced by the authors, DA-LCS, to astrodynamics systems using the
Elliptic-Restricted Three-body Problem (ER3BP) as a test case. We are able to
construct the full, three-dimensional LCS associated with the Sun-Mars ER3BP
directly from the variational theory of LCS even for numerically challenging
initial conditions. The LCS is analysed in detail, showing how it in this case
separates regions of qualitatively different behaviour without any a priori
knowledge. The paper then studies the effect of integration time and the
parameterisation of the initial condition on the LCS found. We highlight how
round-off errors arise from limits of floating-point arithmetic in the most
challenging test cases and provide mitigating strategies for avoiding these
errors practically.

##  1 Introduction

In designing space missions, it can be useful to identify regions of dynamics
which drive system behaviour or separate regions of qualitatively different
dynamics. In time-independent systems, classical dynamical systems theory
identifies the invariant manifolds as separators of qualitatively different
flow in phase space, which are found by studying the system’s behaviour over
long time scales [39]. However, in systems with arbitrary time-dependence,
such structures do not always exist. Instead, one often studies the behaviour
of the system over some practical time period of interest.

Several different methods have been used to attempt to identify analogous
structures to the invariant manifolds in systems with arbitrary time-
dependence. The Finite-time Lyapunov Exponent (FTLE) is a metric that
quantifies the separation of infinitesimally-close trajectories over time
[53]. One might expect that high FTLE values signal the presence of
separatrices in the flow, but studies have shown that the FTLE does not
necessarily indicate the presence of separatrices [22] and in some cases the
hyperbolicity of local dynamics can affect the FTLE value [18]. In
astrodynamics research, the FTLE (and the closely-related Fast Lyapunov
Indicator) has been used to study chaotic astrodynamics systems and profile
behaviour since 1997 [16, 59, 26].

Since the FTLE does not necessarily highlight separatrices, and because it
requires the computationally-expensive evaluation of derivatives, alternative
heuristic measures have been investigated. The Lagrangian Descriptor (LD) has
seen use in astrodynamics to attempt to identify regions of qualitatively
different behaviour [47, 48]. The LD involves the integration of a positive,
bounded quantity along a given trajectory. While for some choices of
integrand, discontinuities in the LD scalar field can identify separatrices
between regions of qualitatively different behaviour, the selection of the
integrand can greatly change the result and must be selected for each
particular flow [35]. Moreover, being a heuristic indicator there is no
rigorous mathematical relation between the LD and any separatrices.

In 2000, Haller and Yuan [20] defined the Lagrangian Coherent Structure (LCS),
surfaces that exert significant influence on nearby trajectories. A particular
type of LCS, the hyperbolic LCS, is a surface that acts as the locally most
repelling or attracting region of flow, playing an analogous role to the
unstable and stable invariant manifolds respectively. Several equivalent
methods for determining LCS exist [19]; original definitions of the hyperbolic
LCS were made with respect to ridges of the FTLE [53], but it was later
discovered that these can only be related rigorously to LCS under certain
additional variational conditions [21].

A more recent formulation in Blazevski and Haller [7] presented a method to
determine LCS in three-dimensions from their variational theory. This method
studies the eigenvectors of the right Cauchy-Green Strain Tensor (CGST), which
quantifies the deformation of a flow over a given trajectory, and seeks
intersections of the repelling and attracting LCS with a set of hyperplanes,
deemed ‘strainlines’ and ‘stretchlines’. This approach objectively and without
any a priori knowledge determines LCS in systems with arbitrary time-
dependence, while respecting their original mathematical definition as the
locally most repelling or attracting surfaces.

The LCS has been used in a wide array of fields for profiling the behaviour of
dynamical systems [42, 51, 57, 15]. LCS derived from the FTLE field has also
been well-studied in astrodynamics. Gawlik et al. [17] investigated LCS and
their representations in the planar Elliptic-Restricted Three-body Problem
(ER3BP), and Short et al. [54] studied the use of LCS in various map
representations in the planar ER3BP. Qingyu et al. [46] extracted the ridges
of the FTLE field using a Particle Swarm optimisation algorithm in the planar
and hyperbolic restricted three-body problem. Both Short et al. [55] and
Oshima and Yanao [43] went on to analyse the behaviour of particles in the
Sun-Earth-Moon bicircular problem using regions of high FTLE as initial
conditions for strainlines and stretchlines.

Two MSc Theses, Ros Roca [50] and Parkash [44], applied the full variational
theory of LCS to inform the identification of ballistic capture orbits in the
planar ER3BP [31]. However, significant numerical difficulties and
computational complexities in the numerical method were encountered which
prevented the determination of the full LCS. This numerical difficulty comes
from the correct approximation of $C^{T}_{t_{0}}$ and the derivatives of its
eigenvalues and eigenvectors. Near the repelling LCS, maximal repulsion is
necessarily expected, and thus the CGST becomes large and potentially ill-
conditioned. Moreover, the eigenvectors of the CGST must be determined
accurately to correctly track along the structure of the LCS. The
approximation of divided differences is often used to determine the
derivatives of the flow and the derivatives of the eigenvectors and
eigenvalues of $C^{T}_{t_{0}}$, but the grid over which to approximate the
derivative must be chosen carefully, often without any a priori insight [14].
Sánchez-Martín et al. [52], who were able to compute a two-dimensional LCS in
a galactic bar model, remarked that it was ‘crucial’ to correctly choose the
grid-sizes used for approximating derivatives to obtain a valid result.

A LCS is also computationally expensive to compute, especially in three-
dimensions [14]. To reduce the computational expense, several different
approaches have been considered. Lipinski and Mohseni [28] used a ridge-
tracking algorithm that prevented the computation of trajectories away from
the LCS, instead constructing each portion of the structure iteratively. Many
different approaches to adaptive mesh refinement have also been used to
decrease run-times and improve accuracy near complex features in the flow
[40]. These methods work in tandem with standard approaches to program
parallelism and high-performance computing techniques. Lin et al. [27]
published work that details the parallelism of LCS computation using the FTLE
field on Graphics Processing Units (GPUs) using NVIDIA’s CUDA programming
paradigm.

To address both the numerical difficulties and computational expense, in Tyler
and Wittig [58] the authors recently outlined an improved numerical method,
DA-LCS, for computing LCS efficiently and automatically in three-dimensions,
without the need to manually determine or approximate derivatives. DA-LCS can
produce insight where traditional methods of approximating derivatives fails
to produce any. Differential Algebra (DA), which provides the full set of
operations to store and manipulate Taylor polynomials in a computer
environment [6, 49, 38], is used to compute derivatives of flows of the
dynamical system and the eigenvectors of $C^{T}_{t_{0}}$ accurate to machine
precision.

This paper will advance the use of LCS in astrodynamics by using this improved
numerical method to compute LCS in the Elliptic-Restricted Three-body Problem
(ER3BP), with applications to arbitrary astrodynamics systems. We will first
introduce the numerical method, and then apply it to a numerically challenging
test case to show how it can be used to accurately and automatically determine
LCS in astrodynamics systems. We analyse the LCS in detail, including
highlighting how the LCS acts in this case as a separatrix between
qualitatively different behaviours. The paper then examines how the
parameterisation of the initial condition affects the LCS found, exploiting
the fact that DA allows us to automatically compute derivatives of
parameterisations without further effort. From this investigation, we
highlight a numerical pitfall that must be considered when computing LCS
practically and provide guiding considerations on how to best choose how to
define the initial conditions in practice.

##  2 Background

###  2.1 Differential Algebra

In the following, we give a very brief introduction to Differential Algebra.
For a comprehensive treatment, the reader is referred to the literature [6].

Differential Algebra (DA) was originally introduced to determine high-order
transfer maps for particle accelerator systems [5], and provides the tools
necessary to compute the derivatives of functions within a computer
environment [8, 6]. DA constructs a Taylor series representation of an
arbitrary map, and has seen widespread use in the study of non-linearities
[34, 33, 12], the management of uncertainties [62, 63, 37, 2], and as a tool
for automatic differentiation [4, 38]. More generally, DA has also been
applied to numerous other astrodynamics problems, such as the Two-Point
Boundary Value Problem; spacecraft guidance and state estimation; trajectory
optimisation; and orbital conjunctions [29, 11, 64, 3, 9].

A full introduction to the inner workings of DA is beyond the scope of this
paper (see Berz [5]), but as an example to introduce DA, consider two real
numbers $a$ and $b\in\mathbb{R}$. The approximation to $a$ and $b$ in a
computational environment is their floating-point representation
$\bar{a},\,\bar{b}\in\mathbb{F}$, which essentially stores a set number of
digits of its binary expansion. Any operation defined in $\mathbb{R}$, $\Box$,
has a corresponding operation in $\mathbb{F},\,\boxtimes$, defined such that
the result is another floating-point approximation of the operation on the
real numbers $a$ and $b$, i.e. $\bar{a}\times\bar{b}$ commutes with the
floating-point representation of $a\times b$, $\overline{a\times b}$.

Similarly, now consider two functions, $c$ and $d$, which are sufficiently
smooth, $k-$differentiable functions of $n$ variables:
$c,\,d:\mathbb{R}^{n}\rightarrow\mathbb{R}$. In the DA framework, a computer
operates on the multivariate Taylor expansion of $c$ and $d$, $[c]$ and $[d]$,
with corresponding operations to those defined in the real function space,
such that the operation of $[c]\cdot[d]$ commutes with the DA representation
of the product $[c\cdot d]$.

Differential Algebra provides a full set of elementary operations to
efficiently operate on these multivariate expansions [61]. This includes
operations for common intrinsic functions such as division, square roots,
trigonometric functions, and exponentials, as well as operations for
differentiation and integration [63]. With these, any arbitrarily complicated
function containing these operations can be coded and evaluated in DA.

Because we are constructing a Taylor expansion to arbitrary order, and we have
access to operators for differentiation and integration, DA can be viewed as a
form of automatic differentiation. Conceptually, this provides similar
functionality as, for example, symbolic math libraries in Matlab or Python
which compute algebraic expressions that can be differentiated analytically.
But instead of representing functions as algebraic expressions, DA represents
them as the coefficients of the Taylor expansion. This makes DA
computationally more efficient, especially for complicated functions.

A key advantage of this automatic differentiation is that the derivatives of
solutions of ODEs with respect to their initial conditions can be taken
automatically, a concept known as flow expansion [9, 64, 45]. Practically,
this is achieved by replacing the floating-point arithmetic in standard
numerical ODE integration schemes by the DA arithmetic. In modern programming
languages with operator overloading such as C++, this process is
straightforward.

The Differential Algebra Computational Engine (DACE) is used to implement DA
in the program [38]. The numerical integration is performed using the 7/8th
Dormand-Prince Runge-Kutta integration as part of the Boost C++ library [1].
Custom norms are defined to connect the step-size control algorithms in the
numerical integrators in Boost and DA objects: since evaluating a vector norm
in DA yields another DA expression, for the purposes of Boost we define the
norm of a DA object to map it into the non-negative real numbers. The norm of
a DA object in this paper is the largest absolute value of any coefficient of
the expansion in any order.

###  2.2 The Cauchy-Green Strain Tensor and Lagrangian Coherent Structures

The numerical method for computing LCS using DA, summarising the main results
of Tyler and Wittig [58], is now introduced. For a more comprehensive
introduction to the mathematical formulation, the reader is referred to
Blazevski and Haller [7].

We study the behaviour of a dynamical system

$\dot{\bm{x}}=f\left(\bm{x},~{}t\right),\bm{x}\in
D\subset\mathbb{R}^{n},\,t\in\left[t_{0},\,T\right]$ (1)

where $f$ is a smooth vector field considered from time $t_{0}$ to time $T$.
Denoting a trajectory of the system starting at position $\bm{x}_{0}$ at time
$t_{0}$ up to a time $T$ as $\bm{x}\left(t_{0},\,\bm{x}_{0};\,T\right)$, the
flow map of Equation 1 is given by

$\bm{F}^{T}_{t_{0}}:\begin{cases}D\rightarrow D\\\
\bm{x}_{0}\mapsto\bm{x}\left(t_{0},\,\bm{x}_{0};\,T\right)\end{cases}$ (2)

which is assumed to be at least twice continuously differentiable. The right
Cauchy-Green Strain Tensor (CGST) $C^{T}_{t_{0}}$ is defined by the Jacobian
of this flow map, $\nabla\bm{F}^{T}_{t_{0}}$, and describes the deformation of
the flow at the end of the given trajectory

$C^{T}_{t_{0}}=\left(\nabla\bm{F}^{T}_{t_{0}}\right)^{\top}\left(\nabla\bm{F}^{T}_{t_{0}}\right)$
(3)

with ⊤ the matrix transpose. $C^{T}_{t_{0}}$ is positive-definite and
symmetric, with real eigenvalues
$0<\lambda_{1}\leq\lambda_{2}\leq\dots\leq\lambda_{n}$ and associated real
eigenvectors $\bm{\zeta}_{1},\bm{\zeta}_{2},\dots,\bm{\zeta}_{n}$. Recalling
that the full, three-dimensional LCS is the surface that is locally maximally
repelling or attractive over a time interval $\left[t_{0},\,T\right]$, this
surface is necessarily everywhere orthogonal to either $\bm{\zeta}_{n}$ or
$\bm{\zeta}_{1}$, respectively.

The LCS is constructed from its intersection with series of reference
hyperplanes $\mathcal{S}$. For repelling LCS, the intersections are termed
reduced strainlines; for attracting LCS, the intersections are reduced
stretchlines. In the following, we show the mathematical formulation for
repulsive LCS, whose structure is derived from the dominant eigenvector
$\bm{\zeta}_{n}$. A similar procedure applies to $\bm{\zeta}_{1}$ to obtain
attracting LCS.

We first sample points on each hyperplane in $\mathcal{S}$ on a uniformly-
spaced grid and compute the helicity $H_{\bm{\zeta}_{n}}$ at each point

$H_{\bm{\zeta}_{n}}=\langle\nabla\times\bm{\zeta}_{n},\,\bm{\zeta}_{n}\rangle.$
(4)

Mathematically, points that satisfy $H_{\bm{\zeta}_{n}}=0$ are maximisers of
repulsion. Numerically, this condition is relaxed such that points with
helicity below some threshold $\alpha>0$ are seed points for an ODE to
propagate the strainline forward. The ODE is tangent to the reduced strainline
and tracks within the hyperplane. In discretised form, the ODE is

$\bm{s}_{i}^{\prime}=\text{sign}\left(\bm{\zeta}_{i,\,n}\cdot\bm{\zeta}_{i-1,\,n}\right)\hat{\bm{n}}_{\mathcal{S}_{i}}\times\bm{\zeta}_{i,\,n}$
(5)

where $\bm{s}_{i}$ is the $i-$th point on the strainline and the term
$\bm{\zeta}_{i,\,n}\cdot\bm{\zeta}_{i-1,\,n}$ is introduced to enforce
continuity in the vector field by selecting the direction most closely aligned
with the previous tangent vector. The unit vector $\hat{\bm{n}}_{\mathcal{S}}$
is the normal to the hyperplane at that point. The numerical integration of
the ODE along the strainline continues until the sum of the helicity at each
$\bm{s}_{i}$ divided by the number of steps performed ($i$) rises above
$\alpha$.

The trajectories of Equation 5 are segments of strainlines forming the LCS.
However, since different initial points can belong to the same strainline, the
trajectories often overlap. They must, therefore, be filtered to provide a
single, continuous curve. If this occurs, the Fréchet distance defined in
Driemel et al. [13] is used to filter strainline segments.

This analysis is repeated for each of the hyperplanes in $\mathcal{S}$. The
strainlines forming part of the LCS on each hyperplane are then interpolated
to produce the full 3D structure of the LCS.

In Tyler and Wittig [58], we show how we use DA to compute Equations 4 and 5
to high accuracy. We first use DA to construct an expansion of $C^{T}_{t_{0}}$
with respect to the initial parameterisation using standard flow expansion
techniques [60, 2, 30, 64], accurate to machine precision. This provides
accurate determination of $\bm{\zeta}_{n}$ for use in Equation 5. Using a
novel modified power law method, another expansion of $\bm{\zeta}_{n}$ with
respect to the initial parameterisation is then constructed to obtain the seed
points in Equation 4 to high accuracy. A key advantage of using DA is that the
relevant quantities for the LCS are computed completely automatically, without
the need to adjust grid sizes, or manually derive and implement explicit
derivatives or variational equations.

##  3 Dynamical Models

###  3.1 Equations of Motion

This paper studies motion in the Sun-Mars Elliptic-Restricted Three-body
Problem (ER3BP). A generalisation of the $n-$body problem for $n=3$, the ER3BP
tracks the motion of an infinitesimally small mass $m_{3}$ under the
gravitational influence of two far larger masses $m_{1}$ (Sun) and $m_{2}$
(Mars), $m_{1}>m_{2}>>m_{3}$, which orbit each other in an ellipse. The system
is parameterised by the mass parameter $\mu=m_{2}/\left(m_{1}+m_{2}\right)$
and the eccentricity of the ellipse of the orbit of $m_{2}$ about $m_{1}$,
$e_{p}$. The special case of $e_{p}=0$ yields an autonomous system where the
classical invariant manifolds partition phase space, but for $e_{p}\neq 0$
these are not as easily defined. For the Sun-Mars case, $\mu=3.227\times
10^{-7}$ and $e_{p}=0.0935$.

The equations of motion for this system for the motion of $m_{3}$ are defined
in a rotating-pulsating reference frame centred on the barycentre of $m_{1}$
and $m_{2}$ where the $m_{1}$-$m_{2}$ distance is normalised to unity. For a
position vector $\bm{x}=\left(x,y,z\right)$ with velocity
$\bm{x}^{\prime}=\left(x^{\prime},\,y^{\prime},\,z^{\prime}\right)$, the
equations of motion are

$\displaystyle x^{\prime\prime}$ $\displaystyle=$ $\displaystyle
2y^{\prime}+\frac{\partial\Omega}{\partial x}$ (6) $\displaystyle
y^{\prime\prime}$ $\displaystyle=$
$\displaystyle-2x^{\prime}+\frac{\partial\Omega}{\partial y}$ (7)
$\displaystyle z^{\prime\prime}$ $\displaystyle=$
$\displaystyle\frac{\partial\Omega}{\partial z}$ (8)

where

$\Omega=\frac{1}{1+e_{p}\cos\nu}\left[\frac{1}{2}\left(x^{2}+y^{2}-z^{2}e\cos\nu\right)+\frac{\mu}{r_{1}}+\frac{1-\mu}{r_{2}}\right]$
(9)

and

$\displaystyle r_{1}$ $\displaystyle=$
$\displaystyle\sqrt{\left(x+\mu\right)^{2}+y^{2}+z^{2}}$ (10) $\displaystyle
r_{2}$ $\displaystyle=$
$\displaystyle\sqrt{\left(x-1+\mu\right)^{2}+y^{2}+z^{2}}.$ (11)

The variable $\nu$ is the true anomaly, the angle $m_{2}$ makes with respect
to $m_{1}$ in an inertial Cartesian coordinate system, and is used as the
independent variable. We use the notation $\dot{\Box}$ to represent
derivatives with respect to time, and $\Box^{\prime}$ to represent derivatives
with respect to the true anomaly $\nu$.

###  3.2 Definition of initial conditions

The algorithm for computing LCS introduced in Section 2.2 is explicitly
designed for a CGST that is $3\times 3$ in dimension and hence represents the
behaviour of a three-dimensional system. However, the ER3BP lives in a phase
space in $\mathbb{R}^{6}$. To simplify visualisation and enable analysis of
the ER3BP, we embed a three-dimensional sub-manifold in the six-dimensional
phase space on which we then compute the LCS. We arrive at this embedding by
first parameterising a region of interest in three dimensions in a non-
rotating frame about Mars, and then defining a transformation
$\Psi:\mathbb{R}^{3}\mapsto\mathbb{R}^{6}$ to complete the full phase space
and transform the initial condition into the ER3BP rotating-pulsating frame.
After propagation under the ER3BP equations of motion, the inverse
transformation $\Pi:\mathbb{R}^{6}\mapsto\mathbb{R}^{3}$ takes the final
condition back into the original parameterisation. This process is shown below
for an arbitrary starting parameterisation
$\left(\kappa,\beta,\,\varsigma\right)\in\mathbb{R}^{3}$:

${\left(x,\,y,\,z,\,x^{\prime},\,y^{\prime},\,z^{\prime}\right)}$${\left(x,\,y,\,z,\,x^{\prime},\,y^{\prime},\,z^{\prime}\right)_{F}}$${\left(\kappa,\,\beta,\,\varsigma\right)}$${\left(\kappa,\,\beta,\,\varsigma\right)_{F}}$$f$$\Psi\circ
f\circ\Pi$$\Psi$$\Pi$

where the subscript F implies final states. In the non-rotating frame, the
distance unit is chosen such that the Sun-Mars distance is unity; the
rationale behind this is given in Section 4.3. For more information on the
conversion between the non-rotating and rotating reference frames, the reader
is referred to Szebehely and Jefferys [56]. We refer to
$\left(\kappa,\,\beta,\,\varsigma\right)$ as the parameterisation of the
manifold and the map $\Psi\cdot f\cdot\Pi$ as the transformation.

Three pairings of $\Psi$ and $\Pi$ are explored in this paper and are
elaborated in the following. A key advantage of DA-LCS is that for all $\Psi$
and $\Pi$ given below, provided the operations are coded using DA arithmetic
the derivatives are computed fully automatically and so there is no need to
manually derive or implement derivatives. This is particularly useful when
considering that these transformations are in some cases non-trivial and
otherwise very difficult to differentiate manually.

####  3.2.1 Orbital elements

In this transformation, the initial position and velocity are given by
standard Keplerian orbital elements in the non-rotating frame, following their
use in the literature to define initial conditions for ballistic capture in
three dimensions [32].

The sub-manifold is parameterised by the radius of periapsis $r_{p}$,
inclination $i$ and argument of periapsis $\omega$ of a purely Keplerian orbit
around Mars. The hyperplanes $\mathcal{S}$ are given by fixed values of $i$
(Figure 1(a)). The full element set is obtained by fixing the remaining
orbital elements $e=0.95$, $\Omega=118^{\circ}$, $M=0$ based on values
suggested in the literature to facilitate ballistic capture [31]. This choice
of $\Omega$ maximises ‘interesting’ dynamical behaviour across all
inclinations, as the high eccentricity leads to orbits which escape into
heliocentric orbits rapidly. This presents a numerically very challenging test
case to compute the LCS for.

The full orbital element set at each point is converted to equivalent
Cartesian position and velocity in the non-rotating frame (subscript I), which
is then transformed into the ER3BP rotating-pulsating frame (subscript R) and
propagated under the full equations of motion (Equation 12).

$\displaystyle\Psi_{OE}:\begin{Bmatrix}r_{p}\\\ i\\\
\omega\end{Bmatrix}_{I}\mapsto\begin{Bmatrix}r_{p}\\\ e\\\ i\\\ \Omega\\\
\omega\\\ M\end{Bmatrix}_{I}\mapsto\begin{Bmatrix}x\\\ y\\\ z\\\ \dot{x}\\\
\dot{y}\\\ \dot{z}\end{Bmatrix}_{I}\mapsto\begin{Bmatrix}x\\\ y\\\ z\\\
x^{\prime}\\\ y^{\prime}\\\ z^{\prime}\end{Bmatrix}_{R}$ (12)
$\displaystyle\Pi_{OE}:\begin{Bmatrix}x\\\ y\\\ z\\\ x^{\prime}\\\
y^{\prime}\\\ z^{\prime}\end{Bmatrix}_{R}\mapsto\begin{Bmatrix}x\\\ y\\\ z\\\
\dot{x}\\\ \dot{y}\\\ \dot{z}\end{Bmatrix}_{I}\mapsto\begin{Bmatrix}r_{p}\\\
e\\\ i\\\ \Omega\\\ \omega\\\ M\end{Bmatrix}_{I}\mapsto\begin{Bmatrix}r_{p}\\\
i\\\ \omega\end{Bmatrix}_{I}.$ (13)

The inverse transformation $\Pi$ is as in Equation 13, and takes the final
condition in the rotating-pulsating frame of the ER3BP, transforms it back
into the non-rotating frame and converts this position and velocity into its
equivalent instantaneous orbital elements. The radius of periapsis,
inclination and argument of periapsis of the final condition are then isolated
from the full set for computing $C^{T}_{t_{0}}$.

The conversion to-and-from orbital elements and Cartesian position is
performed using standard conversion algorithms presented in Curtis [10]
modified to support DA arithmetic.

####  3.2.2 Spherical coordinates

In this transformation, the sub-manifold is parameterised using spherical
coordinates $\bm{\psi}=\left(\rho,\,\theta,\,\phi\right)$ where hyperplanes
are given by fixed $\phi$ (Figure 1(b)). This is a natural extension of the
use of polar coordinates to compute ballistic capture sets in two dimensions
in the literature [23]. The phase space is completed by uniquely associating a
velocity with each point in physical space.

This velocity is defined such that it is the velocity a point in the orbital
element transformation would have at the same Cartesian position. To ensure a
valid comparison and avoid leaking information from one set of coordinates to
the other, we need to define this velocity independently of any of the
original orbital elements. For the initial Cartesian position
$\bm{x}=\left(x,\,y,\,z\right)$ corresponding to $\bm{\psi}$ we have

$\displaystyle x=$ $\displaystyle\rho\cos\theta\sin\phi$ (14) $\displaystyle
y=$ $\displaystyle\rho\sin\theta\sin\phi$ (15) $\displaystyle z=$
$\displaystyle\rho\cos\phi$ (16)

and can define the velocity in the non-rotating frame at this point
$\bm{v}\left(\bm{x}\right)$ as

$\bm{v}\left(\bm{x}\right)=\sqrt{\text{GM}\frac{1+e}{\rho^{3}}}\left[\begin{pmatrix}x\\\
y\\\ z\end{pmatrix}\times\hat{\bm{n}}\left(\bm{x}\right)\right]$ (17)

where GM is the standard gravitational parameter of Mars in the non-rotating
frame, and $e$ has the same meaning and values as before. The quantity
$\hat{\bm{n}}$ is the normal to the orbital plane at that point and can be
computed by considering the vector that points along the line of nodes
$\hat{\bm{\Omega}}=\left(\cos\Omega,\,\sin\Omega,\,0\right)^{\top}$

$\hat{\bm{n}}\left(\bm{x}\right)=\text{sign}\left(z\right)\frac{\hat{\bm{\Omega}}\times\bm{x}}{\lvert\lvert\hat{\bm{\Omega}}\times\bm{x}\rvert\rvert}.$
(18)

The value of $\Omega$ used is the same as in the orbital elements
transformation, and the $\text{sign}\left(z\right)$ term is required to ensure
the correct orientation of the unit normal. Equation 17 is well defined
everywhere except on the line of nodes, where
$\hat{\bm{\Omega}}\times\bm{x}=\bm{0}$. Since there is no one unique
inclination passing through these points, additional information from the
orbital elements transformation would be required to attach an equivalent
velocity at these points. The reference hyperplanes should thus be selected or
sampled in such a way to avoid these ill-defined points.

Once the full phase space is complete, it is converted into the ER3BP
rotating-pulsating frame and propagated under the equations of motion

$\Psi_{S}:\begin{Bmatrix}\rho\\\ \theta\\\
\phi\end{Bmatrix}_{I}\mapsto\begin{Bmatrix}x\\\ y\\\
z\end{Bmatrix}_{I}\mapsto\begin{Bmatrix}x\\\ y\\\ z\\\ \dot{x}\\\ \dot{y}\\\
\dot{z}\end{Bmatrix}_{I}\mapsto\begin{Bmatrix}x\\\ y\\\ z\\\ x^{\prime}\\\
y^{\prime}\\\ z^{\prime}\end{Bmatrix}_{R}.\\\ $ (19)

The spatial dimensions are then converted back into spherical coordinates in
the non-rotating frame to compute $C^{T}_{t_{0}}$

$\Pi_{S}:\begin{Bmatrix}x\\\ y\\\ z\\\ x^{\prime}\\\ y^{\prime}\\\
z^{\prime}\end{Bmatrix}_{R}\mapsto\begin{Bmatrix}x\\\ y\\\ z\\\ \dot{x}\\\
\dot{y}\\\ \dot{z}\end{Bmatrix}_{I}\mapsto\begin{Bmatrix}x\\\ y\\\
z\end{Bmatrix}_{I}\mapsto\begin{Bmatrix}\rho\\\ \theta\\\
\phi\end{Bmatrix}_{I}.$ (20)

####  3.2.3 Cartesian coordinates

The final transformation examined in this paper defines an initial condition
in the non-rotating frame using Cartesian coordinates
$\bm{x}=\left(x,\,y,\,z\right)$ with hyperplanes given by fixed values of $z$
(Figure 1(c)). It is investigated here since it is conceptually simpler to
visualise than the previous two approaches, but has not yet been investigated
in the literature.

The velocity to complete the full phase space is specified in the same manner
as for spherical coordinates, following the procedure from Equation 17 onward.
Overall, this transformation is the $\Psi-\Pi$ pair

$\displaystyle\Psi_{C}:\begin{Bmatrix}x\\\ y\\\
z\end{Bmatrix}_{I}\mapsto\begin{Bmatrix}x\\\ y\\\ z\\\ \dot{x}\\\ \dot{y}\\\
\dot{z}\end{Bmatrix}_{I}\mapsto\begin{Bmatrix}x\\\ y\\\ z\\\ x^{\prime}\\\
y^{\prime}\\\ z^{\prime}\end{Bmatrix}_{R}$ (21)
$\displaystyle\Pi_{C}:\begin{Bmatrix}x\\\ y\\\ z\\\ x^{\prime}\\\
y^{\prime}\\\ z^{\prime}\end{Bmatrix}_{R}\mapsto\begin{Bmatrix}x\\\ y\\\ z\\\
\dot{x}\\\ \dot{y}\\\ \dot{z}\end{Bmatrix}_{I}\mapsto\begin{Bmatrix}x\\\ y\\\
z\end{Bmatrix}_{I}.$ (22)

(a) Graphical depiction of the parameterisation of an initial condition using
standard Keplerian elements. The light green plane represents an example
hyperplane of fixed $i$.

$\phi$
---
$\theta$
---
$\rho$
---
$\hat{x}$
---
$\hat{y}$
---
$\hat{z}$
---
(b) Graphical depiction of the spherical parameterisation of an initial
condition. The light green ‘cone’ represents an example hyperplane of fixed
$\phi$.

$\hat{x}$
---
$\hat{y}$
---
$\hat{z}$
---
$z=c$
---
$z=b$
---
$z=a$
---
(c) Graphical depiction of the Cartesian parameterisation of an initial
condition. The light green planes represents examples of different hyperplanes
of fixed $z$.

Figure 1: The parameterisations explored in this paper to define the
submanifolds. For orbital elements, the embedding is defined by fixing the
remaining orbital elements. For spherical and Cartesian coordinates, the same
unique velocity is associated with each point in space as would be obtained
from the orbital element transformation at that point.

##  4 Results

The results of the application of LCS to the Sun-Mars ER3BP are presented
here. We first analyse in detail the LCS found for the Sun-Mars ER3BP using
the orbital elements transformation and explain how we reconstruct the
surface, then present results on the effect of orbit transformation and the
choice of units that are beneficial for the user when numerically computing
LCS.

###  4.1 The LCS Structure

(a) The first step in reconstructing the full LCS structure is finding the
nearest neighbours to points on adjacent reference planes by computing the
Euclidean distance between them in parameter space.

(b) Triangles are formed between nearest neighbours in each reference plane,
which can be written directly to a STL file and displayed and manipulated as a
surface.

Figure 2: Diagrams to show how the full LCS surface is reconstructed by
interpolating between strainlines on neighbouring reference planes.

In this example, we analyse the Sun-Mars ER3BP from an initial integration
time of $t_{0}=\nu_{0}=0$ to a final integration time of $T=\nu=2\pi$ (1 full
Martian year). $\Psi$ and $\Pi$ are the orbital elements transformation
$\Psi_{OE}$ and $\Pi_{OE}$ and the helicity threshold $\alpha=10^{-9}$. The
set of reference planes $\mathcal{S}$ are defined such that

$\mathcal{S}_{\text{OE}}=\left\\{r_{p}\in\left[r,r_{s}\right],\,\omega\in\left[0,2\pi\right],i\in\left[-85^{\circ},\,-75^{\circ},\dots,85^{\circ}\right]\right\\}.$

The variables $r=1.641\times 10^{-5}$ and $r_{s}=0.00513$ here are the radius
and the Hill sphere of Mars in the scaled units of the non-rotating frame,
respectively. A $100\times 2500$ grid in $r_{p}$ and $\omega$ is used on each
plane in $\mathcal{S}$ to identify seed points. The motivation behind choosing
this grid is elaborated later.

We produce the full structure of the LCS as a STL file, which is a common file
format in 3D modelling applications. The STL file format contains a sequence
of sets of three vertices that form triangles, and the collection of faces of
the triangles are then displayed as a full structure by freely-available
viewers, such as Blender111https://blender.org or
Meshlab222https://meshlab.net.

Since this file format is relatively straightforward to create, we create the
STL file and thus the full structure using a triangulation algorithm. We first
separate the full structure into a collection of strainlines on each reference
plane in parameter space and, by leveraging that we expect there to be a
nearby strainline in the neighbouring reference plane, we find the nearest
neighbour to each point on a strainline in the neighbouring reference plane
(Figure 2(a)). This is performed by considering the Euclidean distance between
points in adjacent reference planes in the parameter space of the submanifold.

Now the two points that should be interpolated between are known, a triangle
is formed using the original point, its nearest neighbour in the next
reference plane, and a point next to the original point in the same strainline
(Figure 2(b)). This procedure is then repeated for every point in every
strainline in the reference plane, and then repeated across all reference
planes. By doing so, we are able to construct a triangulation of the full
structure of the LCS for use in visualisation software.

(a) The full structure of the LCS obtained, which is comprised of an inner
section and two large ‘arms’ which emanate from the inner section. The region
in red is the point-of-view of the section plot shown in more detail in Figure
3(c).

(b) Section view of one of the hyperplanes of fixed inclination with points
that form the strainlines shown in black to highlight more clearly how the
surface is reconstructed.

(c) Inner section of the LCS (highlighted in red in Figure 3(a)) with the
‘arms’ removed and three sample trajectories overlaid to show how it is acting
as a separatrix. The trajectory taken from inside the LCS (red) remains
orbiting Mars while the two trajectories that begin nearby but outside the
inner section (blue) both immediately escape.

(d) Trajectories in the synodic frame of points taken from different parts of
the LCS arms. Each trajectory is qualitatively different and the LCS is again
acting as a separatrix.

Figure 3: A series of representative renderings of the LCS structure
corresponding to the Sun-Mars ER3BP and reference hyperplanes outlined in the
main text. How the LCS is constructed from the strainlines is shown in more
detail, as well as how the LCS is acting a separator between qualitatively
different behaviour.

The full, interpolated repelling LCS structure from the above procedure is
given in Figure 3(a). To show more clearly how the LCS is constructed from the
strainlines on the reference hyperplanes, an example set of strainlines is
shown on a section view of the full LCS in Figure 3(b). The LCS is formed of
two distinct sections: the ‘inner’ portion that encloses Mars (Figure 3(c))
and the two extending arms that emanate from the inner structure (Figure
3(d)); these two regions are now investigated in more detail.

In Figure 3(c) the inner structure is presented with the extending arms
removed and example trajectories corresponding to some initial conditions
inside and outside of the LCS are shown. In red we show a generic trajectory
that begins inside the inner section, which continually orbits Mars over the
integration period considered. Mars lies at the centre of the inner section.
In blue are two trajectories that begin nearby but outside of the inner
section, which both immediately escape from Mars’ sphere of influence. The LCS
here is acting as separator between trajectories which have qualitatively very
different behaviour. This matches the expected structure from investigations
into ballistic capture where desired behaviour is defined beforehand [31], but
the insight here is obtained without any a priori knowledge.

A more detailed presentation of one of the two ‘arms’ is given in Figure 3(d),
again with example trajectories highlighted and the respective initial
conditions indicated by spheres. In green are two trajectories which begin
outside of the LCS structure and are qualitatively similar in immediately
escaping the sphere of influence of Mars, as was the case in Figure 3(c). More
significantly, the remaining blue and red trajectories that are taken from
different parts of the arms are also qualitatively different to both
themselves and the two green trajectories outside of the LCS, highlighting how
effectively the LCS is partitioning dynamical behaviour. In general, these
‘arms’ separate orbits which complete a different number of revolutions from
each other, including those which immediately escape from Mars. This again
matches expectations from the literature [31], which identifies these as
‘stable sets’.

###  4.2 Effect of the choice of orbit transformation and integration time

In the previous Subsection, we showed the result of DA-LCS to the Sun-Mars
ER3BP with a single choice of transformation $\Psi_{OE}$ and $\Pi_{OE}$ and a
single prescribed time period of investigation. However, there are many free
choices available to the mission designer when using LCS, such as the initial
and final times to consider, and the choices of $\Pi$ and $\Phi$. In this
Subsection we present the effect of the choice of orbit transformation and
integration time on the helicity field found. The structure of the LCS can
then be inferred by regions of consistently low helicity, which are coloured
dark blue in the plots. An advantage of DA-LCS is that it makes the
implementation of other transformations straightforward via automatic
differentiation.

(a) Orbital elements transformation.

(b) Spherical transformation.

(c) Cartesian transformation.

(d) Orbital elements transformation.

(e) Spherical transformation.

(f) Cartesian transformation.

Figure 4: The effect of integration time on the helicity field for the ER3BP.
The initial true anomaly is $t_{0}=0$ and the final true anomaly is either
$T=\pi/4$ ($68$ days, top row), or $T=2\pi$ (one Martian year, bottom row). In
general, sufficient time is required to remove the effect of the
transformations on the final derivatives. Regions of low helicity (dark blue)
imply the location of an LCS.

We compare each $\Psi-\Pi$ pair introduced in Section 3.2 for two integration
times. The integration starts at $t_{0}=\nu_{0}=0$ and ends at the true
anomaly corresponding to $T=\nu=\pi/4$ ($68$ days) and $T=2\pi$ (1 Martian
year) on the same set of initial conditions, to enable a fair comparison. The
grid of initial conditions is chosen to be those which lie on the
$i=35^{\circ}$ hyperplane in the orbital elements transformation, again
sampled on a $100\times 2500$ grid with the same parameters as used in the
construction of the full LCS. For each point on the hyperplane in the orbital
elements transformation, we also compute the equivalent Cartesian position.
This Cartesian position forms the initial condition for the Cartesian
transformation, and converting this Cartesian position to spherical
coordinates gives the equivalent initial condition in the spherical
transformation. In this way, we can ensure that each transformation samples
the same points to ensure a valid comparison. However, by choosing to
replicate the points on the $i=35^{\circ}$ plane exactly we will necessarily
pass through the line of nodes, at which point the velocity attachment for
spherical and Cartesian coordinates is ill-defined (Section 3.2.2). To avoid
this we do not evaluate points within $2^{\circ}$ of the line of nodes. When
using hyperplanes in $\phi$ or $z$, as discussed in Sections 3.2.2 and 3.2.3,
this problem does not arise.

Figure 4 presents the helicity fields for each of the transformations for both
integration times. We find that the integration time from $t_{0}$ to $T$ must
be chosen to be sufficiently large to allow the time-dependent derivatives
associated with the dynamical system in $\nabla\bm{F}^{T}_{t_{0}}$ to dominate
the time-independent derivatives associated with the transformations $\Psi$
and $\Pi$. For $T=\pi/4$, sufficient time has not passed for the derivatives
associated with the dynamics to dominate, and the result is inconsistent
between the different transformations. For $T=2\pi$, however, the results are
consistently representing the behaviour of the dynamical system. We find that
$T=\pi$ is sufficient to ensure domination of the dynamics for all points in
the domain for the examples shown here, though we note that many of the points
farthest from Mars have orbital periods of thousands of days and much lower
integration times are required for points closer to Mars.

Since the derivatives of the dynamics dominate the derivatives associated with
the transformations after sufficient time, we would expect that each
transformation yields similar helicity fields. The helicity fields signal
qualitatively the same points as LCS, identifying the arms and the inner
section as among the lowest helicity points for each transformation, despite
differences in the actual value of helicity.

Examining the quantitative agreement of the transformations from visual
inspection of Figure 4 is difficult since the width of the low-helicity
regions is only approximately $40$ km. To highlight such narrow regions, in
Figure 5 we overlay the lowest-helicity points in each field. Accurate
sampling transverse to this very narrow region to ensure we find the minimum
helicity points motivated the use of the asymmetric $100\times 2500$ grid
choice in $r_{p}$ and $\omega$, with the denser grid in $\omega$ designed to
densely sample initial conditions ‘across’ the low-helicity regions.

Figure 5 shows that the orbital elements and spherical transformations also
quantitatively agree in both signalling the ‘arms’ and inner section as the
strictly lowest regions of helicity in each case, as is to be expected.
However, there is additional structure in the spherical transformation at this
low helicity which does not agree with that predicted by the literature or the
orbital elements case. With the Cartesian transformation this additional
structure is signalled before the structure of the arms is completed. This
structure is also discontinuous, despite the helicity being a continuous
expression. By choosing a numerically challenging test case to benchmark DA-
LCS on, we have exposed a numerical error that arises as a limit of floating-
point arithmetic and manifests as ‘false positives’ in the helicity field. We
discuss this error and mitigating strategies in Section 4.3.

(a) Orbital elements transformation.

(b) Spherical coordinates transformation.

(c) Cartesian transformation.

Figure 5: Helicity fields for each of the transformations introduced
previously with the lowest points in the helicity field overlaid. The
integration time is $2\pi$ (1 Martian year.)

We also find other structure at higher helicities in the orbital elements
case, which arises as a result of the transformation also including some
information on the final velocity in the three variables that parameterise the
submanifold. Since obtaining $i$ and $\omega$ depends on the velocity,
additional information on direction of escape and direction of final velocity
is present compared to $\Pi_{S}$ and $\Pi_{C}$. In those, all final velocity
information is discarded (Equations 20 and 22). This additional structure
extends that found in ballistic capture [32, 31], but these points do not
feature in the full LCS: to be consistent with the original definition of the
LCS, the helicity threshold is chosen to find only the most repelling
structures. This yields the inner section and the ‘arms’ before any additional
structure.

###  4.3 Implementation considerations

A potential numerical pitfall exists when numerically constructing LCS in
astrodynamics problems and is elaborated in this Subsection. We also provide
guiding considerations on how to avoid these errors practically. We stress
that since DA gives derivatives accurate to floating-point precision, the
numerical difficulties given below are not as a result of DA-LCS, any
particular software package or the mathematical definition of an LCS. Rather,
they are a result of the finite precision of floating-point numbers and the
large derivatives that arise in numerically challenging test cases.

Since computers approximate the field of real numbers using floating-point
numbers, and these floating-point numbers have finite storage associated with
them, there is a limit on the precision of numbers they can represent. Double-
precision numbers have $64$ bits of storage, which provides a maximum of 16
significant digits that can be represented accurately [41, 24]. Any error
arising from this finite-precision representation compared to the evaluation
in real numbers is known as round-off error [25].

Round-off error can occur in the evaluation of the helicity and construction
of the LCS when $\nabla\bm{F}^{T}_{t_{0}}$, $C^{T}_{t_{0}}$ or
$\nabla\bm{\zeta}_{n}$ become sufficiently ill-conditioned that the ratio of
maximum to minimum absolute values of entries in those quantities approach
$10^{16}$ in any relevant orders of the expansion. At this point, when the
user requests a manipulation of these quantities, information on the smaller
numbers is lost as a result of the finite precision. This error presents as
false positives (false low values) in the helicity field.

(a) Orbital elements transformation, Mars radius.

(b) Spherical coordinates transformation, Mars radius.

(c) Cartesian coordinates transformation, Mars radius.

(d) Orbital elements transformation, Sun-Mars distance.

(e) Spherical coordinates transformation, Sun-Mars distance.

(f) Cartesian coordinates transformation, Sun-Mars distance.

Figure 6: The helicity fields presented for each of the transformations using
two unit lengths: Mars radius and Sun-Mars distance. We can aid the
conditioning of relevant terms by scaling the unit length such that
$C^{T}_{t_{0}}$ is well-conditioned in parameterisations with some bounded
terms.

We find that there are two mitigating strategies to prevent the ill-
conditioning of these terms and prevent false positives in the helicity field
in even the most challenging test cases. The first is to ensure the proper
conditioning of relevant quantities by adjusting the magnitude of the terms of
the orbit parameterisation. In Figure 6 we present the helicity fields
computed for each of the transformations with two different unit lengths: Mars
radius and Sun-Mars distance at Mars periapsis to deliberately induce this
round-off error. When using Mars radius as unit length for the orbital
elements transformation (Figure 6(a)) and the spherical coordinates
transformation (Figure 6(b)) which are both formed of a radius term and two
angles, the helicity field is ‘inverted’, with low-helicity values in
trajectories which all have qualitatively similar behaviour in immediately
escaping, and higher values in the arms and inner section. A trajectory that
immediately escapes from Mars will have a large final radius and associated
derivatives, but with the angle terms still bounded in
$\left[0,\,2\pi\right)$. The terms in $C^{T}_{t_{0}}$ associated with the
radius will thus dominate those associated with derivatives of the angles.
When the helicity at this point is then evaluated, it becomes small not
because of system dynamics, but because of a loss of precision of the
arithmetic. Small numbers that would have otherwise contributed to the
helicity have been lost to round-off error. Condition numbers $\mathcal{C}$
for $\nabla\bm{\zeta}_{n}$ of up to $10^{40}$ have been observed when using
Mars radius as the unit length with the orbital elements transformation.

However, this only occurs in regions of large final radius, i.e. points which
immediately escape. Where the final radius is lower, such as in the arms where
orbits complete several revolutions before escaping, the relevant terms are
sufficiently well-conditioned and thus helicity values are sufficiently free
of round-off errors that the ‘correct’ order of magnitude of helicity is
obtained at this point. Thus, lower helicities are present in regions of
trajectories which all have qualitatively similar behaviour (immediately
escaping) rather than in regions that separate the behaviour of trajectories,
giving rise to the apparent inversion.

We find we can aid the conditioning of relevant terms in transformations where
some of the variables that parameterise the submanifold remain bounded, as is
the case for angles in the orbital elements and spherical coordinates
transformations. For these, the distance units should be chosen such that the
entries of $C^{T}_{t_{0}}$ are all of similar magnitude. This was achieved in
this paper by using the Sun-Mars distance as unit length. Changing only the
length unit in the orbital elements transformation from Mars radius in Figure
6(a) to the Sun-Mars distance in Figure 6(d) reduces the condition number to
approximately $\mathcal{C}=10^{14}$ such that the helicity can be accurately
represented and yields the expected result. We stress that these numerical
effects are not specific to DA-LCS, any particular software package, or the
mathematical definition of LCS. Rather, it is inherent to floating-point
arithmetic and has been highlighted by the choice of a numerically
particularly challenging test case.

Scaling units to ensure the Jacobian is well-conditioned is more difficult in
systems where all variables that parameterise the submanifold are unbounded,
such as the Cartesian transformation. Here, where trajectories typically
escape in the direction of $x$ and $y$ but stay close to the plane in $z$,
terms in $z$ can be dominated. Scaling the length unit uniformly will affect
all terms equally and not change the condition number, which is above
$10^{16}$ after even $68$ days. This causes the poor numerical resolution in
the Cartesian transformation observed earlier: scaling all terms of the
parameterisation does not affect round-off error in regions of escape, causing
the additional structure in trajectories that immediately escape to persist.

The second mitigating strategy is to use a transformation that is naturally
more robust to ill-conditioning. Even when using Sun-Mars distance as the unit
length in spherical coordinates (Figure 6(e)) there remains sufficient ill-
conditioning in the trajectories which escape ($\mathcal{C}\sim 10^{17}$) that
false positives are present at the same helicity as the arms and inner
section. However, for the orbital elements transformation the helicity can be
evaluated accurately at all points in the field. We remark that we have seen
satisfactory performance from spherical coordinates in less challenging test
cases (see Tyler and Wittig [58]) with values of $\Omega$ or $e$ that yield
trajectories which escape less quickly. The numerical stability of orbital
elements has been noted by other authors in similar situations (for example,
Martin and Schaub [36]).

We therefore recommend against the use of the Cartesian transformation, except
in isotropic systems or where one scales each term in the parameterisation
non-uniformly using different unit lengths. In challenging test cases, we
recommend the use of orbital elements to define initial conditions, which have
produced usable insight even for the challenging problem parameters selected
here. They also provide a more natural expression of defining initial
conditions in the context of mission design. Nonetheless, other
transformations may still yield better results depending on the system being
studied.

We again stress that since DA provides derivatives accurate to machine
precision, the above numerical difficulties are a result of double precision
floating-point arithmetic. The use of quadruple precision (allowing up to 36
significant digits) would prevent round-off errors until approaching condition
numbers of $10^{36}$. However, quadruple precision is not currently
implemented natively on most CPUs, and software solutions are generally
computationally very intensive and do not necessarily guarantee
reproducibility between platforms [24, 41].

##  5 Conclusion

We have presented the the results of the application of DA-LCS, a new
numerical method suggested by the authors for the computation of Lagrangian
Coherent Structures in three-dimensions. A full three-dimensional Lagrangian
Coherent Structure in the Sun-Mars ER3BP has been computed accurately with DA-
LCS, even in a numerically challenging test case designed to benchmark the
method against. Detailed analysis of the LCS highlighted how it in this case
acts as a separator between regions of qualitatively different behaviour,
generalising the concept of stable and unstable manifolds. We discovered that
sufficient integration time is required to ensure the LCS found reflects the
system dynamics, and that care must be taken in defining initial conditions to
avoid regions of particularly poor numerical numerical performance due to
limitations in floating-point arithmetic. While the LCS was determined here
for ER3BP, there is no obstacle to computing it in other systems provided the
arithmetic can be readily replaced by DA operations; high-fidelity ephemeris
models have been tested with DA-LCS and shown to yield similar insight as for
the ER3BP case studied here.

We have several suggestions for future work. Firstly, we sampled the helicity
on a $100\times 2500$ grid to find seed points for the strainlines. However,
the major advantage of a LCS is that strainlines can be propagated forward
once a single seed point is known, avoiding the computation of trajectories
that do not otherwise correspond to the structure of the LCS. It would be a
worthwhile pursuit to examine whether more sophisticated grid search
strategies would reduce computation time, potentially exploiting a DA
expansion of the helicity in combination with a gradient descent method to
identify minima through local optimisation rather than grid search. Moreover,
the extension of LCS to $n-$dimensions would avoid the need to define the
dimension reductions introduced in this paper, and is numerically more simple
to undertake than standard approaches if one follows the same methods used in
DA-LCS. However, the interpretation and visualisation of $n-$dimensional
structures remains a challenge, and so significant effort would be required to
store, manipulate and visualise these structures.

To support this work and advance the use of LCS in astrodynamics, the source
code for DA-LCS and the program used to produce this study will be made freely
available under a GNU GPLv3 license.

## Acknowledgments

The authors would like to acknowledge financial support from the EPSRC Centre
for Doctoral Training in Next Generation Computational Modelling grant
EP/L015382/1, and the use of the IRIDIS High Performance Computing Facility
and associated support services at the University of Southampton.

## References

  * Ahnert and Mulansky [2011] Karsten Ahnert and Mario Mulansky.  Odeint – Solving Ordinary Differential Equations in C++.  _AIP Conference Proceedings_ , 1389(1):1586–1589, September 2011.  ISSN 0094-243X.  doi: 10.1063/1.3637934.
  * Armellin et al. [2010] Roberto Armellin, Pierluigi Di Lizia, Franco‘ Bernelli-Zazzera, and Martin Berz.  Asteroid close encounters characterization using differential algebra: The case of Apophis.  _Celestial Mechanics and Dynamical Astronomy_ , 107(4):451–470, 2010.  ISSN 09232958.  doi: 10.1007/s10569-010-9283-5.
  * Armellin et al. [2012] Roberto Armellin, Pierluigi Di Lizia, Alessandro Morselli, and Michelle Lavagna.  An orbital conjunction algorithm based on Taylor models.  In _Advances in the Astronautical Sciences_ , volume 143, pages 2127–2144, 2012.  ISBN 9780877035817.
  * Bartholomew-Biggs et al. [2000] Michael Bartholomew-Biggs, Steven Brown, Bruce Christianson, and Laurence Dixon.  Automatic differentiation of algorithms.  _Journal of Computational and Applied Mathematics_ , 124(1):171–190, December 2000.  ISSN 0377-0427.  doi: 10.1016/S0377-0427(00)00422-2.  URL https://www.sciencedirect.com/science/article/pii/S0377042700004222.
  * Berz [1987] Martin Berz.  The method of power series tracking for the mathematical description of beam dynamics.  _Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment_ , 258(3):431–436, 8 1987.  ISSN 0168-9002.  doi: 10.1016/0168-9002(87)90927-2.  URL https://www.sciencedirect.com/science/article/pii/0168900287909272.
  * Berz [1999] Martin Berz.  _Modern Map Methods in Particle Beam Physics_ , volume 108.  Academic Press, 1999.  ISBN 0120147505.  URL http://bt.pa.msu.edu/cgi-bin/display.pl?name=AIEP108book.
  * Blazevski and Haller [2014] Daniel Blazevski and George Haller.  Hyperbolic and elliptic transport barriers in three-dimensional unsteady flows.  _Physica D: Nonlinear Phenomena_ , 273-274:46–62, 2014\.  ISSN 01672789.  doi: 10.1016/j.physd.2014.01.007.  URL http://dx.doi.org/10.1016/j.physd.2014.01.007.
  * Cavenago et al. [2017] Francesco Cavenago, Pierluigi Di Lizia, Mauro Massari, and Alexander Wittig.  On-board DA-based state estimation algorithm for spacecraft relative navigation.  _7th European conference for aeronautics and space sciences (EUCASS)_ , pages 1–14, 2017.  doi: 10.13009/EUCASS2017-607.  URL https://www.eucass.eu/doi/EUCASS2017-607.pdf.
  * Cavenago et al. [2018] Francesco Cavenago, Pierluigi Di Lizia, Mauro Massari, Simone Servadio, and Alexander Wittig.  DA-based nonlinear filters for spacecraft relative state estimation.  In _2018 Space Flight Mechanics Meeting, 8-12 January 2018_ , number February. American Institute of Aeronautics and Astronautics, 2018.  doi: 10.2514/6.2018-1964.
  * Curtis [2010] Howard D Curtis.  _Orbital mechanics for engineering students_.  Elsevier, Oxford, 2nd editio edition, 2010.  ISBN 9780123747785 (hbk.)\r0123747783 (hbk.)\r9781856179546 (hbk. with on line testing)\r1856179540 (hbk. with on line testing).
  * Di Lizia et al. [2012] P. Di Lizia, R. Armellin, A. Ercoli-Finzi, and M. Berz.  High-order robust guidance of interplanetary trajectories based on differential algebra.  _Journal of Aerospace Engineering, Sciences and Applications_ , 1(1):43–57, 2012.  ISSN 2236577X.  doi: 10.7446/jaesa.0101.05.
  * Di Mauro et al. [2015] Guiseppe Di Mauro, Markus Schlotterer, Stephan Theil, and Michelle Lavagna.  Nonlinear Control for Proximity Operations Based on Differential Algebra.  _Journal of Guidance, Control, and Dynamics_ , 38(11):2173–2187, 2015.  ISSN 0731-5090.  doi: 10.2514/1.g000842.
  * Driemel et al. [2016] Anne Driemel, Amer Krivosija, and Christian Sohler.  Clustering time series under the Fréchet distance.  _Proceedings of the Annual ACM-SIAM Symposium on Discrete Algorithms_ , 2:766–785, 2016.  doi: 10.1137/1.9781611974331.ch55.
  * Farazmand and Haller [2012] Mohammad Farazmand and George Haller.  Computing Lagrangian coherent structures from their variational theory.  _Chaos_ , 22(1), 2012.  ISSN 10541500.  doi: 10.1063/1.3690153.
  * Fiorentino et al. [2012] L. A. Fiorentino, M. J. Olascoaga, A. Reniers, Z. Feng, F. J. Beron-Vera, and J. H. MacMahan.  Using Lagrangian Coherent Structures to understand coastal water quality.  _Continental Shelf Research_ , 47:145–149, 2012.  ISSN 02784343.  doi: 10.1016/j.csr.2012.07.009.  URL http://dx.doi.org/10.1016/j.csr.2012.07.009.
  * Froeschlé et al. [1997] Cl Froeschlé, R. Gonczi, and E. Lega.  The fast Lyapunov indicator : A simple tool to detect weak chaos. Application to the structure of the main asteroidal belt.  _Planetary and Space Science_ , 45(7 SPEC. ISS.):881–886, 1997.  ISSN 00320633.  doi: 10.1016/s0032-0633(97)00058-5.
  * Gawlik et al. [2009] Evan S. Gawlik, Jerrold E. Marsden, Philip C. Du Toit, and Stefano Campagnola.  Lagrangian coherent structures in the planar elliptic restricted three-body problem.  _Celestial Mechanics and Dynamical Astronomy_ , 103(3):227–249, 2009.  ISSN 09232958.  doi: 10.1007/s10569-008-9180-3.
  * Gondelach et al. [2019] David J. Gondelach, Roberto Armellin, and Alexander Wittig.  On the predictability and robustness of Galileo disposal orbits.  _Celestial Mechanics and Dynamical Astronomy_ , 131(12):1–30, 2019.  ISSN 15729478.  doi: 10.1007/s10569-019-9938-9.  URL https://doi.org/10.1007/s10569-019-9938-9.
  * Hadjighasem et al. [2017] Alireza Hadjighasem, Mohammad Farazmand, Daniel Blazevski, Gary Froyland, and George Haller.  A Critical Comparison of Lagrangian Methods for Coherent Structure Detection.  _Chaos_ , 27(5):053104, 2017.  doi: 10.1063/1.4982720.
  * Haller and Yuan [2000] G. Haller and G. Yuan.  Lagrangian coherent structures and mixing in two-dimensional turbulence.  _Physica D: Nonlinear Phenomena_ , 147(3-4):352–370, 2000.  ISSN 01672789.  doi: 10.1016/S0167-2789(00)00142-1.
  * Haller [2011] George Haller.  A variational theory of hyperbolic Lagrangian Coherent Structures.  _Physica D: Nonlinear Phenomena_ , 240(7):574–598, 3 2011.  ISSN 01672789.  doi: 10.1016/j.physd.2010.11.010.  URL https://www.sciencedirect.com/science/article/pii/S0167278910003143.
  * Haller and Sapsis [2011] George Haller and Themistoklis Sapsis.  Lagrangian coherent structures and the smallest finite-time Lyapunov exponent.  _Chaos_ , 21(2):1–7, 2011.  ISSN 10541500.  doi: 10.1063/1.3579597.
  * Hyeraci and Topputo [2010] Nicola Hyeraci and Francesco Topputo.  Method to design ballistic capture in the elliptic restricted three-body problem.  _Journal of Guidance, Control, and Dynamics_ , 33(6):1814–1823, 2010.  ISSN 15333884.  doi: 10.2514/1.49263.
  * Institute for Electronics and Electrical Engineers [2019] Institute for Electronics and Electrical Engineers.  754-2019 - IEEE Standard for Floating-Point Arithmetic.  Technical report, IEEE, 2019.  ISBN: 9781504459242.
  * Lantz [1971] R.B. Lantz.  Quantitative Evaluation of Numerical Diffusion (Truncation Error).  _Society of Petroleum Engineers Journal_ , 11(03):315–320, September 1971.  ISSN 0197-7520.  doi: 10.2118/2811-PA.  URL https://doi.org/10.2118/2811-PA.
  * Lega et al. [2016] Elena Lega, Massimiliano Guzzo, and Claude Froeschlé.  Theory and Applications of the Fast Lyapunov Indicator (FLI) Method.  In _Chaos Detection and Predictability_ , volume 915, pages 35–54. Springer, 2016.  doi: 10.1007/978-3-662-48410-4–“˙˝2.
  * Lin et al. [2017] Mingpei Lin, Ming Xu, and Xiaoyu Fu.  GPU-accelerated computing for Lagrangian coherent structures of multi-body gravitational regimes.  _Astrophysics and Space Science_ , 362(4), 2017.  ISSN 1572946X.  doi: 10.1007/s10509-017-3050-y.  URL http://dx.doi.org/10.1007/s10509-017-3050-y.
  * Lipinski and Mohseni [2010] Doug Lipinski and Kamran Mohseni.  A ridge tracking algorithm and error estimate for efficient computation of Lagrangian coherent structures.  _Chaos_ , 20(1), 2010.  ISSN 10541500.  doi: 10.1063/1.3270049.
  * Lizia et al. [2008] P. Di Lizia, R. Armellin, and M. Lavagna.  Application of high order expansions of two-point boundary value problems to astrodynamics.  _Celestial Mechanics and Dynamical Astronomy_ , 102(4):355–375, 12 2008.  ISSN 0923-2958.  doi: 10.1007/s10569-008-9170-5.  URL http://link.springer.com/10.1007/s10569-008-9170-5.
  * Lunghi et al. [2018] Paolo Lunghi, Michèle Lavagna, and Pierluigi Di Lizia.  Atmospheric entry guidance based on Differential Algebra for high elevation Mars landing.  In _2018 Space Flight Mechanics Meeting, 8-12 January_ , number January. American Institute of Aeronautics and Astronautics, 2018.  doi: 10.2514/6.2018-1458.
  * Luo and Topputo [2015] Z. F. Luo and F. Topputo.  Analysis of ballistic capture in Sun-planet models.  _Advances in Space Research_ , 2015.  ISSN 18791948.  doi: 10.1016/j.asr.2015.05.042.
  * Luo et al. [2014] Z. F. Luo, F. Topputo, F. Bernelli-Zazzera, and G. J. Tang.  Constructing ballistic capture orbits in the real Solar System model.  _Celestial Mechanics and Dynamical Astronomy_ , 120(4):433–450, 2014.  ISSN 15729478.  doi: 10.1007/s10569-014-9580-5.
  * Makino [1998] Kyoko Makino.  _Rigorous Analysis of Nonlinear Motion in Particle Accelerators_.  PhD thesis, Michigan State University, 1998.  URL http://www.bt.pa.msu.edu/pub/papers/makinophd/makinophd.ps.
  * Makino and Berz [1996] Kyoko Makino and Martin Berz.  Remainder Differential Algebras and their Applications.  _Computational Differentiation: Techniques, Applications, and Tools_ , pages 63–74, 1996.
  * Mancho et al. [2013] Ana M. Mancho, Stephen Wiggins, Jezabel Curbelo, and Carolina Mendoza.  Lagrangian Descriptors: A Method for Revealing Phase Space Structures of General Time Dependent Dynamical Systems.  _Communications in Nonlinear Science and Numerical Simulation_ , 18(12):3530–3557, December 2013.  ISSN 10075704.  doi: 10.1016/j.cnsns.2013.05.002.  URL http://arxiv.org/abs/1106.1306.  arXiv:1106.1306 [nlin, physics:physics].
  * Martin and Schaub [2022] John R Martin and Hanspeter Schaub.  Periodic Orbit Discovery Enhanced by Physics-informed Neural Networks.  In _2022 Astrodynamics Specialist Conference, Charlotte, North Carolina, 7-11 August 2022_. American Institute of Aeronautics and Astronautics, August 2022.  In press.
  * Massari et al. [2017] Mauro Massari, Pierluigi Di Lizia, and Mirco Rasotto.  Nonlinear Uncertainty Propagation in Astrodynamics Using Differential Algebra and Graphics Processing Units.  _Journal of Aerospace Information Systems_ , 14(9):493–503, 2017.  doi: 10.2514/1.i010535.
  * Massari et al. [2018] Mauro Massari, Pierluigi Di Lizia, Francesco Cavenago, and Alexander Wittig.  Differential Algebra software library with automatic code generation for space embedded applications.  In _2018 Space Flight Mechanics Meeting, 8-12 January_ , number January. American Institute of Aeronautics and Astronautics, 2018.  doi: 10.2514/6.2018-0398.
  * Meiss [1992] J. D. Meiss.  Symplectic maps, variational principles, and transport.  _Reviews of Modern Physics_ , 64(3):795–848, 1992\.  ISSN 00346861.  doi: 10.1103/RevModPhys.64.795.
  * Miron et al. [2012] Philippe Miron, Jérôme Vétel, André Garon, Michel Delfour, and Mouhammad El.  Anisotropic mesh adaptation on Lagrangian Coherent Structures.  _Journal of Computational Physics_ , 231(19):6419–6437, 2012.  ISSN 0021-9991.  doi: 10.1016/j.jcp.2012.06.015.  URL http://dx.doi.org/10.1016/j.jcp.2012.06.015.
  * Muller et al. [2010] Jean-Michel Muller, Nicolas Brisebarre, Florent de Dinechin, Claude-Pierre Jeannerod, Vincent Lefèvre, Guillaume Melquiond, Nathalie Revol, Damien Stehlé, and Serge Torres.  _Handbook of Floating-Point Arithmetic_.  Birkhäuser Boston, 2010.  doi: 10.1007/978-0-8176-4705-6.  Publication Title: Handbook of Floating-Point Arithmetic.
  * Nolan et al. [2020] Peter J. Nolan, Hosein Foroutan, and Shane D. Ross.  Pollution transport patterns obtained through generalized lagrangian coherent structures.  _Atmosphere_ , 11(2):1–14, 2020.  ISSN 20734433.  doi: 10.3390/atmos11020168.
  * Oshima and Yanao [2014] Kenta Oshima and Tomohiro Yanao.  Applications of gravity assists in the bicircular and bielliptic restricted four-body problem.  In _Advances in the Astronautical Sciences_ , volume 152, pages 503–522, 2014.
  * Parkash [2019] A. S. Parkash.  Application of Lagrangian Coherent Structures to the computation and understanding of ballistic capture trajectories.  Master’s thesis, Delft University of Technology, 2019.
  * Pérez-Palau et al. [2015] Daniel Pérez-Palau, Josep J. Masdemont, and Gerard Gómez.  Tools to detect structures in dynamical systems using Jet Transport.  _Celestial Mechanics and Dynamical Astronomy_ , 123(3):239–262, 2015.  ISSN 15729478.  doi: 10.1007/s10569-015-9634-3.
  * Qingyu et al. [2020] Qu Qingyu, Lin Mingpei, and Xu Ming.  Lagrangian Coherent Structures in the Planar Parabolic/Hyperbolic Restricted Three-Body Problem.  _Monthly Notices of the Royal Astronomical Society_ , 2020.  ISSN 0035-8711.  doi: 10.1093/mnras/staa199.
  * Quinci et al. [2022] Alessio Quinci, Gianmario Merisio, and Francesco Topputo.  Qualitative study of ballistic capture at mars via lagrangian descriptors.  Submitted, under review., 2022.
  * Raffa et al. [2022] Sebastiano Raffa, Gianmario Merisio, and Francesco Topputo.  Finding regions of bounded motion in binary asteroid environment using lagrangian descriptors.  Submitted, under review., 2022.
  * Rasotto et al. [2016] M. Rasotto, A. Moreselli, A Wittig, M. Massari, P. Di Lizia, R. Armellin, C. Y. Valles, and G. Ortega.  Differential Algebra Space Toolbox for Nonlinear Uncertainty Propagation in Space Dynamics.  In _International Conference on Astrodynamics Tools and Techniques_ , number 6, 2016.
  * Ros Roca [2015] Xavier Ros Roca.  Computation of Lagrangian Coherent Structures with Application to Weak Stability Boundaries.  Master’s thesis, Politecnico di Milano, 2015.
  * Rutherford et al. [2012] B. Rutherford, G. Dangelmayr, and M. T. Montgomery.  Lagrangian coherent structures in tropical cyclone intensification.  _Atmospheric Chemistry and Physics_ , 12(12):5483–5507, 2012.  ISSN 16807316.  doi: 10.5194/acp-12-5483-2012.
  * Sánchez-Martín et al. [2018] P. Sánchez-Martín, J. J. Masdemont, and M. Romero-Gómez.  From manifolds to Lagrangian coherent structures in galactic bar models.  _Astronomy & Astrophysics_, 618:A72, 2018.  ISSN 0004-6361.  doi: 10.1051/0004-6361/201833451.
  * Shadden et al. [2005] Shawn C. Shadden, Francois Lekien, and Jerrold E. Marsden.  Definition and properties of Lagrangian coherent structures from finite-time Lyapunov exponents in two-dimensional aperiodic flows.  _Physica D: Nonlinear Phenomena_ , 212(3-4):271–304, 2005.  ISSN 01672789.  doi: 10.1016/j.physd.2005.10.007.
  * Short et al. [2011] Cody R. Short, Kathleen C. Howell, and Xavier M. Tricoche.  Lagrangian coherent structures in the restricted three-body problem.  _Advances in the Astronautical Sciences_ , 140:2111–2126, 2011.  ISSN 00653438.
  * Short et al. [2015] Cody R. Short, Daniel Blazevski, Kathleen C. Howell, and George Haller.  Stretching in phase space and applications in general nonautonomous multi-body problems.  _Celestial Mechanics and Dynamical Astronomy_ , 122(3):213–238, 2015.  ISSN 15729478.  doi: 10.1007/s10569-015-9617-4.
  * Szebehely and Jefferys [1968] Victor Szebehely and William H. Jefferys.  Theory of Orbits: The Restricted Problem of Three Bodies.  _American Journal of Physics_ , 36(4):375–375, 1968.  ISSN 0002-9505.  doi: 10.1119/1.1974535.
  * Tu et al. [2019] Han Tu, Matthew Marzanek, Melissa A. Green, and David E. Rival.  Investigation of accelerating non-slender delta-wing planforms at high angle of attack using lagrangian coherent structures.  In _AIAA Scitech 2019 Forum_ , pages 1–10. American Institute of Aeronautics and Astronautics, 2019.  ISBN 9781624105784.  doi: 10.2514/6.2019-2165.
  * Tyler and Wittig [2022] Jack Tyler and Alexander Wittig.  An Improved Numerical Method for Three-dimensional Hyperbolic Lagrangian Coherent Structures using Differential Algebra.  _Submitted, under review_ , 2022.  URL http://arxiv.org/abs/2204.06236.  arXiv: 2204.06236.
  * Villac [2008] B. F. Villac.  Using FLI maps for preliminary spacecraft trajectory design in multi-body environments.  _Celestial Mechanics and Dynamical Astronomy_ , 102(1-3):29–48, 2008.  ISSN 09232958.  doi: 10.1007/s10569-008-9158-1.
  * Wittig [2012] Alexander Wittig.  _Rigorous High-Precision Enclosures of Fixed Points and their Invariant Manifolds_.  PhD thesis, Michigan State University, 2012.
  * Wittig et al. [2010] Alexander Wittig, Martin Berz, Johnathan Grote, Kyoko Makino, and Sheldon Newhouse.  Rigorous and accurate enclosure of invariant manifolds on surfaces.  _Regular and Chaotic Dynamics_ , 15(2):107–126, 2010.  ISSN 15603547.  doi: 10.1134/S1560354710020024.
  * Wittig et al. [2014] Alexander Wittig, Pierluigi Di Lizia, Roberto Armellin, Franco Bernelli Zazzera, Kyoko Makino, and Martin Berzş.  An automatic domain splitting technique to propagate uncertainties in highly nonlinear orbital dynamics.  _Advances in the Astronautical Sciences_ , 152:1923–1941, 2014.  ISSN 00653438.
  * Wittig et al. [2015] Alexander Wittig, Pierluigi Di Lizia, Roberto Armellin, Kyoko Makino, Franco Bernelli-Zazzera, and Martin Berz.  Propagation of large uncertainty sets in orbital dynamics by automatic domain splitting.  _Celestial Mechanics and Dynamical Astronomy_ , 122(3):239–261, 7 2015.  ISSN 0923-2958.  doi: 10.1007/s10569-015-9618-3.  URL http://link.springer.com/10.1007/s10569-015-9618-3.
  * Zazzera et al. [2012] Franco Bernelli Zazzera, M. Lavagna, Roberto Armellin, Pierluigi Di Lizia, A. Moreselli, J. Olympio, and Dario Izzo.  Trajectory Optimisation under Uncertainties.  Technical report, European Space Agency, 2012.