# Relativity of superluminal observers in $1+3$ spacetime

Andrzej Dragan<EMAIL_ADDRESS>Institute of Theoretical Physics, University
of Warsaw, Pasteura 5, 02-093 Warsaw, Poland  Centre for Quantum Technologies,
National University of Singapore, 3 Science Drive 2, 117543 Singapore,
Singapore  Kacper Dębski  Institute of Theoretical Physics, University of
Warsaw, Pasteura 5, 02-093 Warsaw, Poland  Szymon Charzyński  Department of
Mathematical Methods in Physics, University of Warsaw, Pasteura 5, 02-093
Warsaw, Poland  Krzysztof Turzyński  Institute of Theoretical Physics,
University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland  Artur Ekert  Centre
for Quantum Technologies, National University of Singapore, 3 Science Drive 2,
117543 Singapore, Singapore  Mathematical Institute, University of Oxford,
Woodstock Road, Oxford OX2 6GG, United Kingdom

###### Abstract

We develop an extension of special relativity in $1+3$ dimensional spacetime
to account for superluminal inertial observers and show that such an extension
rules out the conventional dynamics of mechanical point-like particles and
forces one to use a field-theoretic framework. Therefore we show that field
theory can be viewed as a direct consequence of extended special relativity.

##  I Introduction

In a recent paper Dragan and Ekert showed that extending special relativity to
account for superluminal particles and observers does not result in causal
paradoxes [1], as was commonly believed [2]. Instead, such an extension
modifies the notion of causality in the same way quantum theory does. In
particular, it was shown that when Galilean relativity (involving all inertial
observers) is assumed in $1+1$ dimensional spacetime, indeterministic behavior
and motion along multiple paths analogous to quantum mechanical superpositions
becomes inevitable. This result concluded the debate, initiated by Tolman [2],
about alleged problems with causality triggered by the presence of
superluminal particles in relativity [3, 4, 5, 6, 7, 8, 9, 10, 11].

However, generalizing this scheme to a relativistic framework in a $1+3$
dimensional spacetime poses some serious challenges, both mathematical and
interpretational [12]. The main obstacle, pointed out by Marchildon, Antippa
and Everett [13], results from the fact that the smallest group containing
superluminal transformations is $SL(4,\mathbb{R})$. This group, however,
cannot be a symmetry group, because it contains elements such as direction-
dependent dilations, which are not observed. In this work, we propose a way to
overcome these difficulties and consistently extend special relativity in
$1+3$ dimensional spacetime to superluminal frames of reference. We argue that
within such an extension the standard classical dynamics of a point-like
particle cannot be supported and the only relativistically invariant dynamics
of any physical systems involving both subluminal and superluminal observers
requires field-theoretic framework.

##  II Kinematics

An orthodox $1+3$ dimensional subluminal Lorentz transformation between two
mutually unrotated frames $(ct,\bm{r})$ and $(ct^{\prime},\bm{r^{\prime}})$
moving with a relative velocity $\bm{V}$ leaves the transversal spacial
components of four-position unaffected:
$\bm{r^{\prime}}-\frac{\bm{r^{\prime}}\cdot\bm{V}}{V^{2}}\bm{V}=\bm{r}-\frac{\bm{r}\cdot\bm{V}}{V^{2}}\bm{V}$.
The longitudinal component undergoes the standard Lorentz transformation:
$\frac{\bm{r^{\prime}}\cdot\bm{V}}{V^{2}}\bm{V}=\frac{\frac{\bm{r}\cdot\bm{V}}{V^{2}}-t}{\sqrt{1-\frac{V^{2}}{c^{2}}}}\bm{V}$.
Adding these two equations by sides and supplementing them with the Lorentz
transformation for the temporal coordinate $t^{\prime}$ leads to the $1+3$
dimensional Lorentz transformation between a pair of mutually unrotated
observers [14]:

$\displaystyle\bm{r^{\prime}}$ $\displaystyle=$
$\displaystyle\bm{r}-\frac{\bm{r}\cdot\bm{V}}{V^{2}}\bm{V}+\frac{\frac{\bm{r}\cdot\bm{V}}{V^{2}}-t}{\sqrt{1-\frac{V^{2}}{c^{2}}}}\bm{V},$
$\displaystyle ct^{\prime}$ $\displaystyle=$
$\displaystyle\frac{ct-\frac{\bm{r}\cdot\bm{V}}{c}}{\sqrt{1-\frac{V^{2}}{c^{2}}}}.$
(1)

Such a transformation preserves the metric, so that the form of the spacetime
interval remains the same in both frames:
$c^{2}\text{d}t^{2}-\text{d}\bm{r}\cdot\text{d}\bm{r}=c^{2}\text{d}t^{\prime
2}-\text{d}\bm{r^{\prime}}\cdot\text{d}\bm{r^{\prime}}$. The inverse
transformation can be obtained by substituting $\bm{V}\to-\bm{V}$, as well as
$\bm{r}\leftrightarrow\bm{r^{\prime}}$ and $t\leftrightarrow t^{\prime}$,
instead of algebraically reversing the equations. A velocity transformation
formula is obtained by evaluating the derivative
$\bm{v^{\prime}}\equiv\frac{\text{d}\bm{r^{\prime}}}{\text{d}t^{\prime}}$:

$\displaystyle\bm{v^{\prime}}=\frac{\sqrt{1-\frac{V^{2}}{c^{2}}}\left(\bm{v}-\frac{\bm{v}\cdot\bm{V}}{V^{2}}\bm{V}\right)-\left(\bm{V}-\frac{\bm{v}\cdot\bm{V}}{V^{2}}\bm{V}\right)}{1-\frac{\bm{v}\cdot\bm{V}}{c^{2}}}.$
(2)

Let us now carry out a similar procedure for the superluminal transformations
with $V>c$ discussed in [1]. Superluminal observers have been discussed since
the 1960s [15], but even the simplest $1+1$ dimensional transformations
considered in literature were missing the anti-symmetric factor [16, 17]
rendering the whole theory non-covariant. The correct expressions containing
that anti-symmetric term appeared for the first time in a paper by Parker
[12]:

$\displaystyle x^{\prime}$
$\displaystyle=\pm\frac{V}{|V|}\frac{x-Vt}{\sqrt{\frac{V^{2}}{c^{2}}-1}},$
$\displaystyle t^{\prime}$
$\displaystyle=\pm\frac{V}{|V|}\frac{t-\frac{V}{c^{2}}x}{\sqrt{\frac{V^{2}}{c^{2}}-1}}.$
(3)

This transformation changes the sign of the metric:
$c^{2}\text{d}t^{2}-\text{d}x^{2}=-c^{2}\text{d}t^{\prime 2}+\text{d}x^{\prime
2}$, but preserves null intervals, so the speed of light is the same in both
subluminal and superluminal frames of reference, as required. The unspecified
sign $\pm$ cannot be uniquely determined, as no $V\to 0$ exists, but, as
discussed in [1], we will pick the negative sign convention so that the
transformation (II) remains a hyperbolic rotation.

In order to generalize (II) to the $1+3$ dimensional spacetime we introduce an
additional pair of spatial dimensions $y$ and $z$ for the "resting" observer,
as well as an extra pair $\xi^{\prime}$ and $\chi^{\prime}$ for the
superluminal one. We will also assume that the "transversal" components of the
four-position are not affected by the superluminal transformation. It follows
that for the motion along the $x$ axis with a superluminal velocity the
spacetime interval undergoes the following transformation:
$c^{2}\text{d}t^{2}-\text{d}\bm{r}\cdot\text{d}\bm{r}=-c^{2}\text{d}t^{\prime
2}+\text{d}x^{\prime 2}-\text{d}\xi^{\prime 2}-\text{d}\chi^{\prime 2}$, so
the spacetime metric is affected, but the null intervals are still preserved.
This result indicates that the laws of physics in the inertial superluminal
frame of reference are different from those within the orthodox family of
subluminal frames. This agrees with the argument put forth by Machildon et al.
[13]. Furthermore, Sutherland and Shepanski argue [18] that this is due to the
fact the the "spatial" component of the four-position in the superluminal
frame, defined as $(x^{\prime},\xi^{\prime},\chi^{\prime})$, characterizes a
non-Euclidean space. In this work we follow a different interpretation first
mentioned by Demers [19] and developed in [1]: the signs of the metric
components indicate that the extra pair of dimensions $\xi^{\prime}$ and
$\chi^{\prime}$ are in fact temporal, because they enter the metric with the
same sign as the temporal coordinate $t^{\prime}$. According to this
interpretation a superluminal frame of reference is equiped with a single
spatial dimension that we will denote with $r^{\prime}$, as well as three
temporal dimensions $\bm{t^{\prime}}$. In this paper, we will investigate
physical consequences of such a hypothesis.

Let us consider a subluminal frame $(ct,\bm{r})$ and an unrotated,
superluminal frame $(c\bm{t^{\prime}},r^{\prime})$ moving with a superluminal
velocity $\bm{V}$. Here, by "unrotated" we mean that the resting frame should
move relative to the superluminal frame with velocity $-\bm{V}$, so that the
inverse transformation leads to a sign flip in relative velocity. First off,
we note that the superluminal transformation apparently leaves the transversal
spatial components of four-position unaffected:
$c\bm{t^{\prime}}-c\frac{\bm{t^{\prime}}\cdot\bm{V}}{V^{2}}\bm{V}=\bm{r}-\frac{\bm{r}\cdot\bm{V}}{V^{2}}\bm{V}$.
The longitudinal component along the direction $\bm{V}$ undergoes the
transformation given by (II):
$\frac{\bm{t^{\prime}}\cdot\bm{V}}{V^{2}}\bm{V}=-\frac{\bm{V}}{V}\frac{t-\frac{\bm{r}\cdot\bm{V}}{c^{2}}}{\sqrt{\frac{V^{2}}{c^{2}}-1}}$
and the spatial coordinate transforms according to
$r^{\prime}=-\frac{\frac{\bm{r}\cdot\bm{V}}{V}-Vt}{\sqrt{\frac{V^{2}}{c^{2}}-1}}$.
Adding up the first two equations and supplementing the result with the third
one yields [1, 14]:

$\displaystyle r^{\prime}$ $\displaystyle=$
$\displaystyle\frac{Vt-\frac{\bm{r}\cdot\bm{V}}{V}}{\sqrt{\frac{V^{2}}{c^{2}}-1}},$
$\displaystyle c\bm{t}^{\prime}$ $\displaystyle=$
$\displaystyle\bm{r}-\frac{\bm{r}\cdot\bm{V}}{V^{2}}\bm{V}+\frac{\frac{\bm{r}\cdot\bm{V}}{Vc}-\frac{ct}{V}}{\sqrt{\frac{V^{2}}{c^{2}}-1}}\bm{V},$
(4)

which is the coordinate transformation between a subluminal and a superluminal
inertial observer in $1+3$ spacetime. The purpose of this paper is to
investigate physical consequences of (II).

As mentioned before, the $V\to 0$ limit of this transformation does not exist,
but the infinite velocity limit does. For $V\to\infty$, the formulae (II)
reduce to:

$\displaystyle r^{\prime}$ $\displaystyle=$ $\displaystyle ct,$ $\displaystyle
c\bm{t}^{\prime}$ $\displaystyle=$ $\displaystyle\bm{r},$ (5)

regardless of the direction of velocity $\bm{V}$. The spacetime interval
undergoes the following transformation under (II):

$\displaystyle
c^{2}\text{d}t^{2}-\text{d}\bm{r}\cdot\text{d}\bm{r}=-c^{2}\text{d}\bm{t^{\prime}}\cdot\text{d}\bm{t^{\prime}}+\text{d}r^{\prime
2},$ (6)

which does not depend on $\bm{V}$. This signifies the fact that laws of
physics are the same across the whole family of superluminal inertial
observers, although these laws differ from the ones common among subluminal
inertial observers due to a different metric. In the following we will adopt a
convention, in which the spacetime metric in subluminal frames is
$\eta_{\mu\nu}\equiv\text{diag}(1,-1,-1,-1)$, and in superluminal frames it is
defined as $\eta_{\mu\nu}\equiv\text{diag}(1,1,1,-1)$, so that the four-
positions $x^{\mu}\equiv(ct,\bm{r})$ and
$x^{\prime\mu}\equiv(c\bm{t^{\prime}},r^{\prime})$ and all other four-vectors
are related by:

$\displaystyle\eta_{\mu\nu}x^{\mu}x^{\nu}=-\eta^{\prime}_{\mu\nu}x^{\prime\mu}x^{\prime\nu}.$
(7)

The superluminal transformation (II) can be also rewritten in the following
form:

$\displaystyle r^{\prime}$
$\displaystyle=\frac{ct-\frac{c\bm{r}\cdot\bm{V}}{V^{2}}}{\sqrt{1-\frac{c^{2}}{V^{2}}}},$
$\displaystyle c\bm{t}^{\prime}$
$\displaystyle=\bm{r}-\frac{\bm{r}\cdot\bm{V}}{V^{2}}\bm{V}+\frac{\frac{\bm{r}\cdot\bm{V}}{V^{2}}-\frac{c^{2}t}{V^{2}}}{\sqrt{1-\frac{c^{2}}{V^{2}}}}\bm{V},$
(8)

in which it becomes clear that it can be decomposed into the subluminal
Lorentz transformation (II) for velocity
$\bm{V}\to\frac{c^{2}}{V}\frac{\bm{V}}{V}$ followed by the infinite
transformation (II). The first operation takes an object moving with a
superluminal speed $\bm{V}$ to a frame, in which the object moves infinitely
fast, the latter goes to the rest frame of the object, in which $r^{\prime}$
is constant. Any superluminal transformation (II) can be decomposed into a
pair of such operations.

In subluminal relativity the worldline of a moving observer coincides with his
temporal axis $t^{\prime}$ and it uniquely defines the $1+3$ decomposition of
spacetime into time and space, as the three-dimensional space
$\bm{r^{\prime}}$ is orthogonal to the time axis $t^{\prime}$. In order to
uniquely define the orientation of the spacial axes $\bm{r^{\prime}}$
additional information has to be provided. To supplement the missing
information one usually assumes that the two moving frames are mutually
unrotated, which is equivalent of saying that the inverse transformation
involves a simple sign flip in relative velocity.

In the case of superluminal observers the situation is different, but some
analogies remain. A trajectory of the superluminal observer does not
necessarily coincide with neither of the three temporal axes
$\bm{t^{\prime}}$. Also the one-dimensional space $r^{\prime}$ of the
superluminal observer is not uniquely defined by his world line and additional
information has to be provided to define it. Our derivation may be interpreted
as a way to choose the $r^{\prime}$ axis. The missing information is again
provided by assuming that the two considered reference frames are mutually
unrotated. These observations are a part of a bigger picture that will be
better understood once we deal with dynamics. In particular we will show that
point-like trajectories need to be rejected anyway, as incompatible with
relativistic requirements involving superluminal observers.

Our interpretation that the superluminal observer characterizes spacetime
using three temporal dimensions $\bm{t^{\prime}}$ poses several
interpretational challenges. The first question: how to even define a velocity
in a frame that has more than one temporal dimension? In order to answer this
question, let us determine the inverse superluminal transformation to (II). It
can be obtained by an algebraic reversal, which turns out to be equivalent to
substituting $\bm{V}\to-\bm{V}$ and interchanging primes according to
substitutions $r^{\prime}\leftrightarrow ct$ and
$c\bm{t^{\prime}}\leftrightarrow\bm{r}$:

$\displaystyle ct$ $\displaystyle=$
$\displaystyle\frac{\frac{Vr^{\prime}}{c}+\frac{\bm{V}\cdot
c\bm{t^{\prime}}}{V}}{\sqrt{\frac{V^{2}}{c^{2}}-1}},$ $\displaystyle\bm{r}$
$\displaystyle=$ $\displaystyle c\bm{t^{\prime}}-\frac{\bm{V}\cdot
c\bm{t^{\prime}}}{V^{2}}\bm{V}+\frac{\frac{\bm{V}\cdot\bm{t^{\prime}}}{V}+\frac{r^{\prime}}{V}}{\sqrt{\frac{V^{2}}{c^{2}}-1}}\bm{V}.$
(9)

The question that remains is whether $-\bm{V}$ is an actual velocity of the
frame $(ct,\bm{r})$ observed in the frame $(c\bm{t^{\prime}},r^{\prime})$, as
required by the fact that the two frames are mutually unrotated, and how such
a velocity is even defined? Let us take the point of view of the superluminal
primed observer and consider the origin of the "resting", unprimed frame
$\bm{r}=0$. We should expect that according to the right definition of
velocity, that origin should move with velocity $-\bm{V}$. By substituting
$\bm{r}=0$ into (II) we find:

$\displaystyle r^{\prime}$ $\displaystyle=$
$\displaystyle\frac{Vt}{\sqrt{\frac{V^{2}}{c^{2}}-1}},$ $\displaystyle
c\bm{t}^{\prime}$ $\displaystyle=$
$\displaystyle-\frac{\frac{ct}{V}\bm{V}}{\sqrt{\frac{V^{2}}{c^{2}}-1}}.$ (10)

Elimination of $t$ from the equations gives
$r^{\prime}\frac{\bm{V}}{V^{2}}=-\bm{t^{\prime}}$. This result inspires the
definition of a velocity $\bm{v^{\prime}}$ of an object moving in a
superluminal frame such that the following equation is satisfied:

$\displaystyle r^{\prime}\frac{\bm{v^{\prime}}}{v^{\prime 2}}=\bm{t^{\prime}}$
(11)

and such a velocity also has the property that:

$\displaystyle r^{\prime}=\bm{v^{\prime}}\cdot\bm{t^{\prime}}.$ (12)

We can now determine $v^{\prime 2}$ by taking a scalar product of (11) with
itself and plug it back into (11). A rearranged differential form of (11) is
then:

$\displaystyle\bm{v^{\prime}}=\frac{\text{d}r^{\prime}}{\text{d}t^{\prime}}\frac{\text{d}\bm{t^{\prime}}}{\text{d}t^{\prime}},$
(13)

where
$\text{d}t^{\prime}\equiv\sqrt{\text{d}\bm{t^{\prime}}\cdot\text{d}\bm{t^{\prime}}}$.
Equation (13) will provide our definition of an instantaneous velocity in a
superluminal frame.

Let us now use the definition (13) to determine the superluminal velocity
composition formula. Consider an object moving in the subluminal frame with
any velocity $\bm{v}$. Let us transform this velocity to the superluminal
frame. By substituting $\bm{r}=\bm{v}t$ into (II) we obtain:

$\displaystyle r^{\prime}$ $\displaystyle=$
$\displaystyle\frac{Vt-\frac{\bm{v}\cdot\bm{V}}{V}t}{\sqrt{\frac{V^{2}}{c^{2}}-1}},$
$\displaystyle c\bm{t}^{\prime}$ $\displaystyle=$
$\displaystyle\bm{v}t-\frac{\bm{v}\cdot\bm{V}}{V^{2}}\bm{V}t+\frac{\frac{\bm{v}\cdot\bm{V}}{Vc}-\frac{c}{V}}{\sqrt{\frac{V^{2}}{c^{2}}-1}}\bm{V}t$
(14)

and eliminating $t$ from the equations we get:

$\displaystyle
r^{\prime}\frac{\sqrt{\frac{V^{2}}{c^{2}}-1}\left(\bm{v}-\frac{\bm{v}\cdot\bm{V}}{V^{2}}\bm{V}\right)-c\left(1-\frac{\bm{v}\cdot\bm{V}}{c^{2}}\right)\frac{\bm{V}}{V}}{V-\frac{\bm{v}\cdot\bm{V}}{V}}$
$\displaystyle=c\bm{t}^{\prime}.$ (15)

Now we can compare our result (15) with the equation (11) to find the
superluminal velocity addition formula:

$\displaystyle\frac{c^{2}}{v^{\prime
2}}\bm{v^{\prime}}=\frac{\sqrt{1-\frac{c^{2}}{V^{2}}}\left(\bm{v}-\frac{\bm{v}\cdot\bm{V}}{V^{2}}\bm{V}\right)-\left(\frac{c^{2}}{V^{2}}\bm{V}-\frac{\bm{v}\cdot\bm{V}}{V^{2}}\bm{V}\right)}{1-\frac{\bm{v}\cdot\bm{V}}{V^{2}}},$
(16)

where $\bm{v^{\prime}}$ is the velocity of the considered object in the
superluminal frame of reference. Notice that the result (16) can be obtained
directly from (2) by substitution $\bm{V}\to\frac{c^{2}}{V}\frac{\bm{V}}{V}$
and
$\bm{v^{\prime}}\to\frac{c^{2}}{v^{\prime}}\frac{\bm{v^{\prime}}}{v^{\prime}}$.
Finally, we take a square of (16) to compute $v^{\prime 2}$, plug it back and
after some simplifications we find the velocity composition formula for a
superluminal velocity $\bm{V}$, which is a superluminal version of (2):

$\displaystyle\bm{v^{\prime}}$
$\displaystyle=\frac{\sqrt{1-\frac{c^{2}}{V^{2}}}\left(\bm{v}-\frac{\bm{v}\cdot\bm{V}}{V^{2}}\bm{V}\right)-\left(\frac{c^{2}}{V^{2}}\bm{V}-\frac{\bm{v}\cdot\bm{V}}{V^{2}}\bm{V}\right)}{1-\frac{\bm{v}\cdot\bm{V}}{V^{2}}}\left(1-\frac{\left(1-\frac{c^{2}}{V^{2}}\right)\left(1-\frac{v^{2}}{c^{2}}\right)}{\left(1-\frac{\bm{v}\cdot\bm{V}}{V^{2}}\right)^{2}}\right)^{-1}.$
(17)

By squaring the above result we also find the relation between speeds:

$\displaystyle\left(1-\frac{c^{2}}{v^{\prime
2}}\right)=\frac{\left(1-\frac{c^{2}}{V^{2}}\right)\left(1-\frac{v^{2}}{c^{2}}\right)}{\left(1-\frac{\bm{v}\cdot\bm{V}}{V^{2}}\right)^{2}}.$
(18)

The result (18) shows that if $v=c$ then $v^{\prime}=c$. Therefore our
superluminal transformation (II) preserves the speed of light. Moreover, if
$v>c$ then $v^{\prime}<c$ and vice versa. Therefore subluminal objects move
with superluminal speeds according to any superluminal observer and
superluminal objects according to any superluminal observer move with
subluminal speeds. We must remember, however, that since both subluminal and
superluminal families of observers are physically distinguishable, we can
apply these adjectives to observers in an absolute sense. The constancy of the
speed of light across all inertial frames is a property characteristic to the
definition of velocity derived in (13). Notice that alternative, incorrect
definitions, such as one in [20] may lead to different conclusions.

The inverse transformation to (16) can be obtained either by using (II) to
evaluate $\bm{v}\equiv\frac{\text{d}\bm{r}}{\text{d}t}$ or by using (16) with
the substitution: $\bm{V}\to-\bm{V}$ and $\frac{c^{2}}{v^{\prime
2}}\bm{v^{\prime}}\leftrightarrow\bm{v}$. The result is the same:

$\displaystyle\bm{v}$
$\displaystyle=\frac{\sqrt{1-\frac{c^{2}}{V^{2}}}\left(\frac{c^{2}\bm{v^{\prime}}}{v^{\prime
2}}-\frac{c^{2}}{v^{\prime
2}}\frac{\bm{v^{\prime}}\cdot\bm{V}}{V^{2}}\bm{V}\right)+\left(\frac{c^{2}\bm{V}}{V^{2}}+\frac{c^{2}}{v^{\prime
2}}\frac{\bm{v^{\prime}}\cdot\bm{V}}{V^{2}}\bm{V}\right)}{1+\frac{c^{2}}{v^{\prime
2}}\frac{\bm{v^{\prime}}\cdot\bm{V}}{V^{2}}}.$ (19)

Now, that the kinematics is well defined in superluminal frames of reference,
and the corresponding transformation laws are established, let us discuss the
expressions for energy and momentum of superluminal particles.

##  III Energy and momentum

Let us begin by looking for the characterization of energy and momentum of a
point-like particle moving in a subluminal frame with a superluminal velocity
$\bm{v}$. As in conventional special relativity, we require that energy and
momentum of a superluminal particle form a four-vector structure. This will
guarantee that if energy and momentum are conserved for a certain process in
one inertial frame, they will also be conserved for this process in all other
reference frames. Such a process can include interactions between both sub-
and superluminal objects, as well as light.

Similarly to conventional special relativity we will assume that the four-
momentum $p^{\mu}$ of a superluminal particle is proportional to four-velocity
$u^{\mu}$, but this time we will take the latter to be a space-like four-
vector with $u^{\mu}u_{\mu}=-c^{2}$.

Consider a worldline of a superluminal object characterized by the four-
velocity $u^{\mu}\equiv(u^{0},\bm{u})$. We will look for this four-vector in a
form analogous to the subluminal expression for which
$\bm{u}=u^{0}\frac{\bm{v}}{c}$, so that $u^{\mu}=u^{0}(1,\frac{\bm{v}}{c})$.
Notice that the direction of velocity
$\bm{v}=\frac{\text{d}\bm{r}}{\text{d}t}$ and the spacial component $\bm{u}$
are either parallel or antiparallel, depending on the sign of $u^{0}$. Since
the four-vector $u^{\mu}$ is spacelike, even the subluminal Lorentz
transformation (II) can change the sign of $u^{0}$ and reverse the mutual
orientation between $\bm{v}$ and $\bm{u}$. Since we require energy $E$ of the
particle to be proportional to $u^{0}$ and momentum $\bm{p}$ proportional to
$\bm{u}$, therefore we should be aware that the sign of energy will be frame
dependent. Therefore a superluminal particle can either have positive energy
and momentum parallel to velocity, or negative energy and momentum
antiparallel to velocity in another frame of reference. This ambiguity results
from the fact that specifying the mass and velocity of a superluminal object
is insufficient to determine its energy and momentum. We also have to specify,
whether the object in question is a particle, or an antiparticle moving
backwards in time. Furthermore we should expect that in quantum field theory a
tachyonic particle can be transformed into its antiparticle via a mere
subluminal Lorentz transformation (II). This is essentially the same argument,
as the one used by Feynman in his famous lecture on the relativistic reason
for antiparticles [21], in which he was focusing on superluminal virtual
particles.

Let us compute the length of the superluminal four-velocity:

$\displaystyle c=\sqrt{-u^{\mu}u_{\mu}}$
$\displaystyle=\sqrt{\bm{u}\cdot\bm{u}-\left(u^{0}\right)^{2}}=\absolutevalue{u^{0}}\sqrt{\frac{v^{2}}{c^{2}}-1},$
(20)

which gives:

$\displaystyle\absolutevalue{u^{0}}=\frac{c}{\sqrt{\frac{v^{2}}{c^{2}}-1}},$
(21)

that can be used together with $u^{0}=\text{sgn}(u^{0})|u^{0}|$ to define the
the four-momentum in the following way:

$\displaystyle p^{\mu}\equiv
mu^{\mu}=\text{sgn}(u^{0})\left(\frac{mc\,}{\sqrt{\frac{v^{2}}{c^{2}}-1}},\frac{m\bm{v}\,}{\sqrt{\frac{v^{2}}{c^{2}}-1}}\right),$
(22)

where $m$ is the superluminal mass of the particle in analogy to subluminal
particles. An identical procedure in another subluminal reference frame gives:

$\displaystyle p^{\prime\mu}=\text{sgn}(u^{\prime
0})\left(\frac{mc\,}{\sqrt{\frac{v^{\prime
2}}{c^{2}}-1}},\frac{m\bm{v^{\prime}}\,}{\sqrt{\frac{v^{\prime
2}}{c^{2}}-1}}\right)$ (23)

and four-vectors (22) and (23) are related via a subluminal Lorentz
transformation (II). This relation for the temporal component yields:

$\displaystyle u^{\prime 0}$
$\displaystyle=\frac{u^{0}-\frac{\bm{u}\cdot\bm{V}}{c}}{\sqrt{1-\frac{V^{2}}{c^{2}}}},$
(24)

which leads to the transformation properties of $\text{sgn}(u^{0})$:

$\displaystyle\text{sgn}(u^{\prime 0})$
$\displaystyle=\text{sgn}\left(\frac{u^{0}-\frac{\bm{u}\cdot\bm{V}}{c}}{\sqrt{1-\frac{V^{2}}{c^{2}}}}\right)=\text{sgn}\left(u^{0}-\frac{\bm{u}\cdot\bm{V}}{c}\right)$
$\displaystyle=\text{sgn}(u^{0})\,\text{sgn}\left(1-\frac{\bm{v}\cdot\bm{V}}{c^{2}}\right).$
(25)

Finally, the relativistically covariant expressions for energy and momentum of
a superluminal particle in subluminal reference frames are given by:

$\displaystyle E$ $\displaystyle\equiv\frac{\sigma
mc^{2}}{\sqrt{\frac{v^{2}}{c^{2}}-1}},$ $\displaystyle\bm{p}$
$\displaystyle\equiv\frac{\sigma m\bm{v}}{\sqrt{\frac{v^{2}}{c^{2}}-1}},$ (26)

where $\sigma=\pm 1$ carrying the information whether the object is a particle
or an antiparticle transforms according to the equation:

$\displaystyle\sigma^{\prime}=\sigma\,\text{sgn}\left(1-\frac{\bm{v}\cdot\bm{V}}{c^{2}}\right),$
(27)

which was first justified heuristically in [22]. The important observation is
that superluminal objects cannot be slowed down under the speed of light, as
this would require infinite supplies of energy. Also notice that an object
moving with an infinite velocity carries momentum, but no energy, as
$\lim_{v\to\infty}\absolutevalue{\bm{p}}=mc$ and $\lim_{v\to\infty}E=0$. This
leads to an interesting novelty. In conventional relativity a particle cannot
emit another particle without changing its own mass [14]. For example a free
electron cannot emit a photon without violating conservation of energy. The
situation is different if we include superluminal particles in our
considerations. For example, a subluminal particle could emit an infinitely
fast moving particle and simply reverse its own velocity due to recoil. In
such a process, the energy is conserved, because the source does not change
its energy and the emitted particle carries no energy. Momentum conservation
can also be satisfied by an appropriate choice of the initial velocity or mass
of the object. In principle, such an emission could be repeated a number of
times and the object could emit a number of infinitely fast moving particles
[15]. Notice that such an emission is only possible after reaching certain
velocity threshold. If velocity of the "source" is below the threshold, only
the absorption of an infinitely moving particle is permitted unless one
chooses to interpret such a process as an emission of an anti-particle
backwards in time [15].

Energy and momentum can also be transformed to superluminal frames of
reference. Similarly to four-position, the four-momentum in these frames
contains the energy three-vector and a single momentum. Applying the
superluminal transformation (II) to the four-vector (22) we obtain:

$\displaystyle p^{\prime}$ $\displaystyle=$
$\displaystyle\frac{\frac{VE}{c^{2}}-\frac{\bm{p}\cdot\bm{V}}{V}}{\sqrt{\frac{V^{2}}{c^{2}}-1}},$
$\displaystyle\frac{\bm{E}^{\prime}}{c}$ $\displaystyle=$
$\displaystyle\bm{p}-\frac{\bm{p}\cdot\bm{V}}{V^{2}}\bm{V}+\frac{\frac{\bm{p}\cdot\bm{V}}{Vc}-\frac{E}{Vc}}{\sqrt{\frac{V^{2}}{c^{2}}-1}}\bm{V}.$
(28)

Transformation to the rest frame of the superluminal particle, as previously,
can also be decomposed into two steps: first transforming the four-momentum to
a frame, in which the object moves with an infinite velocity in some direction
$\bm{s}$ and then applying the analogue of (II) that interchanges energy with
momentum. This results in the four-momentum taking the form:

$\displaystyle p^{\prime\mu}=(mc\bm{s},0),$ (29)

which is not unique and corresponds to a bundle of infinitely many worldlines
for each choice of $\bm{s}$. This shows that there are infinitely many
observers characterized by non-parallel world lines, each given by a different
$\bm{s}$, that are at relative rest. To confirm this it suffices, to inspect
the definition of velocity in superluminal frames (13) that shows that if
$\text{d}r^{\prime}=0$ then the velocity vanishes. This interesting conundrum
will be clarified once we turn our attention to the dynamical laws including
superluminal extension of relativity.

##  IV Dynamics

In classical mechanics, a dynamical evolution of a system is obtained by
minimizing the action:

$S\equiv\int L\,\text{d}t.$ (30)

Frame-independence of the result is guaranteed by demanding $S$ to be the same
in all (subluminal) frames and for a single point-like particle the result is
a one-dimensional trajectory $\bm{r}(t)$. Now we wish to generalize this
principle of least action to include superluminal frames. As we now have to
deal with three time variables, $\bm{t}^{\prime}$, we hypothesize that the
corresponding variational principle will involve an action of the form:

$S^{\prime}\equiv\int L^{\prime}\,\text{d}^{3}t^{\prime},$ (31)

where $L^{\prime}$ is an appropriately defined superluminal counterpart of the
Lagrangian. Clearly, the resulting "trajectory" of a single particle
$r^{\prime}(\bm{t^{\prime}})$ defines a three-dimensional surface, which is
apparently at odds with subluminal intuitions. It is therefore interesting to
explore this notion in more detail.

We are seeking to write down $L^{\prime}$ such that $S^{\prime}$ is
relativistically invariant. With one time variable, it would be accomplished
with $L$ proportional to the length of the one-dimensional trajectory
$\bm{r}(t)$, i.e. we would take

$L\propto\sqrt{1-\frac{1}{c^{2}}\left(\frac{\mathrm{d}\bm{r}}{\mathrm{d}t}\right)^{2}}.$
(32)

Similarly, with three time variables, the natural choice would be to consider
the three-dimensional volume of the manifold parametrized as $r(\bm{t})$, i.e.
to consider

$L^{\prime}\propto\sqrt{\frac{1}{c^{2}}\left(\bm{\nabla_{\\!t^{\prime}}}r^{\prime}\right)^{2}-1},$
(33)

where $\bm{\nabla_{\\!t^{\prime}}}$ denotes a gradient with respect to the
time variables. The choice of signs in the square root reflects the fact that
we are interested in motions which are superluminal in the chosen reference
frame, i.e. we characterize the example of subluminal particles in a
superluminal frame. With such a choice of $L^{\prime}$, the action (31) is
evaluated over $\mathbb{R}^{3}$ of the time variables. This is strikingly
different from a one-dimensional case, in which instants of time can be
ordered, so we can evaluate the action (30) over a finite interval.

Demanding that the action (31) is stationary with respect to the variations of
the manifold $r^{\prime}(\bm{t^{\prime}})\to
r^{\prime}(\bm{t}^{\prime})+\delta r^{\prime}(\bm{t^{\prime}})$ leads to:

$\displaystyle 0=\delta S^{\prime}$
$\displaystyle\propto\int\frac{\bm{\nabla_{\\!t^{\prime}}}r^{\prime}\cdot\bm{\nabla_{\\!t^{\prime}}}\delta
r^{\prime}}{\sqrt{\frac{1}{c^{2}}\left(\bm{\nabla_{\\!t^{\prime}}}r^{\prime}\right)^{2}-1}}\,\mathrm{d}t^{\prime
3}=\int\frac{\left(\left(\frac{1}{c^{2}}(\bm{\nabla_{\\!t^{\prime}}}r^{\prime})^{2}-1\right)\delta_{ij}-\frac{1}{c^{2}}\frac{\partial
r^{\prime}}{\partial t^{\prime}_{i}}\frac{\partial r^{\prime}}{\partial
t^{\prime}_{j}}\right)\frac{\partial^{2}r^{\prime}}{\partial
t^{\prime}_{i}\partial
t^{\prime}_{j}}}{\left(\frac{1}{c^{2}}\left(\bm{\nabla_{\\!t^{\prime}}}r^{\prime}\right)^{2}-1\right)^{3/2}}\,\delta
r^{\prime}\,\mathrm{d}t^{\prime 3},$ (34)

where in the last equality we assumed that the boundary term vanishes at
infinity and we can evaluate the integral by parts; summation over repeated
indices is understood. Because $\delta r^{\prime}$ is an arbitrary function of
$\bm{t^{\prime}}$, we conclude that the stationary action corresponds to the
equation of motion

$\left(\left(\frac{1}{c^{2}}(\bm{\nabla_{\\!t^{\prime}}}r^{\prime})^{2}-1\right)\delta_{ij}-\frac{1}{c^{2}}\frac{\partial
r^{\prime}}{\partial t^{\prime}_{i}}\frac{\partial r^{\prime}}{\partial
t^{\prime}_{j}}\right)\frac{\partial^{2}r^{\prime}}{\partial
t^{\prime}_{i}\partial t^{\prime}_{j}}=0.$ (35)

Because (35) contains second-order time derivatives, it is solved by functions
of the form:

$r^{\prime}(\bm{t^{\prime}})=r^{\prime}_{0}+\bm{v^{\prime}}\cdot\bm{t}^{\prime},$
(36)

where $r^{\prime}_{0}$ and $\bm{v^{\prime}}$ are constant. Those solutions,
that are compatible with the equation of constant motion (12), can be viewed
as "sheets" of particles traveling with velocity $\bm{v^{\prime}}$
perpendicular to the sheets. This behavior corresponds to the motion of a
collection of point-like particles described in a subluminal frame.
Unfortunately, this is where the similarities of the two points of view end.
In a spherically symmetric case, eq. (35) becomes:

$\frac{\mathrm{d}^{2}r^{\prime}}{\mathrm{d}{t^{\prime}}^{2}}+\frac{2}{t^{\prime}}\frac{\mathrm{d}r^{\prime}}{\mathrm{d}{t^{\prime}}}-\frac{2}{c^{2}t^{\prime}}\left(\frac{\mathrm{d}r^{\prime}}{\mathrm{d}{t^{\prime}}}\right)^{3}=0\,,$
(37)

where $t^{\prime}\equiv|\bm{t^{\prime}}|$. Solutions of eq. (37) include a
constant, a spherical shell moving with the speed of light and in other cases
they are expressed in terms of elliptic functions. In particular, the latter
solution does not correspond to any known representation of a free particle or
a collection thereof.

The problems that we encountered after swapping the temporal and spatial
coordinates by means of a superluminal transformation follow from the
dimensionality of space. Were there only one spatial dimension, the partial
differential equations involving $\bm{\nabla_{\\!t^{\prime}}}$ would become
ordinary differential equations and linear trajectories of point-lime
particles would remain linear trajectories.

Furthermore, the choice of (33) is not unique. However, any action that
respects the shift and rotational symmetry of $\bm{t}^{\prime}$ will depend
only on $(\bm{\nabla_{\\!t^{\prime}}}r)^{2}$, which will lead to the second
order differential equation like (35) for a free particle. In this sense, the
derivation presented above is quite general, as (36) would still be the
solution of the resulting equation of motion.

What is even worse, the formalism described above gives us no chances to
arrive at a relation between the actions (30) and (31) as their solutions are
clearly non-equivalent. We are therefore forced to conclude that our approach
of treating $r^{\prime}(\bm{t}^{\prime})$ as the parametrization of the
configuration space of the system fails badly on many fronts.

Let us take a step back and try to formulate our requirements for a
satisfactory description of a physical system which is relativistically
invariant and, in particular, does not break down for a superluminal
transformation. To this end, let us symbolically denote the state of the
system fully characterizing it in a configuration space in a subluminal frame
by $\psi$ and similarly in the superluminal frame by $\psi^{\prime}$. We will
also assume that the transition to the superluminal frame transforms
$\psi\to\psi^{\prime}$ so that the action $S$ given by (30) can be written
down in both subluminal and superluminal frames as:

$\displaystyle S\equiv\int L[\psi]\text{d}t=\int
L[\psi^{\prime}]\frac{\text{d}r^{\prime}}{c},$ (38)

where in the last equality we applied the transformation (II). We also have
similar relations for the superluminal action $S^{\prime}$ given by (31):

$\displaystyle S^{\prime}\equiv\int
L^{\prime}[\psi^{\prime}]\text{d}^{3}t^{\prime}=\int
L^{\prime}[\psi]\frac{\text{d}^{3}r}{c^{3}}.$ (39)

In order to guarantee that the principle of least action leads to equivalent
solutions in all frames we will enforce that the actions (38) and (39) are
equal: $S=S^{\prime}$. This leads to the following equalities:

$\displaystyle\int L[\psi]\text{d}t=\int
L^{\prime}[\psi]\frac{\text{d}^{3}r}{c^{3}},$ (40)

and

$\displaystyle\int L[\psi^{\prime}]\frac{\text{d}r^{\prime}}{c}=\int
L^{\prime}[\psi^{\prime}]\text{d}^{3}t^{\prime}.$ (41)

It follows that the characterization of the state in terms of a "trajectory"
is no longer possible, because the Lagrangian $L^{\prime}$ must now be
integrated over $\text{d}^{3}r$ and the Lagrangian $L$ must be integrated over
$\text{d}r^{\prime}$, therefore $\psi\neq\bm{r}(t)$ and $\psi^{\prime}\neq
r^{\prime}(\bm{t^{\prime}})$. This circumstance brings us to the conclusion
that it is necessary to introduce a new type of configuration space and a new
notion of a state of the system, $\psi(t,\bm{r})$, that depends on all
spacetime coordinates, over which the action will be minimized. For such
states, the function encoding the state of the system has to be integrated
over both time and space, hence it should be interpreted as a Lagrangian
density of the system. Of course, such a physical picture brings us outside of
the domain of classical mechanics and forces us to consider a familiar
framework involving fields as the only possible descriptions of physical
systems.

Let us now observe that if the Lagrangian $L$ on the left-hand side of (40)
only involves integration over $\text{d}t$, while $L^{\prime}$ on the right-
hand side over $\text{d}^{3}r$, we would face the equality of a pair of
functions of different variables: $\bm{r}$ on the left-hand side and $t$ on
the right-hand side. In order to avoid this without introducing preference for
any of the spacetime point, we are forced to impose that the Lagrangians must
be expressed in terms of Lagrangian densities: $L[\psi]\equiv\int{\cal
L}[\psi]\,\text{d}^{3}r$ and $L^{\prime}[\psi^{\prime}]\equiv\int{\cal
L}[\psi^{\prime}]\,\text{d}r^{\prime}$, so that the additional variables are
also integrated over. This also solves an analogous difficulty with equation
(41). All in all, we find that extending the special relativistic dynamics to
superluminal observers terminates any chance to consider mechanical systems
and forces one to consider field-theoretic configuration space.

At this stage it remains unspecified, whether the resulting "field theory"
should be classical or quantum and in order to determine the answer, one has
to refer to arguments put forth in [1] in favor of the quantum theory. What
remains out of question is that the superluminal extension of special
relativity in $1+3$ spacetime is possible in a fully consistent way. However,
for this task to be covariant it is necessary to abandon the classical picture
of mechanical systems characterized by unique, one-dimensional trajectories
and upgrade physics to the field-theoretic framework.

##  V Example: Maxwell’s theory

In Maxwell’s theory the state of the system is characterized by
electromagnetic potentials and their derivatives:
$\psi\to\\{A^{\mu},\partial_{\nu}A^{\mu}\\}$ and we have the following
Lagrangian density in Lorenz gauge $\partial_{\mu}A^{\mu}=0$:

$\displaystyle{\cal
L}(A^{\mu})\equiv-\frac{1}{4\mu_{0}}F^{\mu\nu}F_{\mu\nu}-A^{\mu}j_{\mu},$ (42)

where $F_{\mu\nu}\equiv\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$ and
$\partial_{\mu}\equiv(\frac{1}{c}\partial_{t},\bm{\nabla})$. The resulting
Euler-Lagrange equations are:

$\displaystyle\partial^{\mu}\partial_{\mu}A^{\nu}=\mu_{0}j^{\nu}.$ (43)

Let us now formulate the theory in superluminal reference frames. This problem
has already been addressed by Dawe and Hines [23, 24], however they adopted
the interpretation in which superluminal observers are characterized by a
single temporal dimension and non-Euclidean space [25]. Our approach, assuming
the three-dimensional time variable, provides a much cleaner mathematical
framework and the task of deriving field equations becomes practically
trivial. Since an arbitrary superluminal transformation (II) can be decomposed
into a subluminal boost (which leaves Maxwell’s equations intact) and the
infinite velocity flip (II), it is sufficient to only consider the latter.

The infinitely fast superluminal transformation (II) imposes that
$\partial^{\prime}_{\mu}=(\frac{1}{c}\bm{\nabla_{\\!t^{\prime}}},\partial_{r^{\prime}})$,
as well as the following set of conditions for the four-potentials $A^{\mu}$
and four-current $j^{\mu}$:

$\displaystyle A^{\prime\mu}$
$\displaystyle\equiv(\frac{1}{c}\bm{\varphi^{\prime}},A^{\prime})=(\bm{A},\frac{1}{c}\varphi),$
$\displaystyle j^{\prime\mu}$
$\displaystyle\equiv(\bm{\varrho^{\prime}}c,j^{\prime})=(\bm{j},\varrho c).$
(44)

The transformation procedure characterized in the previous Section requires
that the Lagrangian density (42) and dynamical equations (43) in the
superluminal frame have the analogous form:

$\displaystyle{\cal
L}(A^{\prime\mu})=-\frac{1}{4\mu_{0}}F^{\prime\mu\nu}F^{\prime}_{\mu\nu}+A^{\prime\mu}j^{\prime}_{\mu},$
(45)

where
$F^{\prime}_{\mu\nu}\equiv\partial^{\prime}_{\mu}A^{\prime}_{\nu}-\partial^{\prime}_{\nu}A^{\prime}_{\mu}$
and the sign change in the source term is due to the equation (7), so that
$\eta_{\mu\nu}A^{\mu}j^{\nu}=-\eta^{\prime}_{\mu\nu}A^{\prime\mu}j^{\prime\nu}$.
The resulting dynamical equations are

$\displaystyle-\partial^{\prime\mu}\partial^{\prime}_{\mu}A^{\prime\nu}=\mu_{0}j^{\prime\nu},$
(46)

which can be also obtained directly by applying (II) to (43). Our construction
of electromagnetic potentials in superluminal frames is manifestly covariant.
We can also introduce electric and magnetic fields in subluminal frames in a
standard way:

$\displaystyle\bm{E}$
$\displaystyle\equiv-\bm{\nabla}\varphi-\partial_{t}\bm{A},$
$\displaystyle\bm{B}$ $\displaystyle\equiv\bm{\nabla}\times\bm{A}$ (47)

and define these fields in superluminal frames using the following analogies:

$\displaystyle\bm{E^{\prime}}$
$\displaystyle\equiv-\bm{\nabla_{\\!t^{\prime}}}A^{\prime}-\partial_{r^{\prime}}\bm{\varphi^{\prime}},$
$\displaystyle\bm{B^{\prime}}$
$\displaystyle\equiv\bm{\nabla_{\\!t^{\prime}}}\times\bm{\varphi^{\prime}}.$
(48)

Then the Maxwell’s equations in superluminal frames take the form:

$\displaystyle\bm{\nabla_{\\!t^{\prime}}}\cdot\bm{E^{\prime}}$
$\displaystyle=-\frac{1}{\varepsilon_{0}c}j^{\prime}$
$\displaystyle\bm{\nabla_{\\!t^{\prime}}}\cdot\bm{B^{\prime}}$
$\displaystyle=0$
$\displaystyle\bm{\nabla_{\\!t^{\prime}}}\times\bm{E^{\prime}}$
$\displaystyle=-\partial_{r^{\prime}}\bm{B^{\prime}}$
$\displaystyle\bm{\nabla_{\\!t^{\prime}}}\times\bm{B^{\prime}}$
$\displaystyle=-\mu_{0}c\bm{\varrho^{\prime}}+\frac{1}{c}\partial_{r^{\prime}}\bm{E^{\prime}}.$
(49)

In the special case of static fields, the Ampère law given by the last
equation (V) reduces to the superluminal form of the "Gauss law":

$\displaystyle\partial_{r^{\prime}}\bm{E^{\prime}}=\frac{\bm{\varrho^{\prime}}}{\varepsilon_{0}}.$
(50)

Our four-dimensional formalism introduced in this work makes the derivation of
the superluminal form of Maxwell’s equations almost effortless.

##  VI Conclusions

We have explicitly shown how to extend special relativity to allow
superluminal inertial observers in $1+3$ dimensional spacetime. Admittedly,
these observers can be physically distinguished from the subluminal ones, but
we show that they are indistinguishable among themselves, just like all the
subluminal inertial observers are. One of the most challenging and counter-
intuitive aspects of our construction is the fact that spacetime metric in
superluminal frames transforms itself into a $3+1$ dimensional one. Therefore
we developed the whole kinematics starting with a derivation of a sensible and
useful definition of velocity in a spacetime involving three-dimensional time
and one-dimensional space. We showed that the speed of light is still
preserved by the superluminal transformations, therefore the task of finding
all possible transformations preserving the speed of light may not be
considered complete with just conventional Lorentz transformations. A subtlety
that is easy to overlook is that superluminal observers require their own,
unorthodox definition of velocity. Furthermore, we characterized expressions
for energy and momentum of superluminal particles and discussed some of their
novel properties, absent in conventional subluminal special relativity.

Finally, we showed that dynamical theories based on a relativistically
invariant principle of least action including superluminal observers cannot be
based on a mechanical paradigm of a classical, point-like trajectory.
Inclusion of superluminal observers leads to a principle of least action based
on a field-theoretical framework with a notable example of the Maxwell’s
theory. Therefore one of the most interesting aspects of extending special
relativity to superluminal frames of reference is the emergence of field
theory. The need for field-theoretic description of matter in terms of "matter
waves" came out of the blue in the early era of quantum theory. Until recently
it was generally believed that postulates underlying quantum theory are
fundamental and cannot be derived from anything more basic. In this work we
showed that the justification of quantum theory using extended relativity,
presented in [1], can be naturally generalized to $1+3$ spacetime and such an
extension leads to the field-theoretic formulation of the quantum theory. This
justifies, or at least provides a plausibility argument, why this extension is
not just an eccentric thought exercise, but reflects something fundamental
about symmetries of laws of physics.

A relevant question that remains to be answered is whether tachyons –
understood as localized lumps of energy traveling with supeluminal speeds –
can physically exist [26]. As was shown here, and earlier in [1], special
relativity does not exclude such a possibility. In quantum field theory,
mechanisms of spontaneous symmetry breaking, such as the Higgs mechanism,
involve fields whose mass squared is negative in the unbroken phase [27], so
excitations of those fields can be regarded as tachyons [28, 29]. In the phase
with broken symmetry, those field are expanded around a degenerate local
minimum, so their excitations have positive mass squared and are are known as
Higgs particles behaving as regular, subluminal particles. However the basic
principle of the mechanism always involves tachyons as a starting point. Given
rather interesting kinematics and dynamics of superluminal particles discussed
in this work, it would be interesting to explore this initial phase before
symmetry breaking in more detail and the role of tachyons therein. We believe
our current work should stimulate new research in this direction.

###### Acknowledgements.

We would like to thank Iwo Białynicki-Birula and Tomasz Miller for useful
comments. K. D. is financially supported by the (Polish) National Science
Center Grant 2021/41/N/ST2/01901. K. T. is partially supported by grant
2018/30/Q/ST9/00795 from the (Polish) National Science Centre.

## References

  * Dragan and Ekert [2020] A. Dragan and A. Ekert, New J. Phys. 22, 033038 (2020).
  * Tolman [1917] R. C. Tolman,  _The theory of the relativity of motion_ , Vol. 44 (University of California Press, 1917).
  * Pirani [1970] F. A. Pirani, Phys. Rev. D 1, 3224 (1970).
  * Parmentola and Yee [1971] J. A. Parmentola and D. D. Yee, Phys. Rev. D 4, 1912 (1971).
  * Landsberg [1972] P. T. Landsberg,  _Time in Statistical Physics and Special Relativity_ (Springer, Berlin, 1972).
  * Goldhaber and Smith [1975] A. S. Goldhaber and J. Smith, Rep. Prog. Phys. 38, 731 (1975).
  * Barrowes [1977] S. C. Barrowes, Foundations of Physics 7, 617 (1977).
  * Basano [1977] L. Basano, International Journal of Theoretical Physics 16, 715 (1977).
  * Maund [1979] J. B. Maund, Found. Phys. 9, 557 (1979).
  * Basano [1980] L. Basano, Found. Phys. 10, 937 (1980).
  * Antippa [1972] A. F. Antippa, Nuov. Cim. A 10, 389 (1972).
  * Parker [1969] L. Parker, Phys. Rev. 188, 2287 (1969).
  * Marchildon  _et al._ [1983] L. Marchildon, A. F. Antippa, and A. E. Everett, Phys. Rev. D 27, 1740 (1983).
  * Dragan [2021] A. Dragan,  _Unusually Special Relativity_ (World Scientific, London, 2021).
  * Bilaniuk  _et al._ [1962] O. M. P. Bilaniuk, V. K. Deshpande, and E. C. G. Sudarshan, Am. J. Phys. 30, 718 (1962).
  * Recami and Mignani [1974] E. Recami and R. Mignani, Riv. Nuovo Cim. 4, 209 (1974).
  * Recami [1986] E. Recami, Riv. Nuovo Cim. 9, 1 (1986).
  * Sutherland and Shepanski [1986] R. I. Sutherland and J. R. Shepanski, Phys. Rev. D 33, 2896 (1986).
  * Demers [1975] P. Demers, Can. J. Phys. 53, 1687 (1975).
  * Ziino [1979] G. Ziino, Phys. Lett. A 70, 87 (1979).
  * Feynman [1987] R. P. Feynman,  _The reason for antiparticles_ (Cambridge University Press, 1987).
  * Terletskii [1968] Y. P. Terletskii,  _Paradoxes in the Theory of Relativity_ (Springer, New York, 1968).
  * Dawe and Hines [1994] R. L. Dawe and K. C. Hines, Aust. J. Phys. 47, 431 (1994).
  * Dawe and Hines [1998] R. L. Dawe and K. C. Hines, Aust. J. Phys. 51, 477 (1998).
  * Dawe and Hines [1992] R. L. Dawe and K. C. Hines, Aust. J. Phys. 45, 591 (1992).
  * Ehrlich [2022] R. Ehrlich, Symmetry 14, 1198 (2022).
  * Srednicki [2007] M. Srednicki,  _Quantum Field Theory_ (Cambridge University Press, 2007).
  * Feinberg [1967] G. Feinberg, Phys. Rev. 159, 1089 (1967).
  * Dhar and Sudarshan [1968] J. Dhar and E. C. Sudarshan, Phys. Rev. 174, 1808 (1968).


# Synthesis of Ballistic Capture Corridors at Mars via Polynomial Chaos
Expansion

M. Liotta111MSc graduate<EMAIL_ADDRESS>Politecnico di
Milano, Via La Masa 34, Milano, 20156, Italy  G. Merisio222Post Doctoral
Research Fellow, Department of Aerospace Science and Technology,
<EMAIL_ADDRESS>Politecnico di Milano, Via La Masa 34, Milano,
20156, Italy  C. Giordano333Post Doctoral Research Fellow, Department of
Aerospace Science and Technology<EMAIL_ADDRESS>AIAA Member. and
F. Topputo444Professor, Department of Aerospace Science and Technology,
<EMAIL_ADDRESS>AIAA Senior Member. Politecnico di Milano, Via La
Masa 34, Milano, 20156, Italy

##  1 Introduction

The space sector is experiencing a flourishing growth. Evidence is mounting
that the near future will be characterized by a large amount of deep-space
missions [1, 2, 3, 4]. In the last decade, CubeSats have granted affordable
access to space due to their reduced manufacturing costs compared to
traditional missions. Nowadays, most miniaturized spacecraft have thus far
been deployed into near-Earth orbits, but soon a multitude of interplanetary
CubeSats will be employed for deep-space missions as well [5]. Nevertheless,
the current paradigm for deep-space missions strongly relies on ground-based
operations [6]. Although reliable, this approach will rapidly cause saturation
of ground slots, thereby hampering the current momentum in space exploration.
At the actual pace, human-in-the-loop, flight-related operations for deep-
space missions will soon become unsustainable.

Self-driving spacecraft challenge the current paradigm under which spacecraft
are piloted in interplanetary space. They are intended as machines capable of
traveling in deep space and autonomously reaching their destination. In
EXTREMA (short for _Engineering Extremely Rare Events in Astrodynamics for
Deep-Space Missions in Autonomy_) [5, 7], these systems are used to engineer
ballistic capture (BC) [8, 9, 10], thereby proving the effectiveness of
autonomy in a complex scenario. The BC mechanism allows capture about a planet
exploiting the natural dynamics, thus without requiring maneuvers [11, 12, 13,
14]. At the expense of longer transfer times, BC orbits are cheaper, safer,
and more versatile from the operational perspective than Keplerian solutions
[11]. Furthermore, BC is a desirable solution for limited-control platforms,
which cannot afford to enter into orbits about a planet due to a lack of
significant control authority. The key is to accomplish low-thrust orbits
culminating in BC. For this, a bundle of BC orbits named ballistic capture
corridor (BCC) can be targeted far away from a planet [15, 16, 17]. Mars is
chosen in this work due to its relevance in the long-term exploration [5].

BC is an extremely rare event, thus massive numerical simulations are required
to find the specific conditions supporting capture. On average, only 1 out of
$10\,000$ conditions explored by algorithms grants capture [18]. To achieve BC
at Mars without any a priori instruction, an inexpensive and accurate method
to construct BCC directly on board is required. Therefore, granting spacecraft
the capability to manipulate stable sets in order to compute autonomously BCC
is crucial [15]. The goal of the paper is to numerically synthesize a corridor
exploiting the polynomial chaos expansion (PCE) technique, thereby applying a
suited uncertainty propagation technique to BC orbit propagation. In the
literature, PCE was introduced and exploited for uncertainty quantification
[19, 20]. PCE can be used as an effective interpolation method, avoiding a
priori definition of interpolant functions but selecting them automatically
starting from the input samples, so that they possess spectral convergence
with respect to the input variables. Furthermore, proving it to be successful
in propagating all-in-once a bundle of trajectories in a deterministic setting
opens the door to a wider use for PCE in data-driven approaches [21]. The
proposed approach is validated against Monte Carlo (MC) simulation. The heavy
computational loads derived by multiple point-wise propagations of BC orbits
when performing guidance and control tasks during the low-thrust
interplanetary cruise are unburdened by the devised methodology. Broadly, this
will facilitate the paradigm shift towards autonomy, thereby favoring the
reduction of mission load on ground stations by decreasing the demand foreseen
in the next years.

The remainder of the paper is organized as follows. In Section 2, the
background is introduced. Then, the proposed methodology follows in Section 3.
Results are presented in Section 4. Eventually, conclusions are drawn in
Section 5.

##  2 Background

###  2.1 Dynamical model

Following the nomenclature in [9], a target and a primary are defined. The
target is the body around which BC is studied. The primary is the main body
around which the target revolves. Target and primary masses are $m_{t}$ and
$m_{p}$, respectively. The mass ratio of the system is
$\mu=m_{t}/(m_{t}+m_{p})$. This work focuses on BC having Mars as target and
the Sun as primary. Reference frames used in this work are the J2000 and
RTN@$t_{i}$ [22].

The precise states of the Sun and the major planets are retrieved from the Jet
Propulsion Laboratory (JPL)’s planetary ephemerides de440s.bsp555Data publicly
available at:
https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de440s.bsp
[retrieved Sep 1, 2023]. (or DE440s) [23]. Additionally, the ephemerides
mar097.bsp of Mars (the target) and its moons are employed666Data publicly
available at: ~/spk/satellites/mar097.bsp [retrieved Sep 1, 2023].. The
following generic leap seconds kernel (LSK) and planetary constant kernel
(PCK) are used: naif0012.tls, pck00010.tpc, and gm_de440.tpc777Data publicly
available at: ~/lsk/naif0012.tls, ~/pck/pck00010.tpc, and ~/pck/gm_de440.tpc
[retrieved Sep 1, 2023]..

The Equations of motion (EoM) of the restricted $n$-body problem are
considered. The gravitational attractions of the Sun, Mercury, Venus,
Earth–Moon (B888Here B stands for barycenter.), Mars (central body), Jupiter
(B), Saturn (B), Uranus (B), and Neptune (B) are taken into account.
Additionally, solar radiation pressure (SRP), Mars’ non-spherical gravity
(NSG), and relativistic corrections [24] are considered. Table 1 collects the
assumed spacecraft parameters needed to evaluate the SRP perturbation. They
are compatible with those of a 12U deep-space CubeSat [25]. Terms of the
infinite series modeling NSG are considered up to degree
$n_{\mathrm{deg}}=$20\text{\,}$$ and order $n_{\mathrm{ord}}=$20\text{\,}$$
[17]. The coefficients to evaluate the NSG perturbation are retrieved from the
MRO120F gravity field model of Mars. Data are publicly available in the file
jgmro_120f_sha.tab, archived in the Geosciences Node of NASA’s Planetary Data
System999Data publicly available at: https://pds-
geosciences.wustl.edu/mro/mro-m-rss-5-sdp-v1/mrors_1xxx/data/shadr/ [retrieved
Sep 1, 2023].. Far from Mars, when in heliocentric motion, the NSG
perturbation is neglected. EoM are integrated in the J2000 inertial frame.

Table 1: Spacecraft parameters for SRP evaluation [25]. Parameter | Unit | Value
---|---|---
Mass–SRP area ratio $m/A$ | $\mathrm{kg}\text{\,}{\mathrm{m}}^{-2}$ | $75\text{\,}$
Coefficient of reflectivity $C_{r}$ | - | $1.3\text{\,}$

The EoM in a non-rotating Mars-centered reference frame are [9, 26, 24]

$\begin{split}\ddot{\mathbf{r}}=&-\frac{\mu_{t}}{r^{3}}\mathbf{r}-\sum\limits_{i\in\mathbb{P}}\mu_{i}\left(\frac{\mathbf{r}_{i}}{r_{i}^{3}}+\frac{\mathbf{r}-\mathbf{r}_{i}}{\left\lVert\mathbf{r}-\mathbf{r}_{i}\right\rVert^{3}}\right)+\frac{QA}{m}\frac{\mathbf{r}-\mathbf{r}_{\odot}}{\left\lVert\mathbf{r}-\mathbf{r}_{\odot}\right\rVert^{3}}-\mathcal{R}\frac{\mu_{t}}{r^{2}}\left(\Lambda\frac{\mathcal{R}^{\top}\mathbf{r}}{r}-\left[J^{\top}\
K^{\top}\ H^{\top}\right]^{\top}\right)+\\\
&+\frac{\mu_{t}}{c^{2}r^{3}}\left[\left(4\frac{\mu_{t}}{r}-v^{2}\right)\mathbf{r}+4\left(\mathbf{r}\cdot\dot{\mathbf{r}}\right)\dot{\mathbf{r}}\right]+2\left(\mathbf{\Omega}\times\dot{\mathbf{r}}\right)+2\frac{\mu_{t}}{c^{2}r^{3}}\left[\frac{3}{r^{2}}\left(\mathbf{r}\times\dot{\mathbf{r}}\right)\left(\mathbf{r}\cdot\mathbf{J}\right)+\left(\dot{\mathbf{r}}\times\mathbf{J}\right)\right]\end{split}$
(1)

where $\mu_{t}$ is the gravitational parameter of the target body (i. e., Mars
in this work); $\mathbf{r}$ and $\dot{\mathbf{r}}=\mathbf{v}$ are the position
and velocity vectors of the spacecraft with respect to the target,
respectively, being $r$ and $v$ their magnitudes; $\mathbb{P}$ is a set of
$n-2$ indexes (where $n$ concerns the $n$-body problem) each one referring to
the perturbing bodies; $\mu_{i}$ and $\mathbf{r}_{i}$ are the gravitational
parameter and position vector of the $i$-th body with respect to the target,
respectively; $A$ is the Sun-projected area on the spacecraft for SRP
evaluation; $m$ is the spacecraft mass; $\mathbf{r}_{\odot}$ is the position
vector of the Sun with respect to the target; $\mathcal{R}$ is the time-
dependent matrix transforming vector components from the Mars-fixed frame to
the non-rotating frame in which the EoM are written; $\Lambda$, $J$, $K$, and
$H$ are defined as in [27];
$c=$299\,792\,458\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$$ (from SPICE
[28, 29]) is the speed of light in vacuum; $\mathbf{J}$ is the rotating
central body’s angular momentum per unit mass in the J2000 frame. Then,
$Q=LC_{r}/\left(4\pi c\right)$ where $C_{r}$ is the spacecraft coefficient of
reflectivity, and $L=S_{\odot}4\pi d_{\mathrm{AU}}^{2}$ is the luminosity of
the Sun. The latter is computed from the solar
constant101010https://extapps.ksc.nasa.gov/Reliability/Documents/Preferred_Practices/2301.pdf
[last accessed Sep 1, 2023].
$S_{\odot}=$1367.5\text{\,}\mathrm{W}\text{\,}{\mathrm{m}}^{-2}$$ evaluated at
$d_{\mathrm{AU}}=$1\text{\,}\mathrm{AU}$$. Lastly,
$\mathbf{\Omega}=\tfrac{3}{2}\dot{\mathbf{r}}_{\odot/t}\times\left(-\mu_{\odot}\mathbf{r}_{\odot/t}\right)/\left(c^{2}r_{\odot/t}^{3}\right)$
where $\mu_{\odot}$ is the gravitational parameter of the Sun;
$\mathbf{r}_{\odot/t}$ and $\dot{\mathbf{r}}_{\odot/t}=\mathbf{v}_{\odot/t}$
are the position and velocity vectors, respectively, of the target body with
respect to the Sun, being $r_{\odot/t}$ and $\dot{r}_{\odot/t}=v_{\odot/t}$
their magnitudes.

The EoM in Eq. (1) are integrated with the GRavity TIdal Slide (GRATIS) tool
[30] in their nondimensional form to avoid ill-conditioning (see normalization
units in Table 2) [9]. Numerical integration is carried out with the DOPRI8
propagation scheme [31, 32]. The dynamics are propagated with relative and
absolute tolerances set to $\mathrm{10^{-13}}$ [9].

Table 2: Nondimensionalization units.

Unit | Symbol | Value | Comment
---|---|---|---
Gravity parameter | $\mathrm{MU}$ | $42\,828.376\text{\,}{\mathrm{km}}^{3}\text{\,}{\mathrm{s}}^{-2}$ | Mars’ gravity parameter $\mu_{t}$
Length | $\mathrm{LU}$ | $3396.0000\text{\,}\mathrm{km}$ | Mars’ radius $R_{\mars}$
Time† | $\mathrm{TU}$ | $956.281\,42\text{\,}\mathrm{s}$ | $(\mathrm{LU}^{3}/\mathrm{MU})^{0\\!.5}$
Velocity | $\mathrm{VU}$ | $3.551\,255\,8\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}$ | $\mathrm{LU}/\mathrm{TU}$

  * †

Time unit chosen such that the nondimensional period of a circular orbit of
radius $\mathrm{LU}$ equals $$2$\pi$.

###  2.2 Ballistic capture corridors

BC orbits are characterized by initial conditions (ICs) escaping the target
when integrated backward and performing $n$ revolutions about it when
propagated forward, neither impacting nor escaping the target. In forward
time, particles flying on BC orbits approach the target coming from outside
its sphere of influence and remain temporarily captured about it. After a
certain time, the particle escapes if an energy dissipation mechanism does not
take place to make the capture permanent. To dissipate energy either a
breaking maneuver or the target atmosphere (if available) could be used [33,
34]. Effects on BC by gravitational attractions of many bodies besides the
primaries and SRP have been investigated in previous works [35, 17, 36, 18].

A particle stability is inferred using a plane in the three-dimensional
physical space [12], this according to the spatial stability definition
provided in [9]. Based on its dynamical behavior, a propagated trajectory is
said to be: i) weakly stable (sub-set $\mathcal{W}_{i}$) if the particle
performs $i$ complete revolutions around the target neither escaping nor
impacting with it or its moons; ii) unstable (sub-set $\mathcal{X}_{i}$) if
the particle escapes from the target before completing the $i$th revolution;
iii) target–crash (sub-set $\mathcal{K}_{i}$) if the particle impacts with the
target before completing the $i$th revolution; iv) moon–crash (sub-set
$\mathcal{M}_{i}$) if the particle impacts with one of the target’s moons
before completing the $i$th revolution; v) acrobatic (sub-set
$\mathcal{D}_{i}$) if none of the previous conditions occurs within the
integration time span. Conditions ii)-v) apply after the particle performs
$(i-$1\text{\,}$)$ revolutions around the target (see Fig. 1). The sub-sets
are defined for $i\in\mathbb{Z}\textbackslash\\{0\\}$, where the sign of $i$
informs on the propagation direction. When $i>0$ ($i<0$) the IC is propagated
forward (backward) in time. A capture set is defined as
$\mathcal{C}_{-1}^{n}:=\mathcal{W}_{n}\cap\mathcal{X}_{-1}$. Therefore, it is
the intersection between the stable set in forward time $\mathcal{W}_{n}$ and
the unstable set in backward time $\mathcal{X}_{-1}$ [9].

Figure 1: Illustration of particle stability. Weakly stable set
$\mathcal{W}_{n}$ in red, unstable set $\mathcal{X}_{n}$ in green, target-
crash set $\mathcal{K}_{n}$ in purple, moon-crash set $\mathcal{M}_{n}$ in
blue, and acrobatic set $\mathcal{D}_{n}$ in yellow. Example with $n=$1$$. The
points on the intersection semi-plane are the ICs.

BCCs are time-varying manifolds supporting capture [15] obtained backward
propagating ICs belonging to a capture sets $\mathcal{C}_{-1}^{n}$, where
$n>0$ is the number of revolutions after capture. They are defined in what
follows. Firstly, a trajectory is defined as:

###### Definition 1

Let $(\mathbf{x}_{0},t_{0})\in\mathbb{R}^{6}\times\mathbb{R}$ and
${\bm{\varphi}}(\mathbf{x}_{0},t_{0};t)$ be the starting point and the
solution at time $t$, respectively, of the state-space representation
$\dot{\mathbf{x}}=\mathbf{f}(\mathbf{x},t)$ of the EoM in Eq. (1). Then, a
trajectory $\gamma$ is defined as
$\gamma(\mathbf{x}_{0},t_{0}):=\\{\bm{\varphi}(\mathbf{x}_{0},t_{0};t)\
\forall t\in\mathbb{R}\\}$. Similarly, backward and forward legs $\gamma_{b}$
and $\gamma_{f}$, respectively, are defined as
$\gamma_{b}(\mathbf{x}_{0},t_{0}):=\\{\bm{\varphi}(\mathbf{x}_{0},t_{0};t)\
\forall t\in[t_{0}-$10$T_{\mars},\,t_{0}]\\}$ and
$\gamma_{f}(\mathbf{x}_{0},t_{0}):=\\{\bm{\varphi}(\mathbf{x}_{0},t_{0};t)\
\forall t\in[t_{0},\,t_{0}+$10$T_{\mars}]\\}$, where
$T_{\mars}=2\pi\sqrt{a_{\mars}^{3}/\mu_{\odot}}=$687\text{\,}\mathrm{days}$$
is the revolution period of Mars, with
$a_{\mars}=$2.279\text{\times}{10}^{8}\text{\,}\mathrm{km}$$ and
$\mu_{\odot}=$1.327\text{\times}{10}^{11}\text{\,}{\mathrm{km}}^{3}\text{\,}{\mathrm{s}}^{-2}$$
the semi-major axis of the Sun–Mars system and the gravitational parameter of
the Sun, respectively.

Sets $\Gamma_{\mathcal{W}_{n}}$, $\Gamma_{\mathcal{X}_{-1}}$, and
$\Gamma_{\mathcal{C}_{-1}^{n}}$ of trajectories $\gamma(\mathbf{x}_{0},t_{0})$
whose ICs $(\mathbf{x}_{0},t_{0})$ belong to weakly-stable set
$\mathcal{W}_{n}$, escape set $\mathcal{X}_{-1}$, and capture set
$\mathcal{C}_{-1}^{n}$, respectively, are
$\Gamma_{\mathcal{W}_{n}}=\\{\gamma(\mathbf{x}_{0},t_{0})\
\forall\,(\mathbf{x}_{0},t_{0})\in\mathcal{W}_{n}\\}$,
$\Gamma_{\mathcal{X}_{-1}}=\\{\gamma(\mathbf{x}_{0},t_{0})\
\forall\,(\mathbf{x}_{0},t_{0})\in\mathcal{X}_{-1}\\}$, and
$\Gamma_{\mathcal{C}_{-1}^{n}}=\\{\gamma(\mathbf{x}_{0},t_{0})\
\forall\,(\mathbf{x}_{0},t_{0})\in\mathcal{C}_{-1}^{n}\\}$. Similarly to a
capture set $\mathcal{C}_{-1}^{n}$, a corridor is designated is defined as
$\mathcal{\check{B}}_{-1}^{n}=\\{\gamma_{b}(\mathbf{x}_{0},t_{0})\
\forall\,(\mathbf{x}_{0},t_{0})\in\mathcal{C}_{-1}^{n}\\}$. An exterior
corridor $\mathcal{\check{E}}_{-1}^{n}$ is a subset of a corridor
$\mathcal{\check{B}}_{-1}^{n}$ including pre-capture trajectories having
heliocentric semi-major axis $a_{\odot}$ greater than the target body’s one
(i. e., Mars, whose semi-major axis
$a_{t}=a_{\mars}=$1.5237\text{\,}\mathrm{AU}$$). It is defined as
$\mathcal{\check{E}}_{-1}^{n}=\\{\gamma_{b}(\mathbf{x}_{0},t_{0})\in\mathcal{\check{B}}_{-1}^{n}:a_{\odot}\left(\bm{\varphi}(\mathbf{x}_{0},t_{0};t)\right)>a_{t}\
\forall t\in[t_{0}-$10$T_{\mars},\,\hat{t}]\\}$ where $\hat{t}<t_{0}$ is a
certain time before capture epoch $t_{0}$ when the escape (or pre-capture) leg
ends in backward time. Contrarily, an interior corridor
$\mathcal{\check{I}}_{-1}^{n}$ is the subset of a corridor
$\mathcal{\check{B}}_{-1}^{n}$ including all trajectories having semi-major
axis smaller than the central body’s one (i. e., Mars). It is defined as
$\mathcal{\check{I}}_{-1}^{n}=\\{\gamma_{b}(\mathbf{x}_{0},t_{0})\in\mathcal{\check{B}}_{-1}^{n}:a_{\odot}\left(\bm{\varphi}(\mathbf{x}_{0},t_{0};t)\right)<a_{t}\
\forall t\in[t_{0}-$10$T_{\mars},\,\hat{t}]\\}$. Consequently,
$\mathcal{\check{B}}_{-1}^{n}=\mathcal{\check{E}}_{-1}^{n}\cup\mathcal{\check{I}}_{-1}^{n}$.

In this work, the interior corridor is of interest because it extends between
Mars and Earth’s orbits. A subcorridor $\mathcal{\check{S}}_{-1}^{n}$, a
generic subset of a corridor $\mathcal{\check{B}}_{-1}^{n}$, is defined as
$\mathcal{\check{S}}_{-1}^{n}=\\{\gamma_{b}(\mathbf{x}_{0},t_{0})\
\forall\,(\mathbf{x}_{0},t_{0})\in\mathcal{D}_{-1}^{n}\\}$, where the generic
domain is
$\mathcal{D}_{-1}^{n}=\\{\mathbf{x}_{0}:\mathbf{x}_{0}\in\mathcal{C}_{-1}^{n}\wedge\mathbf{g}(\mathbf{x}_{0})\leq\mathbf{0}\wedge\mathbf{h}(\mathbf{x}_{0})=\mathbf{0}\\}$,
with $\mathbf{g}(\mathbf{x}_{0})$ and $\mathbf{h}(\mathbf{x}_{0})$ being two
sets of $m\geq 0$ inequality constraints and $n\geq 0$ equality constraints,
respectively, with $m$ and $p$ finite. Finally, the envelope
$\mathcal{\partial\check{S}}_{-1}^{n}$ of a subcorridor
$\mathcal{\check{S}}_{-1}^{n}$ is constructed backward propagating the
subcorridor domain border $\partial\mathcal{D}_{-1}^{n}$. Therefore, it is
defined as
$\partial\mathcal{\check{S}}_{-1}^{n}=\\{\gamma_{b}(\mathbf{x}_{0},t_{0})\
\forall\,(\mathbf{x}_{0},t_{0})\in\partial\mathcal{D}_{-1}^{n}\\}$. An
illustration of prior definitions is proposed in Fig. 2.

Figure 2: Illustration of BCC definitions.

###  2.3 Polynomial chaos expansion

PCE is an uncertain quantification method able to provide an efficient mean
for long term propagation of non-Gaussian distributions. PCE approximates the
stochastic solution of the governing dynamics as a weighted sum of
multivariate spectral polynomials, function of the input random variables [20]

$\bm{x}(t,\bm{\xi})\approx\hat{\bm{x}}(t,\bm{\xi})=\sum_{\bm{\alpha}\in\Lambda_{p,d}}\bm{c}_{\bm{\alpha}}(t)\psi_{\bm{\alpha}}(\bm{\xi})$
(2)

where $\bm{\xi}\in\mathbb{R}^{d}$ is the random input vector,
$\bm{c}_{\bm{\alpha}}(t)$ is the vector of polynomial chaos coefficients
(PCCs), and $\Lambda_{p,d}$ is the set of multi-indices of size $d$ and order
$p$, having so a total dimension equal to [20]

$P=|\Lambda_{p,d}|=\dfrac{(p+d)!}{p!\,d!}.$ (3)

Namely, $\bm{c}_{\bm{\alpha}}(t)\in\mathbb{R}^{n}$, where $n$ is the number of
quantities of interest. If considering the full state propagation, $\bm{x}$
represents both the position and velocity and, hence, $n=6$.

The basis functions $\psi_{\bm{\alpha}}(\bm{\xi})$ are multidimensional
spectral polynomials, orthonormal with respect to the joint probability
measure of the vector $\bm{\xi}$. Thus, the basis functions choice depends
only on the properties of the input variables. For example, Hermite and
Legendre polynomials are the basis for the Gaussian and uniform distributions,
respectively [37].

Generating a PCE means computing the PCCs by projecting the exact solution
onto each basis function $\psi_{\bm{\alpha}}(\bm{\xi})$ [20]

$\bm{c_{\alpha}}(t)=\mathit{E}[\bm{x}(t,\cdot)\psi_{\bm{\alpha}}(\cdot)]=\int_{\Gamma^{d}}\bm{x}(t,\bm{\xi})\psi_{\bm{\alpha}}(\bm{\xi})\rho(\bm{\xi})d\bm{\xi}$
(4)

where $\Gamma^{d}$ is the $d$-dimensional hypercube where the random input
variables are defined. Once the PCCs are computed, the system state associated
to a specific sample within the capture subset can be retrieved at time $t$
effortlessly using Eq. (2). The PCCs estimation methods fall into two
categories: intrusive and non-intrusive. While the first requires laborious
modifications in the governing equations, the latter treat the dynamics as a
black box, thus they are better suited for problems with high-fidelity complex
dynamics [19].

##  3 Methodology

###  3.1 Problem statement

The goal of the BCC synthesis is to produce a numerical approximation of a
subcorridor. The approximation is later made available to the autonomous
guidance and control unit implemented onboard the limited-capability
spacecraft. Ideally, the evaluation of the synthetic subcorridor shall be fast
and inexpensive for spacecraft having limited onboard resources. In
mathematical terms, the general subcorridor synthesis problem can be thus
stated as follows [15]:

###### Problem 1

Numerically synthesize the subcorridor
$\mathcal{\check{S}}_{-1}^{n}=\\{\gamma_{b}(\mathbf{x}_{0},t_{0})\
\forall\,(\mathbf{x}_{0},t_{0})\in\mathcal{D}_{-1}^{n}\\}$ as a function
$\mathbf{x}=\bm{\psi}(\mathbf{p})$ of some parameters $\mathbf{p}$ such that,
given the parameters $\mathbf{p}^{*}$, the state
$\mathbf{x}^{*}=\bm{\psi}(\mathbf{p}^{*})$ is retrieved. In particular, the
state $\mathbf{x}^{*}$ must be targeted by the spacecraft to be temporarily
captured by the central body at time epoch $t_{0}$, so performing at least $n$
revolutions about it.

The choice of the parameters $\mathbf{p}$ to be selected as support of the
subcorridor is of paramount importance, since they define the input for the
synthesis. The aim is to find a set of coordinates for which the regular
capture sub-region is sufficiently large. Indeed, the wider the capture subset
considered, the bigger the region to target in the interplanetary leg. For the
purpose of this work, capture sets at Mars are built following the methodology
depicted in [38]. According to the methodology the initial computational grid
of ICs is bidimensional, thus it is reasonable to use two parameters to define
the subcorridor. The most significant results were achieved in Cartesian and
Keplerian coordinates. To be compliant with the approach discussed in [39],
two components of the Cartesian coordinates (namely, $x$ and $y$) have been
chosen to properly represent the capture set. Additionally, since the BCC is
designed propagating a conveniently selected capture subset, which ensures
ideal post-capture behavior [39, 16], the probability distribution of initial
condition in the capture subset can be considered uniform, i.e., each
condition within the subset boundaries leads to capture. This assumption paves
the way to the use of PCE as an efficient synthetization method, since initial
BC conditions can be treated as stochastic variables. Hence, a PCE-based
subcorridor synthesis problem can be stated as follows:

###### Problem 2

Find the polynomial chaos coefficients $\bm{c}_{\bm{\alpha}}(t)$, such that
$\bm{x}\left(t,\bm{\xi}\right)=\sum_{\bm{\alpha}\in\Lambda_{p,d}}\bm{c}_{\bm{\alpha}}(t)\psi_{\bm{\alpha}}(\bm{\xi})$,
with $\bm{\xi}=(x_{0},y_{0})$ numerically synthesize the subcorridor
$\mathcal{\check{S}}_{-1}^{n}=\\{\gamma_{b}(\mathbf{x}_{0},t_{0})\
\forall\,(\mathbf{x}_{0},t_{0})\in\mathcal{D}_{-1}^{n}\\}$.

###  3.2 Polynomial chaos expansion application

Among the different non-intrusive PCE techniques, given the low number of
input parameters (i.e., $x$ and $y$ coordinates of the capture set ICs), the
pseudospectral collocation approach with full tensor grid [20] is selected. As
a matter of fact, bi-dimensional quadrature schemes (i.e., with $d=2$) suffer
in a limited way of the curse of dimensionality, while pseudospectral
collocation approaches guarantee a simple and accurate method for the PCCs
computation. In this case, PCCs can be computed by solving the stochastic
integral [20]

$\bm{c}_{\bm{\alpha}}(t)\asymp\mathcal{Q}\left[\bm{x}(t,\cdot)\psi_{\bm{\alpha}}(\cdot)\right]=\sum_{i=1}^{d}\sum_{q_{i}=1}^{m_{i}}\bm{x}\left(t,\bm{\xi}_{i}^{q_{i}}\right)\psi_{\bm{\alpha}}(\bm{\xi}_{i}^{q_{i}})\omega_{q_{i}}$
(5)

where $\mathcal{Q}$ is the quadrature integration, $\\{\bm{\xi}_{q}\\}$ is the
set of quadrature nodes and $\\{\omega_{q}\\}$ are the quadrature weights.
Several quadrature rules are available. Gaussian quadrature is exploited in
this work due to their high degree of accuracy [40]. Consequently, nodes
$\bm{\xi}_{q}$ and corresponding weights $\omega_{q}$ are selected depending
on the orthogonal polynomials associated to the probability density function
$\tilde{\rho}(\xi_{i})$ related to each random input $\xi_{i}$. In this case,
since the input $\xi_{i}=(x_{0}^{i},y_{0}^{i})$ is described as a uniform
random variable, the nodes $\\{\xi_{i}^{q_{i}}\\}$ and the weights
$\\{\omega_{q_{i}}\\}$ with $q_{i}=1,\dots,m_{i}$ are the zeros of the
Legendre polynomial of order $m_{i}$ and its quadrature weights.

##  4 Results

###  4.1 Capture subset dimension

Without any loss of generality, the analysis has been conducted on the capture
set depicted in Fig. 3. First, a squared capture subset should be identified
according to the generic domain $\mathcal{D}_{-1}^{n}$ definition provided in
Section 2.2. This is done by identifying one point, which is taken to be the
bottom-left vertex of the subdomain, and a square side length $l$. The bottom-
left vertex of the domain is associated to the following $(r_{0},\omega_{0})$
pair : $r_{0}=R_{t}+$3400\text{\,}\mathrm{km}$$ and
$\omega_{0}=$250\text{\,}\mathrm{deg}$$, where $R_{t}$ is the target (i. e.,
Mars) mean radius. The square side $l$ range is selected to span the cluster
width completely. The lower limit is set to $10\text{\,}\mathrm{km}$, while
the upper one to $500\text{\,}\mathrm{km}$, yielding
$l\in[$10$,$500$]\,$\mathrm{km}$$, as shown in Fig. 3.

Figure 3: Capture set $\mathcal{C}_{-1}^{2}$ at $t_{0}=$ December 9, 2023 at
12:00:00 (UTC). Branches developing exterior BBCs are colored in gray.
Regularity index [10] of ICs belonging to branches developing interior BCCs.
Nondimensional coordinates on the orbital plane
$i_{0}=\Omega_{0}=$0.2$\pi\,$\mathrm{rad}$$ defined in the Mars-centered
RTN@$t_{0}$ frame. Mars is the gray circle with black surround. In the
magnifications, details of chosen capture subsets.

To verify the applicability of PCE to the subcorridor synthesis, the influence
of the subset dimensions on the expansion parameters (i. e., $m_{i}$ and $p$)
is investigated. Firstly, the correlation between the quadrature nodes number
$m_{i}$ and the side length $l$ is assessed. For this analysis, the polynomial
basis order is kept constant to $p=6$. The range [5, 25] is assumed for
$m_{i}$. The outcome is displayed in Fig. 4. The accuracy is evaluated at the
beginning of the pre-capture trajectory (i. e., epoch
$t_{f}=t_{0}-$400\text{\,}\mathrm{days}$$). For $5\leq m_{i}\leq 10$, the
level curves show a ramp-like behavior. Intuitively, the accuracy increases
with the density of quadrature nodes. For instance, consider a side length
$l=$150\text{\,}\mathrm{km}$$, by progressively increasing $m_{i}$ from 5 to
10, it improves by three digits, both in position and velocity. However, for
$m_{i}>10$, a plateau is reached. This is because theoretical results suggests
that, for an increasing number of nodes, the estimation accuracy does not
improve significantly once the estimation is exact. With exact estimation, it
is intended the Gaussian quadrature exactness for polynomial of degree
$p\leq(2m_{i}-1)$ [20]. On the other hand, the method seems to fail for
$l=$300\text{\,}\mathrm{km}$$. The error in position is stuck in the order of
$\mathrm{km}$ and accuracy stops improving with quadrature nodes number.

(a)

(b)

Figure 4: Approximation error at $t_{f}$ as a function of quadrature nodes
number $m_{i}$ and square side length $l$. Polynomial basis order fixed to
$p=6$. Specific subsets having constant $l$ highlighted with dashed horizontal
lines.

Next, the correlation between the polynomial basis order $p$ and the side
length $l$ is investigated. For this analysis, the quadrature nodes number is
kept constant to $m_{i}=10$. Results are shown in Fig. 5. Remarkably, the
method accuracy is greatly affected by the polynomial basis order. Indeed, the
method accuracy improves by progressively increasing $p$. For
$l=$100\text{\,}\mathrm{km}$$, the approximation errors decrease by one order
of magnitude when $p$ increases from 6 to 10. Nevertheless, accuracy is lost
when increasing indefinitely the polynomial basis order. As clearly visible on
the right of plots in Fig. 5, increasing $p$ results in being a disadvantage
for the method accuracy. This behavior may be justified considering that
accuracy of Gaussian quadrature rule is exact for polynomials of degree
$p\leq(2m_{i}-1)$ [20]. The latter translates into the necessity of rising the
quadrature nodes density as $p$ increases. From Fig. 5, it is clear how for
large $p$ values, the quadrature rule fails in estimating accurately the
integral in Eq. (5), thereby growing the approximation error.

(a)

(b)

Figure 5: Approximation error at $t_{f}$ as a function of polynomial basis
order $p$ and square side length $l$. Quadrature nodes number fixed to
$m_{i}=10$.

###  4.2 Simulation parameters

The correlation between $m_{i}$ and $p$ is now investigated by fixing the
square side length of the chosen capture subset to
$l=$10\text{\,}\mathrm{km}$$. In this analysis, the quadrature nodes number
$m_{i}$ and the polynomial basis function $p$ are let vary in ranges [3, 15]
and [4, 12], respectively. The outcome is proposed in Fig. 6. As already
noted, the method accuracy is almost unchanged for large quadrature nodes
numbers. Furthermore, for $m_{i}<9$, accuracy is lost as the polynomial basis
order $p$ increases. The latter confirms what was previously discussed in
Section 4.1.

(a)

(b)

Figure 6: Approximation error at $t_{f}$ as a function of quadrature nodes
number $m_{i}$ and polynomial basis order $p$. Square side length fixed to
$l=$10\text{\,}\mathrm{km}$$.

###  4.3 Polynomial chaos coefficients distribution

The PCCs distribution for several $(m_{i},p)$ pairs, and with
$l=$10\text{\,}\mathrm{km}$$ are shown in Fig. 7. In the plot, dots denote
distributions $1\text{\,}\mathrm{day}$ before capture, while crosses the
distributions $400\text{\,}\mathrm{days}$ before capture (corresponding to
epoch $t_{f}$). Coefficients are normalized with respect to coefficient
$c_{0}$. The higher the decay rate, the more accurate is the expansion result.
This is because higher order terms become less relevant [20].

Results show that for the (10, 6) pair convergence is achieved at both epochs.
Namely, at $t_{0}-$1\text{\,}\mathrm{day}$$ the distribution converges to
within machine precision ($\sim\\!10^{-15}$). Differently, the normalized
coefficients converge at a lower pace to a higher asymptotic value
($\sim\\!10^{-10})$ at $t_{f}$. This is consistent with results presented in
[20], according to which the convergence value can be correlated to the digit
precision in the state estimate. Indeed, numerical errors increase with the
propagation time, leading to lower digit precision. Nonetheless, the digit
precision is very high, both in nondimensional position and velocity,
therefore results are satisfactory.

(a)

(b)

(c)

(d)

Figure 7: Normalized PCE coefficients distributions for 4 out of the 6 state
variables. Dots and crosses indicate coefficients distribution $1\text{\,}$
and $400\text{\,}\mathrm{days}$ before capture, respectively. Note that
$400\text{\,}\mathrm{days}$ before capture corresponds to the final backward
propagation epoch $t_{f}$. Several $(m_{i},p)$ pairs shown.

Varying the quadrature nodes number, the expansion fails in achieving
convergence for $m_{i}=3$. According to results in Fig. 6, $m_{i}>4$ is
recommended. Contrarily, no significant variations are observed between the
(10, 6) pair and the case with $m_{i}=15$. This is consistent with results in
Figs. 4 and 6. Differently, by letting $p$ vary, convergence is assessed for
(10, 4) and (10, 12) pairs. However, all high-order polynomials contribute
marginally to the expansion accuracy for the (10, 12) pair. The latter
representation is usually referred to as sparse [20]. Thus, only the most
significant expansion coefficients are retained when $p=4$, thereby providing
a good estimate of the system state. Finally, the top-right region of Fig. 6
is investigated. Convergence is achieved for the (6, 4) pair. Differently, the
(6, 12) pair fails to converge since the quadrature rule is unsuccessful in
accurately estimating the PCCs for $p>(2m_{i}-1)$ [20].

###  4.4 Validation

The PCE technique is validated against MC analysis. In the following, the
convergence analysis focuses on a specific case study, which develops the
corridor upon the same capture subset with $l=$10\text{\,}\mathrm{km}$$
described in Section 4.1. The selected simulation parameters are $m_{i}=10$
and $p=6$.

The convergence analysis is evaluated at fixed time epoch. Therefore, the
epoch is suitably selected to assess the PCE validity in approximating the
pre-capture trajectory. To this aim, the PCE and MC approaches are compared
considering their evolution in time. Namely, the PCE estimated mean is
compared with the mean estimate computed from $10^{2}$ samples. Samples are
drawn with the Latin hypercube sampling (LHS) technique. In Fig. 8, the
approximation error is evaluated as $\bm{e}(t)=|\bm{c}_{0}(t)-\bm{\mu}(t)|$,
where $\bm{c}_{0}(t)$ is the first PCC, being the mean according to the PCE,
and $\bm{\mu}(t)$ is the mean value retrieved with the MC analysis. From
results, it is observed that the error increases with the propagation time
(from $t_{0}$ to $t_{f}=$-400\text{\,}\mathrm{days}$$). It comes naturally to
assess the method convergence at the final epoch $t_{f}$, which is associated
with the largest approximation error.

(a)

(b)

Figure 8: Method accuracy in estimating the mean. Evolution of first PCC
compared to estimated mean drawn from capture subset with LHS technique and
computed from $10^{2}$ samples.

The convergence is evaluated on the mean and standard deviation. Fig. 9
compares the results for the MC simulation and the PCE technique as a function
of the samples $n$. Convergence errors are computed as

$\varepsilon_{\mu_{r}}=\max_{i=x,y,z}\left|\dfrac{\mu^{\mathrm{MC}}_{i}(t_{f})-\mu^{\mathrm{PCE}}_{i}(t_{f})}{\mu^{\mathrm{PCE}}_{i}(t_{f})}\right|,\qquad\
\varepsilon_{\mu_{v}}=\max_{i=v_{x},v_{y},v_{z}}\left|\dfrac{\mu^{\mathrm{MC}}_{i}(t_{f})-\mu^{\mathrm{PCE}}_{i}(t_{f})}{\mu^{\mathrm{PCE}}_{i}(t_{f})}\right|,$
(6)
$\varepsilon_{\sigma_{r}}=\max_{i=x,y,z}\left|\dfrac{\sigma^{\mathrm{MC}}_{i}(t_{f})-\sigma^{\mathrm{PCE}}_{i}(t_{f})}{\sigma^{\mathrm{PCE}}_{i}(t_{f})}\right|,\qquad\
\varepsilon_{\sigma_{v}}=\max_{i=v_{x},v_{y},v_{z}}\left|\dfrac{\sigma^{\mathrm{MC}}_{i}(t_{f})-\sigma^{\mathrm{PCE}}_{i}(t_{f})}{\sigma^{\mathrm{PCE}}_{i}(t_{f})}\right|.$
(7)

Referring to Fig. 9, the mean and standard deviation convergence errors for
$n=10^{3}$ are $\varepsilon_{\mu_{r}}=$5.682\text{\times}{10}^{-8}$$,
$\varepsilon_{\mu_{v}}=$2.924\text{\times}{10}^{-12}$$,
$\varepsilon_{\sigma_{r}}=$4.943\text{\times}{10}^{-5}$$, and
$\varepsilon_{\sigma_{v}}=$5.4235\text{\times}{10}^{-7}$$. Eventually, the MC
analysis converges to the reference solution with a large degree of accuracy.
Consistently with theoretical results predicted in the literature, the
standard deviation converges at a lower rate with respect to the mean [40].
Therefore, the PCE technique is validated against the MC approach for the
problem at hand. Remarkably, the MC approach with LHS requires the propagation
of $10^{3}$ samples, whereas only $m_{i}^{d}=10^{2}$ propagations are required
by the PCE method when applied to this case study.

(a)

(b)

(c)

(d)

Figure 9: Estimated mean and standard deviation at
$t_{f}=t_{0}-$400\text{\,}\mathrm{days}$$ for 2 out of the 6 state variables.
Results from MC simulation (solid blu line) and PCE (dashed red line) as a
function of samples number. Nondimensional quantities (see Table 2) in the
Mars-centered J2000 frame.

##  5 Conclusion

In this paper, a procedure to accurately and inexpensively synthesize
ballistic capture corridors exploiting the polynomial chaos expansion
technique is discussed and validated against Monte Carlo simulations. Results
prove the convergence of the method, assess the feasibility of ballistic
capture corridor numerical synthesis, and highlight its convenience in terms
of computational efficiency. Remarkably, as the capture subset dimension
increases, the method accuracy is preserved by properly tuning the quadrature
nodes number and the polynomial basis order. For constant polynomial basis
order, the method accuracy improves as the number of quadrature nodes
increases up to the point a plateau is reached. Denser quadrature nodes imply
higher computational costs, reducing the computational efficiency and making
the ballistic capture corridor construction more expensive. On the other hand,
for fixed quadrature nodes number, the method accuracy does not improve by
increasing the polynomial basis order. Indeed, the method accuracy decreases
because the quadrature nodes are insufficient, thereby poorly estimating high-
order polynomials. Overall, results show that a convenient combination of
quadrature nodes number and polynomial basis order improves the accuracy of
the method at limited computational costs.

## Funding Sources

The authors would like to acknowledge the European Research Council (ERC)
since part of this work has received funding from the ERC under the European
Union’s Horizon 2020 research and innovation programme (Grant Agreement No.
864697).

## References

  * Poghosyan and Golkar [2017] Poghosyan, A., and Golkar, A., “CubeSat evolution: Analyzing CubeSat capabilities for conducting science missions,” _Progress in Aerospace Sciences_ , Vol. 88, 2017, pp. 59–83.  10.1016/j.paerosci.2016.11.002.
  * Bandyopadhyay et al. [2016] Bandyopadhyay, S., Foust, R., Subramanian, G. P., Chung, S.-J., and Hadaegh, F. Y., “Review of formation flying and constellation missions using nanosatellites,” _Journal of Spacecraft and Rockets_ , Vol. 53, No. 3, 2016, pp. 567–578.  10.2514/1.a33291.
  * Kalita et al. [2017] Kalita, H., Asphaug, E., Schwartz, S., and Thangavelautham, J., “Network of nano-landers for in-situ characterization of asteroid impact studies,” _68th International Astronautical Congress_ , Adelaide, Australia, 2017.  IAC-17-D3.3.2.
  * Hein et al. [2018] Hein, A. M., Saidani, M., and Tollu, H., “Exploring potential environmental benefits of asteroid mining,” _69th International Astronautical Congress_ , Bremen, Germany, 2018.  IAC-18,D4,5,11,x47396.
  * Di Domenico et al. [2022] Di Domenico, G., Andreis, E., Morelli, A. C., Merisio, G., Franzese, V., Giordano, C., Morselli, A., Panicucci, P., Ferrari, F., and Topputo, F., “The ERC-funded EXTREMA project: Achieving self-driving interplanetary CubeSats,” _Modeling and Optimization in Space Engineering –- New Concepts and Approaches_ , Springer, Cham, Switzerland, 2022, pp. 167–199.  10.1007/978-3-031-24812-2_6.
  * Turan et al. [2022] Turan, E., Speretta, S., and Gill, E., “Autonomous navigation for deep space small satellites: Scientific and technological advances,” _Acta Astronautica_ , Vol. 193, 2022, pp. 56–74.  10.1016/j.actaastro.2021.12.030.
  * Andreis et al. [2022] Andreis, E., Franzese, V., and Topputo, F., “Onboard orbit determination for deep-space CubeSats,” _Journal of Guidance, Control, and Dynamics_ , Vol. 45, No. 8, 2022, pp. 1466–1480.  10.2514/1.G006294.
  * Hyeraci and Topputo [2010] Hyeraci, N., and Topputo, F., “Method to design ballistic capture in the elliptic restricted three-body problem,” _Journal of Guidance, Control, and Dynamics_ , Vol. 33, No. 6, 2010, pp. 1814–1823.  10.2514/1.49263.
  * Luo et al. [2014] Luo, Z.-F., Topputo, F., Bernelli Zazzera, F., and Tang, G. J., “Constructing ballistic capture orbits in the real Solar System model,” _Celestial Mechanics and Dynamical Astronomy_ , Vol. 120, No. 4, 2014, pp. 433–450.  10.1007/s10569-014-9580-5.
  * Dei Tos et al. [2018] Dei Tos, D. A., Russell, R. P., and Topputo, F., “Survey of Mars ballistic capture trajectories using periodic orbits as generating mechanisms,” _Journal of Guidance, Control, and Dynamics_ , Vol. 41, No. 6, 2018, pp. 1227–1242.  10.2514/1.G003158.
  * Topputo and Belbruno [2015] Topputo, F., and Belbruno, E., “Earth–Mars transfers with ballistic capture,” _Celestial Mechanics and Dynamical Astronomy_ , Vol. 121, No. 4, 2015, pp. 329–346.  10.1007/s10569-015-9605-8.
  * Belbruno and Miller [1993] Belbruno, E., and Miller, J., “Sun-perturbed Earth-to-Moon transfers with ballistic capture,” _Journal of Guidance, Control, and Dynamics_ , Vol. 16, No. 4, 1993, pp. 770–775.  10.2514/3.21079.
  * Belbruno and Carrico [2000] Belbruno, E., and Carrico, J., “Calculation of weak stability boundary ballistic lunar transfer trajectories,” _Astrodynamics Specialist Conference_ , Denver, CO, 2000.  10.2514/6.2000-4142, aIAA 2000-4142.
  * Quinci et al. [2023] Quinci, A., Merisio, G., and Topputo, F., “Qualitative study of ballistic capture at Mars via Lagrangian descriptors,” _Communications in Nonlinear Science and Numerical Simulation_ , Vol. 123, 2023, p. 107285.  10.1016/j.cnsns.2023.107285.
  * Merisio [2023] Merisio, G., “Engineering ballistic capture for autonomous interplanetary spacecraft with limited onboard resources,” Ph.D. thesis, Politecnico di Milano, Milan, Italy, 2023.  URL https://hdl.handle.net/10589/196152.
  * Morelli et al. [2022] Morelli, A. C., Merisio, G., Hofmann, C., and Topputo, F., “A Convex Guidance Approach to Target Ballistic Capture Corridors at Mars,” _44th AAS Guidance, Navigation and Control Conference_ , Breckenridge, CO, 2022.  AAS 22-083.
  * Aguiar and Topputo [2018] Aguiar, G., and Topputo, F., “A Technique for Designing Earth–Mars Low-Thrust Transfers Culminating in Ballistic Capture,” _7th International Conference on Astrodynamics Tools and Techniques (ICATT)_ , Oberpfaffenhofen, Germany, 2018.
  * Luo and Topputo [2015] Luo, Z.-F., and Topputo, F., “Analysis of ballistic capture in Sun–planet models,” _Advances in Space Research_ , Vol. 56, No. 6, 2015, pp. 1030–1041.  10.1016/j.asr.2015.05.042.
  * Giordano [2021] Giordano, C., “Analysis, Design, and Optimization of Robust Trajectories for Limited-Capability Small Satellites,” Ph.D. thesis, Politecnico di Milano, Milan, Italy, 2021.  URL http://hdl.handle.net/10589/177695.
  * Jones et al. [2013] Jones, B. A., Doostan, A., and Born, G. H., “Nonlinear propagation of orbit uncertainty using non-intrusive polynomial chaos,” _Journal of Guidance, Control, and Dynamics_ , Vol. 36, No. 2, 2013, pp. 430–444.  10.2514/1.57599.
  * Pugliatti et al. [2023] Pugliatti, M., Giordano, C., and Topputo, F., “The image processing of Milani: challenges after DART impact,” _12th International Conference on Guidance, Navigation and Control Systems (GNC)_ , Sopot, Poland, 2023.
  * Caleb et al. [2022] Caleb, T., Merisio, G., Di Lizia, P., and Topputo, F., “Stable sets mapping with Taylor differential algebra with application to ballistic capture orbits around Mars,” _Celestial Mechanics and Dynamical Astronomy_ , Vol. 134, No. 5, 2022, pp. 1–22.  10.1007/s10569-022-10090-8.
  * Park et al. [2021] Park, R. S., Folkner, W. M., Williams, J. G., and Boggs, D. H., “The JPL Planetary and Lunar Ephemerides DE440 and DE441,” _The Astronomical Journal_ , Vol. 161, No. 3, 2021, p. 105.  10.3847/1538-3881/abd414.
  * Huang et al. [1990] Huang, C., Ries, J. C., Tapley, B. D., and Watkins, M. M., “Relativistic effects for near-Earth satellite orbit determination,” _Celestial Mechanics and Dynamical Astronomy_ , Vol. 48, No. 2, 1990, pp. 167–185.  10.1007/BF00049512.
  * Topputo et al. [2021] Topputo, F., Wang, Y., Giordano, C., Franzese, V., Goldberg, H., Perez-Lissi, F., and Walker, R., “Envelop of reachable asteroids by M-ARGO CubeSat,” _Advances in Space Research_ , Vol. 67, No. 12, 2021, pp. 4193–4221.  10.1016/j.asr.2021.02.031.
  * Dei Tos and Topputo [2019] Dei Tos, D. A., and Topputo, F., “High-fidelity trajectory optimization with application to saddle-point transfers,” _Journal of Guidance, Control, and Dynamics_ , Vol. 42, No. 6, 2019, pp. 1343–1352.  10.2514/1.G003838.
  * Gottlieb [1993] Gottlieb, R. G., “Fast gravity, gravity partials, normalized gravity, gravity gradient torque and magnetic field: Derivation, code and data,” Tech. rep., 1993.  188243, prepared for Lyndon B. Johnson Space Center under contract NAS9-17885.
  * Acton [1996] Acton, C., “Ancillary data services of NASA’s navigation and ancillary information facility,” _Planetary and Space Science_ , Vol. 44, No. 1, 1996, pp. 65–70.  10.1016/0032-0633(95)00107-7.
  * Acton et al. [2018] Acton, C., Bachman, N., Semenov, B., and Wright, E., “A look towards the future in the handling of space science mission geometry,” _Planetary and Space Science_ , Vol. 150, 2018, pp. 9–12.  10.1016/j.pss.2017.02.013.
  * Topputo et al. [2018] Topputo, F., Dei Tos, D. A., Mani, K. V., Ceccherini, S., Giordano, C., Franzese, V., and Wang, Y., “Trajectory design in high-fidelity models,” _7th International Conference on Astrodynamics Tools and Techniques (ICATT)_ , Oberpfaffenhofen, Germany, 2018, pp. 1–9.
  * Montenbruck and Gill [2000] Montenbruck, O., and Gill, E., _Satellite Orbits Models, Methods and Applications_ , Springer, Heidelberg, Germany, 2000.  10.1007/978-3-642-58351-3, pages 117–156.
  * Prince and Dormand [1981] Prince, P. J., and Dormand, J. R., “High order embedded Runge–Kutta formulae,” _Journal of Computational and Applied Mathematics_ , Vol. 7, No. 1, 1981, pp. 67–75.  10.1016/0771-050x(81)90010-3.
  * Luo and Topputo [2021] Luo, Z.-F., and Topputo, F., “Mars orbit insertion via ballistic capture and aerobraking,” _Astrodynamics_ , Vol. 5, No. 2, 2021, pp. 167–181.  10.1007/s42064-020-0095-4.
  * Giordano and Topputo [2022] Giordano, C., and Topputo, F., “Aeroballistic Capture at Mars: Modeling, Optimization, and Assessment,” _Journal of Spacecraft and Rockets_ , 2022, pp. 1–15.  10.2514/1.A35176.
  * Merisio and Topputo [2021] Merisio, G., and Topputo, F., “Characterization of ballistic capture corridors aiming at autonomous ballistic capture at Mars,” _2021 AAS/AIAA Astrodynamics Specialist Conference_ , Big Sky, Virtual, 2021.  AAS 21-677.
  * Luo and Topputo [2017] Luo, Z.-F., and Topputo, F., “Capability of satellite-aided ballistic capture,” _Communications in Nonlinear Science and Numerical Simulation_ , Vol. 48, 2017, pp. 211–223.  10.1016/j.cnsns.2016.12.021.
  * Xiu and Karniadakis [2002] Xiu, D., and Karniadakis, G. E., “The Wiener–Askey Polynomial Chaos for Stochastic Differential Equations,” _SIAM Journal on Scientific Computing_ , Vol. 24, No. 2, 2002, pp. 619–644.  10.1137/S1064827501387826.
  * Luo [2020] Luo, Z.-F., “The role of the mass ratio in ballistic capture,” _Monthly Notices of the Royal Astronomical Society_ , Vol. 498, No. 1, 2020, p. 1515–1529.  10.1093/mnras/staa2366.
  * Merisio and Topputo [2022] Merisio, G., and Topputo, F., “An algorithm to engineer autonomous ballistic capture at Mars,” _73rd International Astronautical Congress_ , Paris, France, 2022.  IAC-22,C1,9,10,x73057.
  * Xiu [2010] Xiu, D., _Numerical methods for stochastic computations: a spectral method approach_ , Princeton University Press, Princeton, NJ, 2010.  10.1515/9781400835348, Chap. 1, 4.