# Relativistic GPS in 3-dimensions

S.Carlonia,b, L.Fatibenec,d,e, M.Ferrarisc, R.G.McLenaghane, A.Orizzontec,d

a DIME Sez. Metodi e Modelli Matematici, Università di Genova, (Italy).
b Institute of Theoretical Physics Faculty of Mathematics and Physics, Charles
University (Czech Republic)
c Department of Mathematics, University of Torino (Italy)
d INFN - Sezione Torino (Italy)
e Department of Applied Mathematics, University of Waterloo (ON-Canada)

###### Abstract

We extend to three dimensions the proposal of a completely relativistic
positioning system (rPS). The system does not rely on approximations, in fact,
it works at a few Schwarzschild radii from a black hole, and it does not rely
on Newtonian physics or special relativity. Since general relativity (GR)
claims to be our fundamental framework to describe classical physics, it must
provide tools to bootstrap physics within the theory itself, without relying
on previous approximated frameworks.

The rPS is able to self-diagnose, that is, it detects deviations from
assumptions about the gravitational field and consequently stops operations;
in addition it is robust, i.e., it is able to autonomously restore operations
when assumptions are restored.

From a more general viewpoint, the rPS is equivalent to geodesy in spacetime,
which establishes a (conventional) coordinate system on a surface by means of
measurements within the surface itself, as well as allowing it to extract
information about the intrinsic geometry of the same surface. In other words,
the positioning system is potentially able to extract information about the
gravitational field (which in fact is identified with the geometry of
spacetime) in addition to the gravitational theory, which describes its
dynamics. Thus, it becomes a framework within which one can operationally
distinguish different theories of gravitation.

##  1 Introduction

In a previous paper [1] we proposed a setting for a toy model relativistic
positioning system (rPS), which is able to determine the position of a user
(in spacetime), together with the orbital parameters of the constellation of
satellites (up to isometries) as well as to check the assumptions made about
the gravitational field and other forces acting on the satellites. The idea is
not new since there are many similar though not identical proposals; (see [2],
[3], [4], [5], [6], [7], [8], [9]).

The whole project 1s based on a rather naive initial question, namely, why
should one need a ground control, also named control segment in the GPS
literature, when, in principle, one could use the GPS infrastructure itself to
determine the position of satellites’ constellation? We also aim to
investigate a system which can function in a strong field regime (namely
something which can work well at a few multiples of the Schwarzschild radius
of a compact object) producing, in principle, exact results. Moreover, we aim
for a fully relativistic system that does not resort in any way to Newtonian
or special relativistic physics. We desire a system which does not assume
structures fixed on spacetime other than a given gravitational field, and
requires no external synchronizations, no external validation, and no
calibration.

Of course, in [1], we did some theoretical simplifications: first of all, we
assumed satellites to carry proper clocks and to be free falling and,
secondly, we restricted to 1+1 dimensions, namely a system of one spatial
dimension. As a result, the system presented was a toy model, a proof of
concept for a system that needed further development.

In 1+1 dimension, unlike in dimension 2+1, one has no orbits around a compact
object, all satellites either fly away or into the central compact objects, in
the second case soon disappearing over an event horizon. The problem is also
greatly simplified in view of the fact that any light ray emitted, sooner or
later reaches any satellite moving on the same side. In 1+1 there is no
control problem to be solved, namely one has only to select the initial
direction to hit the satellite. The analogous problem in 2+1 is much more
complicated since the satellite is moving and the gravitational field, at the
same time, is bending the light ray trajectories.

Working in dimension 2+1, we show that the toy model can be extended and that
these problems do not impact in a fundamental way on an rPS. At a minimum, we
collect mathematical methods to treat the problem in GR, in all its glory. We
present a rPS which does not make any weak field approximation and does not
require an external synchronization or calibration, in other words, it does
not need a control segment. It determines all orbital parameters of the space
segment, namely the satellites and the clocks they carry under the assumption
of a static, rotationally symmetric gravitational field. The system is able to
diagnose deviations from prior assumptions (e.g., the free falling of
satellites, or about the gravitational field), and notify the users to prevent
ill positioning. The system is also able to autonomously restore operations
when these anomalies disappear. In principle, it can also model deviations
perturbatively by regarding some perturbations as part of the variables to fit
together with the orbital parameters of the space segment.

The procedure is not, in principle, different from the one used in geodesy:
one measures relative distances (or angles) between a number of points fixed
on a surface (e.g., the surface of the Earth) from which measurements one can
extract the intrinsic geometry of the surface, e.g., determine the curvature
of the Earth’s surface from intrinsic measurements [10]. If one of the points
is the position of the user and the surface is replaced by spacetime, that is
exactly what an rPS does; it determines the position of the user with respect
to the other points in the constellation as well as the intrinsic geometry of
the surface (which in the case of spacetime is the gravitational field).

A similar problem has been studied by other authors from a more theoretical
viewpoint; (see [11], [12]). In fact, a rPS is a recognized way of measuring
the gravitational field and defining observables in GR, which in vacuum is
well known to have no non-trivial observables.

Moreover, by claiming GR to be a fundamental theory, specifically, our
framework to understand classical physics, one expects to able to provide,
within GR, a way of bootstrap fundamental experiments to describe classical
phenomena. The first step is to be able to describe the motion of test
particles in a spacetime neighbourhood without relying on a Newtonian (or SR)
scheme in which coordinates are endowed with a direct meaning of distances,
and one selects a class of observers (namely inertial observers) which gives
us a description of physics which defines structures on spacetime (namely,
affine structure or absolute time and space) which we know to be unphysical.

A relativistic positioning system is a necessary step that enables us to
describe the motion of test particles around an observer, that is, a first
step towards the extraction of the laws of motion without reliance on the
Newtonian ones that we know are only approximately true.

The paper is organized as follows: in Section 2, we review how one may
describe analytically a satellite constellation. In Section 3, we consider the
description of light rays, in particular how to find a ray emitted from a
satellite which hits a target event.

In Appendix A, we discuss the structure of the light conoid, i.e., the “light
conoid” in spacetime (as opposed to the light cone in the tangent space). This
study allows us to divide light rays into classes depending on whether they
hit the BH or are simply scattered, or whether they go around the BH clockwise
or counterclockwise and, if needed, to divide each ray into branches which can
be suitably parameterized. We note that each class corresponds to a region on
the light conoid which also exactly accounts for all details of strong lensing
[13].

In Appendix B, we present a symplectic framework to describe light rays which
is equivalent to the Lagrangian formalism used in Section 3. However, it
provides a framework to solve in general the control problem of finding a
solution of a Hamiltonian system passing through an initial and a final
position.

In Section 4, we set up a rPS and discuss how one can parameterize the signals
received at a target event $p_{0}$ on a satellite $\chi_{0}$. In Section 5, we
give an example of the simulation of the model in the case of a rPS around a
BH in the strong regime. In Section 6, we similarly describe the parameters of
a rPS around the Earth, with orbital parameters similar to the actual NAVSTAR-
GPS. These two cases present different challenges, and although they are
particular cases of a general framework, they can be optimized in different
directions. Finally, we briefly discuss conclusions and perspectives in
Section 7.

##  2 Satellite constellation

We have to deal with two problems: describing the orbits of satellites
$\chi_{i}$ and the light rays exchanged between them. We start with the first
problem, since the second will be a modification of it. We consider the
Schwarzschild metric in dimension 2+1 which, with respect to the standard
coordinate system $(t,r,\theta)$, has the form

$g=-A(r)c^{2}dt^{2}+B(r)dr^{2}+r^{2}d\theta^{2},$ (1)

where $A(r)=1-\hbox{$\alpha\over r$}$, $B=\hbox{$1\over A$}$, $\alpha$ is a
constant related to the mass of the central object, and $c$ denotes the speed
of light in the limit $r\rightarrow\infty$. This is the most general
stationary, spherically symmetric solution of the vacuum Einstein equations.
Here, unlike the 1+1 case, we have real bounded orbits, which we can extend to
the 3+1 situation, by restricting the constellation (and the user) to be on
the equatorial plane.

The basic idea is to avoid approximations and instead, use the Lagrangian
formalism to reduce the problem to Weierstrass equations which can be
integrated by semi-analytic methods. We use the parameter invariant Lagrangian
for describing timelike geodesics, since it is well known that it describes
the worldlines of test particles and, since it is invariant with respect to
parameterizations, it can also be used later for light rays on which, of
course, we cannot define proper time. The invariant arc-length associated with
the Lagrangian defines the geodesic distance between two points in spacetime.
The standard definition of the _distance-function_ $\Gamma$ is the square of
the geodesic distance with the sign according to whether the non-light-like
geodesic is time-like or space-like. (See [14], [15], [16], [17], [18], [19],
[20], [21], [22]). Synge [23] calls $\Omega=\hbox{$1\over 2$}\Gamma$ the world
function.

He argues that one can deduce from it the geometric structure of spacetime,
which makes sense since the geometric structure is determined by the
gravitational field via Eienstein’s equations, and the gravitational field, as
any other field, is encoded and observable only through the effects it has on
test particles, the motion of which is described by the invariant Lagrangian.

We can choose any parameterization gauge, but here we shall initially use
coordinate time $t$, not the proper time; we shall try to highlight why later
on. Hence the Lagrangian for timelike worldlines is

$\tilde{L}=\sqrt{Ac^{2}(t^{\prime})^{2}-{(r^{\prime})^{2}\over
A}-r^{2}(\theta^{\prime})^{2}}\>ds$ (2)

where primes denote derivatives with respect to the arbitrary parameter $s$.
Choosing $t$ as the parameter we obtain

$L=\sqrt{Ac^{2}-\hbox{$\dot{r}^{2}\over A$}-r^{2}\dot{\theta}^{2}}\>dt,$ (3)

where the dots denote derivatives with respect to $t$.

Since the Lagrangian $L$ does not depend explicitly on $t$ or on $\theta$ we
have the following two first integrals:

${c^{2}}\epsilon=-{{c^{2}}A\over\sqrt{Ac^{2}-\hbox{$\dot{r}^{2}\over
A$}-r^{2}\dot{\theta}^{2}}},\qquad\qquad
k=-{r^{2}\dot{\theta}\over\sqrt{Ac^{2}-\hbox{$\dot{r}^{2}\over
A$}-r^{2}\dot{\theta}^{2}}},$ (4)

where $\epsilon$ and $\theta$ are constant. Solving for the Lagrangian
velocities we get

$\dot{r}^{2}={A^{2}\over
r^{2}\epsilon^{2}}((\epsilon^{2}c^{2}-A)r^{2}-Ak^{2}),\qquad\qquad\dot{\theta}^{2}={k^{2}A^{2}\over
r^{4}\epsilon^{2}}.$ (5)

These can be cast into the form of Weierstrass equations

$\displaystyle\left({dr\over dt}\right)^{2}=$ $\displaystyle{A^{2}\over
r^{2}\epsilon^{2}}((\epsilon^{2}c^{2}-A)r^{2}-Ak^{2})=:\Phi(r;\epsilon,k)$ (6)
$\displaystyle\left({dr\over d\theta}\right)^{2}=$ $\displaystyle{r^{2}\over
k^{2}}((\epsilon^{2}c^{2}-A)r^{2}-Ak^{2})=:\Psi(r;\epsilon,k)$

The functions $\Phi$ and $\Psi$ are called Weierstrass functions or also
effective potentials. The value of $(\epsilon,k)$ parameterizes the solutions
since they can be obtained from initial conditions. Both of them diverge to
$-\infty$ for light rays. The motion of the system takes place where $\Phi\geq
0$, which is called the allowed region.

We remark that $\epsilon$ and $k$ are not exactly physical energy and angular
momentum (as dimensional analysis shows directly). If one wants to map first
integrals to physical quantities, one can compare them to what happens in the
case of Keplerian motions.

The next step will be to set $A=1-\hbox{$\alpha\over r$}$, i.e. the
Schwarzschild solution rather than working on a generic static, spherically
symmetric metric.

If one uses $A=1-\hbox{$\alpha\over r$}+\lambda r^{2}$, one would have AdS-
Schwarzschild. This case will be left for future investigation, e.g., to
determine how long we should observe the system to detect a given non-zero
$\lambda$.

By the substitution $A=1-\hbox{$\alpha\over r$}$ into the definitions, the
Weierstrass functions take the following form:

$\displaystyle\Phi(r;\epsilon,k)=$ $\displaystyle{(r-\alpha)^{2}\over
r^{5}\epsilon^{2}}((\epsilon^{2}c^{2}-1)r^{3}+\alpha
r^{2}-k^{2}r+k^{2}\alpha)$ (7) $\displaystyle\Psi(r;\epsilon,k)=$
$\displaystyle{r\over
k^{2}}((\epsilon^{2}c^{2}-1)r^{3}+r^{2}\alpha-k^{2}r+k^{2}\alpha)$

Figure 1: we show the qualitative allowed region of polynomial $p(r)$ for
massive particles and $\epsilon^{2}c^{2}-1>0$. The dashed line is at the
Schwarzschild radius. The red thin line is the limiting case with an
asymptotic circular orbit.
A thick blue line on the left shows one allowed region; the material point
falls into the horizon.
The thick blue line on the right shows two allowed regions; the material point
either falls into the horizon or is scattered to infinity.

The allowed regions are determined by the sign of the polynomial

$p(r)=(\epsilon^{2}c^{2}-1)r^{3}+\alpha r^{2}-k^{2}r+k^{2}\alpha$ (8)

we shall focus on the external region $r>\alpha$. For bounded orbits, there
should be an allowed region $[r_{-},r_{+}]$ in the external region, which
implies that $r=r_{\pm}>\alpha$ be a simple root of the polynomial,
$p(r_{\pm})=0$. Since this is a polynomial of degree 3, it has at least one
and at most three real roots. Note that
$p(\alpha)=\alpha^{3}c^{2}\epsilon^{2}\geq 0$. As $r\rightarrow\infty$ ,
$p(r)$ diverges to infinity with a sign given by the leading coefficient
$\epsilon^{2}c^{2}-1$.

If one has $\epsilon^{2}c^{2}-1>0$ (see Fig.1), the polynomial $p(r)$ has
either 0 or 2 roots in the external region. If there are zero roots, the whole
exterior region is allowed, and the test particle either escapes to infinity
or falls into the BH. If there are two roots, one has two allowed regions, one
near the horizon where the test particle falls and one near infinity where the
test particle escapes to infinity. In both cases there are no bounded orbits.
If $\epsilon^{2}c^{2}-1=0$, the leading term becomes $r^{2}$, which implies
either 0 or 2 solutions, and in both cases no bounded orbits.

Figure 2: Qualitative allowed region of polynomial $p(r)$ for massive
particles and $\epsilon^{2}c^{2}-1<0$. The dashed line is at the Schwarzschild
radius. The thin red line is the limiting case with an asymptotic circular
orbit.
The thick blue line on the left shows one allowed region; the material point
falls into the horizon.
The thick blue line on the right shows two allowed regions; the material point
either falls into the horizon or goes along a bounded orbit.

If $\epsilon^{2}c^{2}-1<0$ (see Fig.2), the polynomial has either 1 or 3 roots
in the external region. Only one solution means an allowed region near the
horizon, across which all particles fall. which again implies no bounded
orbits. We are left with the last option (3 real roots) which yields two
allowed regions: one close to the horizon across which test particles fall and
the other which describes bounded orbits. As a consequence of the above
analysis we are lead to consider the case of three roots where we want to
recast the polynomial $p(r)$ in the form:

$(\epsilon^{2}c^{2}-1)(r-r_{0})(r-r_{1})(r-r_{2})$ (9)

By comparison, we obtain three equations for the unknowns
$(\epsilon^{2},k^{2},r_{0},r_{1},r_{2})$ which can be solved for
$(\epsilon^{2},k^{2},r_{0})$ as functions of $(r_{1},r_{2})$. Since expressing
the roots as functions of $(\epsilon^{2},k^{2})$ is difficult, we prefer to
write $(\epsilon^{2},k^{2})$ (and $r_{0}$) as a function of the other two
roots which form the boundary of the allowed region
$[r_{1}=r_{-},r_{2}=r_{+}]$ associated to the bounded orbit. In this way, we
are parameterizing solutions by perihelion (perigee) $r_{-}$ and aphelion
(apogee) $r_{+}$, instead of $(\epsilon^{2},k^{2})$ or initial conditions.
This is convenient since $r_{\pm}$ are also integration limits of improper
integrals one which may be used to determine the orbital period. In this case
the appropriate Weierstrass functions in terms of $(r_{-},r_{+})$ are given by

$\displaystyle\Phi(r;r_{-},r_{+})=$
$\displaystyle{c^{2}\alpha(r-\alpha)^{2}(r_{+}-r)(r-r_{-})\over(r_{+}-\alpha)(r_{-}-\alpha)(r_{+}+r_{-})r^{5}}\left((r_{+}r_{-}-\alpha(r_{+}+r_{-}))r-\alpha
r_{+}r_{-}\right)$ (10) $\displaystyle\Psi(r;r_{-},r_{+})=$
$\displaystyle{r(r_{+}-r)(r-r_{-})\over
r_{+}^{2}r_{-}^{2}}\left((r_{+}r_{-}-\alpha(r_{+}+r_{-}))r-\alpha
r_{+}r_{-}\right)$

We can utilize any computer algebra program to evaluate the integral, namely

$t(r;r_{+},r_{-})=\int_{r_{-}}^{r}{dR\over\sqrt{\Phi(R;r_{+},r_{-})}}$ (11)

We stress that Maple or Mathematica can evaluate this integral analytically
precisely because the polynomial $p(r)$ is factorized, which results from the
fact that the orbits are parameterized with respect to their aphelion and
perihelion values $r_{\pm}$.

We also remark that even though we started with a parameterization in terms of
the relative time $t$, we eventually get orbits parameterized by $r$. We argue
that we do not need or should try to resist this. After all, we are interested
in the orbits, and the parameterization is just a means to describe them.

In order to follow the satellites for as many orbits as we wish, the branches
which describe an orbit must be glued together. Thus we require a standardized
way of denoting the different branches.

Hence, we start at the perihelion $r=r_{-}$ at $t=t_{0}$. On branch 0, $r$
increases to aphelion $r_{+}$ in time $t(r;b=0)=t_{0}+t(r;r_{+},r_{-})$ which
it reaches at time $t_{1}=t_{0}+T$, where we set $T:=t(r_{+};r_{+},r_{-})$ for
half of the $t$-period. Branch 1 starts at aphelion $r_{+}$ at $t_{1}$ and $r$
decreases to $r_{-}$. Branch 1 is described as
$t(r;b=1)=t_{1}+T-t(r;r_{+},r_{-})$. It reaches $r=r_{-}$ at time
$t_{2}=t_{0}+2T$. Repeating the process for branch 2 we have $r$ increasing
from $r=r_{-}$ to $r=r_{+}$ in time $t(r;b=2)=t_{2}+t(r;r_{+},r_{-})$. It ends
the branch at time $t_{3}=t_{2}+T=t_{0}+3T$. On branch 3, $r$ decreases from
$r_{+}$ to $r_{-}$ with $t(r;b=3)=t_{3}+T-t(r;r_{+},r_{-})$. It ends when
$r=r_{-}$ and $t_{4}=t_{3}+T=t_{0}+4T$. This procedure may be continued. We
also have branch -1 in which $r$ decreases from $r_{+}$ to $r_{-}$ where
branch 0 starts. During the branch $t(r;b=-1)=t_{0}-t(r;r_{+},r_{-})$.
Accordingly, branch -1 starts at $r_{+}$ and $t_{-1}=t_{0}-T$. Similarly,
during branch -2, $r$ increases from $r_{-}$ to $r_{+}$ with
$t(r;b=-2)=t_{0}-2T+t(r;r_{+},r_{-})$. And so on.

At the same time, we can determine the angle $\theta$ by

$\theta(r;\theta_{0},r_{+},r_{-})=\theta_{0}\pm\int_{r_{-}}^{r}{dR\over\sqrt{\Psi(R;r_{+},r_{-})}}$
(12)

For definiteness we shall consider counterclockwise orbits, which correspond
to the $+$ sign on outgoing (even) branches ($\theta$ increasing with $r$) and
the $-$ on ingoing (odd) branches ($\theta$ increases as $r$ decreases).
During a branch, the angle increases by the amount

$\Theta=\int_{r_{-}}^{r_{+}}{dR\over\sqrt{\Psi(R;r_{+},r_{-})}}$ (13)

If this quantity is $\pi$ (as in Keplerian motions), the orbit is closed. The
deviation $\delta=2(\Theta-\pi)$ measures the precession in radians per orbit.

For the proper time, we have

$d\tau=c^{-1}\sqrt{Ac^{2}(t^{\prime})^{2}-\hbox{$1\over
A$}(r^{\prime})^{2}-r^{2}(\theta^{\prime})^{2}}\>ds=c^{-1}\sqrt{\hbox{$Ac^{2}\over\Phi$}-\hbox{$1\over
A$}-\hbox{$r^{2}\over\Psi$}}\>dr$ (14)

Thus, we can define

$\tau(r;\tau_{0},r_{+},r_{-})=\tau_{0}\pm
c^{-1}\int_{r_{-}}^{r}\sqrt{\hbox{$Ac^{2}\over\Phi$}-\hbox{$1\over
A$}-\hbox{$R^{2}\over\Psi$}}\>dR$ (15)

This result allows us to keep track of proper time along branches so that
everything is locally parameterized in $r$. By this means, we can compute the
orbit of a satellite in spacetime back and forth for as many branches as we
need. In addition, the whole orbit of a satellite is determined by
$(t_{0},\theta_{0},r_{+},r_{-})$. The proper time is known once we set
$\tau_{0}$, namely the proper time at perihelion on the branch $b=0$. See the
Maple drawing of the orbit in Figure 3.

BH Earth

Figure 3: BH. We set $\alpha=1$ and $c=1$. Orbital parameters are
$r_{-}=20\alpha$, $r_{+}=28\alpha$, $t_{0}=10c^{-1}\alpha$,
$\theta_{0}=\hbox{$7\pi\over 6$}$, and $\tau_{0}=-c^{-1}\alpha$. The
continuous line is the (exact analytical) satellite orbit for the branches -1,
0, and 1. The thick part is branch 0. The cylinder at the center is the
Schwarzschild event horizon at $\alpha=1$. The orbit precesses by
$\delta=2\Theta-2\pi=0.4477$ (radians per orbit). The precession can be
computed to any desired precision (e.g. in this case, we have
$\delta=0.4476574054738341735420832219320682$ which is clearly more than we
can hope to measure). The orbital period is $2T=1117.0454$ (in units fixed by
setting $c=1$) in relative time, $2\tau=1081.5289$ in satellite proper time.
Hence, we see that the satellite clock slows down (with respect to coordinate
time) by a factor of $0.9682$ as an effect of the gravitational field of the
BH. Again we can compute these quantities to any desired precision.

Earth. Here we use SI units. We set Earth parameters $\alpha=0.008870355\>m$
and $c=299792458\>ms^{-1}$. Orbital parameters are $r_{-}=24289000\>m$,
$r_{+}=25089000\>m$, $t_{0}=-1600\>s$, and $\theta_{0}=-\hbox{$\pi\over 6$}$
$\tau_{0}=-1600\>s$, which are comparable to those of a GPS satellite. The
continuous line is the (exact analytical) satellite orbit for the branches
from -1 to 1. The thick part is branch 0. The cylinder at the center of the
orbit is the worldsheet of Earth’s surface at $3189\>km$. The orbital
precession is $\delta=2\Theta-2\pi\simeq 3.4\cdot 10^{-9}$ (radians per
orbit). The orbital period is $2T=38606.364226049\>s$ (i.e. about $10.7h$,
which is comparable with what is expected for NAVSTAR-GPS satellites) in
Earth’s relative time, $2\tau=38606.364215646\>s$ in satellite proper time.
The slow down of proper time cannot be observed at this level of precision.
However, the absolute slowing down can be computed to be $2T-2\tau=1.040\times
10^{-5}\>s$. Hence, we see that in these conditions the satellite precession
is small, and the satellite clock slows down by about $10\>\mu s$ per orbit.

It is not difficult to add further satellites (two more, in fact) on similar
orbits. For the sake of simplicity, we select the orbital parameters so that
the orbits never cross, each satellite being in a ring disjoint from the
others, starting from $\chi_{0}$ in the smaller ring, $\chi_{1}$ being in the
middle, and $\chi_{2}$ being the furthest out. In this way, we have a
constellation of 3 satellites, which (at least when they do not eclipse each
other) is a minimal constellation as discussed in [1].

Before considering light rays, we remark that if we had utilized the quadratic
Lagrangian, which is known to provide geodesics parameterized by proper time,
we would, accordingly, be looking for solutions in the form $\tau(r;r_{\pm})$.

In view of our solution, we can now guess how to produce the solution
parameterized by proper time: we should have to write $\tau(r(t);r_{\pm})$,
where the function $r(t)$ would need to be obtained by solving our solution
$t=t(r;r_{\pm})$. In our case, though, the function $t(r;r_{\pm})$ is
explicitly computed, and we see it is a (rather complicated) combination of
elliptic functions. The inverse cannot be easily expressed in terms of special
functions, as it would not be for the solutions of the quadratic Lagrangian.
For this reason, it is convenient to use relative time rather than proper time
(and eventually the coordinate $r$), even if proper time may have a better
physical meaning.

##  3 Light rays description

We now turn to the determination of the orbits of light rays in the
Schwarzschild spacetime. This task is relatively easy, since we are able to
use the formalism developed in the previous section to determine the satellite
orbits. This is possible because for the earlier calculation we used the
parameter invariant Lagrangian (2) instead of the more usual quadratic one
based on proper time parameterization.

It is well known (see [24], [25], [26]) that light rays can be approached
arbitrarily closely by the orbits of material particles. Thus, we can obtain a
description of light rays as the limit of a sequence of orbits closer and
closer to the light conoid. Regarding terminology , we call light cone with
vertex $p_{0}$ the subset of the tangent space at $p_{0}$ defined by the
lightlike lines through $p_{0}$ , while we call light conoid with vertex
$p_{0}$, the envelope of lightlike worldlines through $p_{0}$ as a subset of
spacetime. An event $p$ belongs to the light conoid with vertex $p_{0}$ iff
there exists a light ray connecting $p$ to $p_{0}$.

Both the first integrals $\epsilon$ and $k$ diverge for light rays, as one can
see from equations (4). However, the quotient $K=\hbox{$k\over\epsilon$}$ is
finite. As a consequence, the light rays are obtained by replacing
$k=K\epsilon$ in the Weierstrass equations (7), then letting
$\epsilon\rightarrow-\infty$. This procedure gives the Weierstrass equations
directly for light rays as

$\displaystyle\left({dr\over dt}\right)^{2}=$
$\displaystyle{(r-\alpha)^{2}\over
r^{5}}(c^{2}r^{3}-K^{2}r+K^{2}\alpha)=:\Phi(r;K)$ (16)
$\displaystyle\left({dr\over d\theta}\right)^{2}=$ $\displaystyle{r\over
K^{2}}(c^{2}r^{3}-K^{2}r+K^{2}\alpha)=:\Psi(r;K)$

We can now use the techniques that we employ for the satellites with one less
parameter. The real roots of the Weierstrass functions are governed by the
polynomial $p(r)=c^{2}r^{3}-K^{2}r+K^{2}\alpha$. They can be $1$ or $3$ real
roots (see Fig.4).

Figure 4: Qualitative allowed region for the polynomial $p(r)$ for light
rays.The dashed line is at the Schwarzschild radius. The red thin line is the
limiting case of an asymptotic circular orbit at $r=\hbox{$3\over 2$}\alpha$.
The thick blue line on the left illustrates one possible allowed region where
the light ray crosses the horizon.
The thick blue line on the right illustrates two allowed regions: the light
ray either crosses the horizon or it scatters to infinity.

The case of 1 real root occurs when $0\leq K^{2}\leq\hbox{$27\over
4$}\alpha^{2}c^{2}=:K^{2}_{cr}$, in this case the root $r_{0}$ satisfies the
inequality $r_{0}\leq-3\alpha$. The outer region $(\alpha,+\infty)$ is
entirely allowed, which means that a light ray either falls into the BH or
escapes to infinity. These light rays, called in-falling rays, form a single
branch, either ingoing or outgoing. The worldline of an in-falling light ray
in spacetime is defined by the parameterized curve

$t(r;t_{0},K)=t_{0}+e\int_{r_{0}}^{r}{dR\over\sqrt{\Phi(R;K)}}\qquad\theta(r;\theta_{0},K)=\theta_{0}+\sigma
e\int_{r_{0}}^{r}{dR\over\sqrt{\Psi(R;K)}}$ (17)

where $e=\pm 1$ determines whether the branch is ingoing ($e=-1$) or outgoing
($e=+1$), and $\sigma=\pm 1$ denotes whether the light ray goes around the BH
clockwise ($\sigma=-1$) or counterclockwise ($\sigma=+1$).

In the boundary case $K^{2}=K^{2}_{cr}$, one has an asymptotic goal for
$r=r_{cr}=\hbox{$3\over 2$}\alpha$ which corresponds to a double root of
$p(r)$. That is a single ray, we ignore it since the integrals diverge in this
case.

In the case of 3 distinct real solutions (i.e. $K^{2}>K^{2}_{cr}$), it is
convenient to write the polynomial in factored form:
$p(r)=c^{2}r^{3}-K^{2}r+K^{2}\alpha=c^{2}(r-r_{1})(r-r_{m})(r+r_{1}+r_{m})$,
with $0<r_{1}\leq r_{m}$ where $r_{m}$ is the minimal approach radius of the
light ray. In this case, we can express $K^{2}$ and $r_{1}$ as a function of
$r_{m}$, which becomes the only parameter. Thus we can express the Weierstrass
functions as $\Phi(r;r_{m})$ and $\Psi(r;r_{m})$. It follows that there are
two allowed regions: one is confined in the region $\alpha\leq r\leq
r_{1}<\hbox{$3\over 2$}\alpha$, the other is $r\in[r_{m},\infty)$ which
corresponds to a light ray that arrives from and escapes to spatial infinity.
Accordingly, such rays are called scattered rays.

A scattered light ray is entirely determined by $r_{m}$. It consists of two
branches, an ingoing and an outgoing branch. By writing the Weierstrass
functions in terms of $r_{m}$ we can use Maple or Mathematica to integrate the
equations analytically

$\displaystyle t(r;t_{0},r_{m})=$ $\displaystyle
t_{0}+e\int_{r_{0}}^{r}{dR\over\sqrt{\Phi(R;r_{m})}}$ (18)
$\displaystyle\theta(r;\theta_{0},r_{m})=$ $\displaystyle\theta_{0}+\sigma
e\int_{r_{0}}^{r}{dR\over\sqrt{\Psi(R;r_{m})}}$

in terms of elliptic functions as we did for the satellites.

It follows that we can fix an event $p_{0}$, choose a value for $r_{m}$, and
draw exactly one light ray through $p_{0}$. The light ray reaches a minimum
$r=r_{m}$ at some point on its worldline in the past or in the future of
$p_{0}$ and escapes again to spatial infinity. In Appendix $A$ we use Maple to
draw explicitly the exact past conoid through an event $p_{0}$. The surface we
obtain accounts for all multiple images of a particle that an observer at
$p_{0}$ can observe around the BH (namely, its strong lensing effect[13]). We
also see that the surface is regular in a neighbourhood of $p_{0}$ except at
$p_{0}$ itself, namely, it is locally a conoid. Further away from the vertex,
the surface becomes multi-sheeted and self-intersecting, which is precisely
what determines the strong lensing effect. An important feature of the light
conoid we obtain is that we are able to draw it exactly and analytically.

To summarize the situation, for each satellite on a bounded orbit, we have
$(\theta_{0},t_{0})$ the initial position given at the perihelion $r_{-}$ on
branch 0, as well as $(r_{-},r_{+})$ to identify the orbit. Since the
satellite carries a proper clock, we also need an initial value $\tau_{0}$ for
the proper time. Then one can compute the half-periods $T$, $\Theta$, and
$\tau$. Each orbit then splits into a countable number of branches
parameterized by $b\in\mathbb{Z}$, alternating between ingoing (odd $b$) and
outgoing (even $b$) half-orbits.

For a light ray, we have in-falling (one for each sub-critical $K^{2}$) and
scattered (one for each $r_{m}>\hbox{$3\over 2$}\alpha$) ray. In-falling light
rays can be ingoing or outgoing and clockwise or counterclockwise. Each in-
falling light ray is made of a single branch.

On the other hand, each scattered light ray consists of two branches, one
ingoing and one outgoing. They are classified by the minimal approach
parameter $r_{m}$, and by whether the target is on the ingoing or outgoing
branch and whether they go around the BH clockwise or counterclockwise.

However, for positioning, we need instead to solve a control problem. Imagine
we fix an event $p$ (possibly, but not necessarily, on a satellite) and
consider a(nother) satellite (namely, a generic timelike worldline
$(t_{s}(r),r,\theta_{s}(r))$, $s$ standing for satellite (or source of the
message). We need to find a light ray through $p$ and some point $p_{i}$ of
the source worldline corresponding to some value of the parameter $r=r_{i}$.
If this is possible, we can determine the proper time $\tau_{i}$ of the clock
on the satellite at emission.

It is now clear how to do this. The satellite worldline is described by

$(r=r_{s},\theta_{s}(r_{s};\theta_{0},r_{\pm}),t_{s}(r_{s};t_{0},r_{\pm}))$
(19)

while the light ray through the target $p=(r_{t},\theta_{t},t_{t})$ is
described by

$(r,\theta_{l}(r;t_{t},\theta_{t},r_{m}),t_{l}(r;r_{t},t_{t},r_{m}))$ (20)

(or, equivalently, by
$(r,\theta_{l}(r;r_{t},\theta_{t},K^{2}),t_{l}(r;r_{t},t_{t},K^{2}))$ for an
in-falling ray). Thus a point $r_{s}$ on both the ray and the source worldline
must obey the equations

$\theta_{l}(r_{s};t_{t},\theta_{t},r_{m})=\theta_{s}(r_{s};\theta_{0},r_{\pm})+2n\pi\qquad
t_{l}(r_{s};r_{t},t_{t},r_{m})=t_{s}(r_{s};t_{0},r_{\pm})$ (21)

even though the first does not depend on $t_{t}=t_{0}$ the second does not
depend on $\theta_{t}=\theta_{0}$. We have two technical complications to
overcome: there are different types of rays, each of which is differently
parameterized. Indeed, the signal source, in the case of scattered rays, could
be on the same branch of the target or on the other. These possibilities yield
3 (in-falling, scattered on the same branch, scattered on the opposite branch)
$\times 2$ (ingoing or outgoing at the target) $\times 2$ (counterclockwise or
clockwise) light rays, for a total of 12 types of light rays, each of which is
characterized by an analytic expression that can be examined for possible
intersections with the source worldline. However, the parameterizations of the
trajectories are complicated combinations of elliptic functions which are
difficult to be solved analytically. Although we have analytic descriptions of
rays and orbits, from the control problem on, we need to proceed numerically,
which introduced some concerns about precision of computed values that need to
be considered along the way. For these reasons, the search for the
intersection is a bit convoluted. In practice, it is convenient to eliminate
the angle overcounting (encoded by the integer $n\in\mathbb{Z}$) by solving
numerically two of the equations among

$\cos(\theta_{l}(r_{s}))=\cos(\theta_{s}(r_{s}))\qquad\sin(\theta_{l}(r_{s}))=\sin(\theta_{s}(r_{s}))\qquad
t_{l}(r_{s})=t_{s}(r_{s})$ (22)

We then verify that the 3rd is also satisfied. This approach eliminates
spurious solutions which correspond to the right cosine and the wrong sine, or
vice versa.

We fix the target (so that $(r_{t},\theta_{t},t_{t})$ become known parameters)
and the source satellite (so that also $(\theta_{0},t_{0},r_{\pm})$ become
known parameters). Thus, we are left with two unknowns to solve for,
$(r_{s},r_{m})$ for the scattered light rays, $(r_{s},K^{2})$ for the in-
falling rays. The solution value $r_{m}$ (or $K^{2}$) determines the ray,
while $r_{s}$ determines the point of intersection of the ray and the
satellite worldline, that is, the source event.

We should also mention that, on scattered rays, it is convenient to compare
the values $(\theta_{m},t_{m})$ of the minimal approach event on the two
branches of a ray passing through the source and the target so that one treats
the two scattered cases on an equal footing (and also because prolonging a ray
to another branch accumulates numerical errors).

If a solution is not found, one needs to explore the previous branch of the
source orbit, as happens in particular when the source satellite is near its
aphelion or perihelion. Since the search is carried out on a type by type
basis for rays, if a solution is not found on a specific type, one should
search on the other types.

One can speed up the search by noticing that $r_{t}>r_{s}$ and following
clockwise the source along its orbit. We cycle among the rays of the following
types: in-falling-counterclockwise-outgoing, scattered on the same branch-
counterclockwise-outgoing, scattered on other branch-counterclockwise-
outgoing, scattered on other branch-clockwise-outgoing, scattered on the same
branch-clockwise-outgoing, in-falling, clockwise-outgoing.

Analogously, if $r_{s}>r_{t}$, we cycle among the rays of the following types;
in-falling-counterclockwise-ingoing, scattered on the same branch-
counterclockwise-ingoing, scattered on other branch-counterclockwise-outgoing,
scattered on other branch-clockwise-outgoing, scattered on the same branch-
clockwise-ingoing, in-falling-clockwise-ingoing.

Two types of rays are not realized with all counterclockwise satellite orbits
(scattered on other branch-clockwise-ingoing and scattered on other branch-
counterclockwise-ingoing). Either way, we need to explore only six types of
rays in each case instead of 12.

Accordingly, we can start from the event on the source orbit, which has the
same $t$ as the target, and then retrace backwards the orbit looking for
solutions. If we find a solution, we may continue to look for a better
solution (one with a greater $r_{m}$). Note that changing the orientation of
the rotation corresponds to the satellites passing behind the BH, and,
consequently, to an event of the strong regime. Strong field events are harder
to compute, although they carry more information about the system, since the
emission time is more sensitive to the characteristic of the system.

The cases are undoubtedly complicated to treat thoroughly. However, we
eventually, find several intersections corresponding to messages sent from
satellites to the events, each of which has its own emission proper time
$\tau_{i}$.

Depending on the situation, we can be satisfied with one message from each
satellite (e.g., the one with a greater $r_{m}$) or decide to consider
multiple messages of the same satellites. The choice depends essentially on
how close a ray can go to the horizon and still be eventually detected. In an
Earth setting, any ray hitting the surface at $r\sim 7\times 10^{7}\alpha$
cannot be later detected. In a BH setting (with no accretion disk at
$r\sim\alpha$) we can approach the horizon very closely and be bent at will.
Consideration of satellites around the Earth is easier, but redundancy in the
number of satellites is required to overcome the fact that the central mass
can eclipse a pair of satellites. This situation cannot arise in a BH setting,
since one can use strong lensing to avoid eclipses.

In any event, once we have determined the first two messages, which can be
referred to as first generation messages, their emission is an event itself
(on a satellite’s wordlines), and one can look back for messages sent there
from other satellites. We call these second generation messages, which, are 4
in number. It follows that, given each generation with $2^{g}$ emission
events, we can go back to the previous generation’s $2^{g+1}$ messages. It is
thus evident that this number increases exponentially with the number of
generations.

There is also an equivalent formulation of the control problem based instead
on the Hamiltonian formalism. We presented it in Appendix B, since it is a
more general formulation of the control problem for a general Hamiltonian
system, which can be utilized whenever we have a complete integral of the
Hamilton-Jacobi equation.

##  4 Determining constellation geometry

We set up the positioning system in a region around the central mass by
considering 3 freely falling satellites, each one equipped with an atomic
clock. They are able to exchange signals and mirror the signals they receive
to the other satellites. We assume that the signals contain information about
the clock reading at emission. Accordingly, at some point, a satellite
receives a set of signals that have bounced back and forth among the
satellites. In principle, each satellite receives an infinite sequence of
proper times, organized in generations.

We still need to discuss how we can extract, out of the infinite sequence of
proper times $\tau_{i}$ available at a base event $p_{0}$ on a satellite,
e.g., on satellite $\chi_{0}$, enough information to determine all orbital
parameters of the satellite constellation, possibly together with some
gravitational parameters (as, e.g., $\alpha$ or perhaps some parameterization
of the perturbations of the gravitational field to account for deviations from
the Schwarzschild metric). We would also like to determine some additional
orbital parameters (e.g., to model non-purely free falling satellites, which
take into account the forces acting on the satellites due to solar wind or
nuclear reactor thermal leaks as occurred with Pioneer’s acceleration
anomaly), and even some parameters for the modification of the gravitational
theory (as, for example, happens if we use Palatini $f({\hbox{\cal
R}})$-theories and employ Taylor coefficients of the function $f({\hbox{\cal
R}})$ expansion to parameterize the theory).

One can regard these situations as discretizing the classical degrees of
freedom for a generic gravitational field in a generic gravitational theory.
If we keep the number of unknown parameters to be determined finite, there is
no real difference between positioning, measuring the gravitational field, or
discussing the observability of different gravitational theories.

We reason backwards. The positioning of a user can easily be determined once
the positions at the emission of satellites are known.

The user receives a signal from each satellite which contains information
about the satellites’ positions (in spacetime) at emission. By knowing the
position of the satellites in spacetime, the user may conclude that he is on
the future light conoid at the emission position. In dimension 2+1, this is a
2-dimensional surface, and the user has three such surfaces. The intersection
of 3 surfaces in a 3d manifold is a discrete set of points. If we add
redundancy (e.g. increasing the number of satellites, restricting the domain
of the chart, or remembering some information from previous positioning), the
user can determine his position in spacetime, up to an isometry.

We remark that the user does not need to use the clock readings it receives
with the signal. This means that we do not need the clocks on satellites to be
synchronized in any way at this point.

Of course, for this scenario to work, each satellite should be able to
determine its position in spacetime. This requires that the satellites be able
to exchange signals.

Figure 5: The first four generations of signals exchanged by three satellites
in dimension 3, considering a base event on $\chi_{0}$. The event $p_{A}$ is
on satellite $\chi_{A\>mod\>3}$. Thus the signal from $p_{19}$ to $p_{9}$ is
sent from satellite $\chi_{1}$ to $\chi_{0}$.

Suppose that a satellite receives the sequence

$(\tau_{0};\tau_{1},\tau_{2};\tau_{3},\tau_{4},\tau_{5},\tau_{6};\tau_{7},\dots,\tau_{14};\tau_{15}\dots)$
(23)

which is the data from which it needs to determine the orbital parameters of
all satellites. On the other hand, the orbital parameters correspond to the
degrees of freedom of the satellite constellation. As we discussed in Section
2, the orbit of each satellite $\chi_{i}$ is described by by five parameters
$(t_{i},\theta_{i},\tau_{i},r^{i}_{\pm})$, namely its initial position
$(t_{i},\theta_{i})$ on its current orbit branch and clock reading at
perihelion $\tau_{i}$, as well as the perihelion and aphelion of the orbit
$r_{\pm}$.

If the clock reading $\tau_{e}$ at emission is given, then these orbital
parameters completely determine the satellite’s position at emission. For
example, given its orbital parameters, one can solve
$\tau(r_{e};r_{\pm})=\tau_{e}$ to obtain $r_{e}$ at emission. Then
$(r_{e},\theta(r_{e};r_{\pm}),t(r_{e};r_{\pm}))$ is the position at emission.
In other words, each satellite can compute its position by reading its clock,
once the orbital parameters are known.

However, because of isometries, we cannot obtain all of the orbital parameters
from the signals. Since all signals are left invariant by rotations and time
translations, we can fix the gauge by choosing $t_{0}=0$ and $\theta_{0}=0$,
i.e., giving all other positions with respect to the position of $\chi_{0}$.
We can also decide to make absolute clock readings irrelevant by setting the
initial condition $\tau_{i}$ so that the proper time of $\chi_{i}$ at its
first event $p_{i}$, is $\tau_{i}=0$. By doing this, we make any initial
synchronizations irrelevant.

A particular sequence is entirely determined by the ($10=15-3-2$) orbital
parameters as well as the $r_{0}$ of the base event on $\chi_{0}$. Thus we
have 11 parameters as unknowns to parameterize the possible sequences of
signals received at $p_{0}$.

Given the orbital parameters $(r_{0};\theta_{i},t_{i},r^{i}_{\pm})$, we have
the base event
$p_{0}=(r_{0},\theta(r_{0};\theta_{0},r^{0}_{\pm}),t(r_{0};t_{0},r^{0}_{\pm}))$.
As described in the previous Section, we find the emission events $p_{1}$ on
$\chi_{1}$ and $p_{2}$ on $\chi_{2}$ which can send a light signal to $p_{0}$.
They are determined by $r_{1}$ and $r_{2}$ which completely determine
$p_{k}=(r_{k},\theta(r_{k};\theta_{k},r^{k}_{\pm}),t(r_{k};t_{k},r^{k}_{\pm}))$,
where $r_{k}=r_{k}(r_{0};\theta_{i},t_{i},r^{i}_{\pm})$. Thus we have
determined the first-generation signals.

Given the events $p_{1}$ and $p_{2}$, we can look for signals
$\chi_{0}\rightarrow\chi_{1}$ and $\chi_{2}\rightarrow\chi_{1}$ to
$p_{1}\in\chi_{1}$, which are emitted at $r_{3}$ and $r_{5}$, i.e. at the
events $p_{3}$ and $p_{5}$, as well as signals $\chi_{0}\rightarrow\chi_{2}$
and $\chi_{1}\rightarrow\chi_{2}$ to $p_{2}\in\chi_{2}$, which are emitted at
$r_{6}$ and $r_{4}$, i.e., at the events $p_{6}$ and $p_{4}$. Indeed, if we
write $r_{A}(r_{0};\theta_{i},t_{i},r^{i}_{\pm})$ for each event, we can
express the clock reading at emission as
$\tau_{A}(r_{0};\theta_{i},t_{i},r^{i}_{\pm})=\tau_{A}(r_{A}(r_{0};\theta_{i},t_{i},r^{i}_{\pm};\theta_{i},t_{i},r^{i}_{\pm})$.
In other words, we can proceed backwards at will to compute the
$\tau_{A}(r_{0};\theta_{i},t_{i},r^{i}_{\pm})$ which are specific functions of
the unknowns.

In this way, we obtain an infinite sequence of equations which are invariant
with respect to isometries and clock resettings, namely we have

$\begin{cases}\tau_{3}(r_{0},\theta_{i},t_{i},r^{i}_{\pm})=\tau_{0}-\tau_{3}\cr\tau_{4}(r_{0},\theta_{i},t_{i},r^{i}_{\pm})=\tau_{1}-\tau_{4}\cr\tau_{5}(r_{0},\theta_{i},t_{i},r^{i}_{\pm})=\tau_{2}-\tau_{5}\cr\dots\cr\tau_{A}(r_{0},\theta_{i},t_{i},r^{i}_{\pm})=\tau_{A\>mod\>3}-\tau_{A}=:\Delta\tau_{A}\cr\dots\cr\end{cases}$
(24)

Finally, we can imagine that each satellite broadcasts to the other satellites
its clock reading as well as mirrors the signals which it is receiving at the
same time from the other satellites. From the point of view of a satellite,
e.g., $\chi_{0},$ which receives this transmission, then it is receiving
precisely the sequence $\tau_{A}$ at the base event $p_{0}$. Since all the
information $\tau_{A}$ is available to $\chi_{0}$ at $p_{0}$, it is able to
write down the system (24), solve it and determine the orbital parameters for
all satellites as well as the base point $p_{0}$.

Each satellite, following the same procedure simultaneously, can compute its
position (together with other information, e.g., multiple past positions of
the other satellites and its own). This can be done at any moment by each
satellite which can then broadcast to the user a signal encoding its position
at emission. In this way the user, upon receiving such messages from each
satellite, can determine its own position.

One does not need an unlimited number of equations to determine a finite
number of unknown parameters. Thus, we can select some $N$ initial equations,
and keep the rest of the sequence for later use.

Unfortunately, the system (24) cannot be easily solved; thus, to proceed, we
turn it into a minimization problem for the function

$\chi^{2}(r_{0},\theta_{i},t_{i},r^{i}_{\pm})=\sum^{N+2}_{A=3}(\tau_{A}(r_{0},\theta_{i},t_{i},r^{i}_{\pm})-\Delta\tau_{A})^{2}$
(25)

The solutions of the system (24) correspond to absolute minima of the function
$\chi^{2}$, for which one has $\chi^{2}=0$. This result may be helpful to
exclude local minima which are not of interest here.

A solution $(r_{0},\theta_{i},t_{i},r^{i}_{\pm})$ of (25) not only satisfies
the first $N$ equations in system (24), but also, whether the initial
assumptions we made (about the gravitational field, the motion of the
satellites, and the gravitational theory we used) are indeed satisfied. We
note that, in addition, all the other equations from the system (24) have to
be identically satisfied along the solution.

If this if is not the case, the system will stop functioning as a positioning
system to avoid broadcasting false information to the user. Meanwhile, the
satellites keep exchanging signals, computing orbital parameters, and checking
whether the assumptions are eventually restored. If they are, each satellite
restarts broadcasting data to the user, which is able to determine its
position again. In this sense, the rPS is self-locating and robust.

Alternatively, satellites may test different hypotheses (e.g., transient non-
gravitational forces acting on satellites causing them to deviate from
timelike geodesics or transient deviations from the Schwarzschild metric
caused by, for example, an incident gravitational wave). In this way, the
system can autonomously try to model more general situations and self-
determine extra parameters to model perturbations.

Of course, this is only a concept toy model which needs to be made more
precise. However, it shows how the rPS can be used as a gravitational detector
and as a test for gravitational theories. We also remark that we have been
rather vague about which parameters we fit because we wish to preserve
flexibility. If we are working around the Earth, we can assume that
$\alpha\sim 0.008870\>m$. However, if we dropped out of hyperspace near a
black hole and we set up a rPS there to perform a physical experiment, we
probably would like to measure $\alpha$ as well.

In that case, we would consider $\alpha$ to be an unknown and, if necessary,
increase the number of generations that we use for determining the system
parameters, and fit it with the orbital parameters. We would proceed in a
similar fashion, if we are expecting perturbations of the gravitational field
or want to test some non-standard gravitational theory. The essential idea,
however, is that we split the signals into two sequences: a first sequence
used for fitting and a second containing the remaining signals to be used for
checking the assumptions. The splitting time is called the lookback time,
since only signals after the lookback time are used for the fitting. What
happens before the lookback time is irrelevant to the fitting. Indeed, we are
assuming that deviations from assumptions are negligible after the lookback
time.

##  5 Test the model: a BH case

Using the techniques presented above, we are able to predict (exactly) the
signals generated by a given constellation of satellites. For example, we can
consider the three satellites $\chi_{i}$.

Each satellite carries a proper clock, $\chi_{0}$ has orbital parameters
$r_{+}=10\>\alpha$ and $r_{-}=19\>\alpha$. Its branch 0 starts at
$(t_{0}=0\>\alpha c^{-1},r_{0}=r_{-}=10\>\alpha,\theta_{0}=0)$ with
$\tau_{0}=1$. One can compute that the $t$-coordinate period of $\chi_{0}$ is
$2T_{0}=551.8586\>\alpha c^{-1}$, while its proper period is
$2\tau_{0}=522.2874\>\alpha c^{-1}$, which shows the clock slows down due to
the gravitational field. It precesses by $\delta=0.8736$ radians per orbit.

We stress these quantities are calculated, not measured. Accordingly, we can
provide a prediction as precise as we wish, e.g., with 120 digits. Of course,
a prediction that is too precise to be experimentally confirmed is physically
useless. This does not mean that it does not have other uses. For example,
data are obtained in floating point with a given precision, e.g., 40 digits
and we have no idea on how reliable are the last digits. However, it is simple
to repeat the simulation at 120 digits and check how the digits neighbouring
40 behave. Moreover, we note that, in the case of the Earth, we can discuss
situations ranging from one satellite orbit (about $10^{4}\>s$) down to the
accuracy of an atomic clock, that for the new optical clocks, can be expected
to be less than $10^{-20}\>s$. In other words, the situation we are
considering naturally spans a large number of orders of magnitude.
Consequently, we need to guarantee at least the precision to describe this
situation.

The satellite $\chi_{1}$ has orbital parameters $r_{+}=28\>\alpha$ and
$r_{-}=20\>\alpha$ with $(t_{1}=10\>\alpha c^{-1},\theta_{1}=-\hbox{$7\pi\over
6$},\tau_{1}=-\>\alpha c^{-1})$. One can compute that the coordinate period of
$\chi_{1}$ is $2T_{1}=1117.0454\>\alpha c^{-1}$, while its proper period is
$2\tau_{1}=1081.5289\>\alpha c^{-1}$. It precesses by $\delta=0.4477$ radians
per orbit.

The satellite $\chi_{2}$ has orbital parameters $r_{+}=40\>\alpha$ and
$r_{-}=31\>\alpha$ with $(t_{2}=-10\>\alpha c^{-1},\theta_{2}=\hbox{$3\pi\over
4$},\tau_{2}=-3\>\alpha c^{-1})$. One can compute that the coordinate period
of $\chi_{2}$ is $2T_{2}=1964.3794\>\alpha c^{-1}$, while its proper period is
$2\tau_{2}=1922.4102\>\alpha c^{-1}$. It precesses at $\delta=0.2886$ radians
per orbit.

Then it is just a matter of computation (and the selection of one light ray
exchanged by two satellites when multiple ones are available) to complete the
table for signals described in Figure 2. In this case, we have (in units of
$\alpha c^{-1}$):

Generation 2
$\displaystyle\Delta\tau_{3}=39.17850096193592366660462091592062373$ (26)
$\displaystyle\Delta\tau_{4}=80.11715717097871213365890494429573535$
$\displaystyle\Delta\tau_{5}=35.03574261053054630657492872958108159$
$\displaystyle\Delta\tau_{6}=72.57341977793111680898418627081248044$
Generation 3
$\displaystyle\Delta\tau_{7}=113.17356079751935436022513046580032987$
$\displaystyle\Delta\tau_{8}=35.40991007306522652038724051092347128$
$\displaystyle\Delta\tau_{9}=103.18980669701380595609210925548892818$
$\displaystyle\Delta\tau_{10}=43.31148869147579252763214908406813583$
$\displaystyle\Delta\tau_{11}=65.43360770882868000737676330764741734$
$\displaystyle\Delta\tau_{12}=125.55127652440722347053216752374541211$
$\displaystyle\Delta\tau_{13}=80.13407626061348502470439015672952679$
$\displaystyle\Delta\tau_{14}=112.24603172686866256539758944644350801$
Generation 4
$\displaystyle\Delta\tau_{15}=159.84959178824144499228158441696501744$
$\displaystyle\Delta\tau_{16}=113.17357459679319449200407586632423227$
$\displaystyle\Delta\tau_{17}=93.10407540385622828546311154167323011$
$\displaystyle\Delta\tau_{18}=86.43058203591576981983712540155117890$
$\displaystyle\Delta\tau_{19}=113.52586678836403966572945824855443803$
$\displaystyle\Delta\tau_{20}=143.91779026987465003621630041257875864$
$\displaystyle\Delta\tau_{21}=103.52001932858460948977835457701477422$
$\displaystyle\Delta\tau_{22}=141.74429768374661292012973911416795685$
$\displaystyle\Delta\tau_{23}=76.85574698554547517515015464726496156$
$\displaystyle\Delta\tau_{24}=125.56902461398873586013930329151515120$
$\displaystyle\Delta\tau_{25}=185.53542787116689390175923961978786119$
$\displaystyle\Delta\tau_{26}=112.26226941709415618622735669007918004$
$\displaystyle\Delta\tau_{27}=130.22536702210812661237475328349493992$
$\displaystyle\Delta\tau_{28}=136.84337647699999488633948073154505905$
$\displaystyle\Delta\tau_{29}=113.45578988671430465571931586496739842$
$\displaystyle\Delta\tau_{30}=172.61747721153048259619184107188535649$
$\displaystyle\dots$

This is the observed sequence, the shown digits of which, are a truncation of
the actual computed values to 35 decimal places. We remark that these are
observables of the strictest kind in GR; they correspond to coincidences in
the theory, of the sort that one can use to model observations at the most
fundamental level.

It remains to show that we can actually solve for the orbital parameters, and
determine the precision and accuracy of the solution. This is an issue
precisely because we do not have much control over the expressions of the
functions $\tau_{A}(r_{0};\theta_{i},t_{i},r^{i}_{\pm})$ or on their
differentiability. We shall return to these questions for the case of the
Earth considered below. We continue our analysis by determining how the
function $\chi^{2}$ behaves near a solution. We need to guarantee that the
solution is an isolated minimum of $\chi^{2}$ and that the function is not too
flat around it so that the minimizing algorithm can converge to the solution.

As an initial test, perturb the unknowns by a Gaussian noise signal
($\sigma=10^{-10}$) and check to see how the function $\chi^{2}$ behaves. We
observe that, as expected, it responds by $\Delta\chi^{2}\sim 10^{-17}$. If we
perturb only the base event ($\sigma=10^{-10}$) on $\chi_{0}$, then we get
$\Delta\chi^{2}\sim 10^{-21}$. If we perturb ($\sigma=10^{-10}$) only
satellite $\chi_{0}$, $\Delta\chi^{2}\sim 10^{-16}$. If we perturb
($\sigma=10^{-10}$) only satellite $\chi_{1}$, $\Delta\chi^{2}\sim 10^{-17}$.
If we perturb ($\sigma=10^{-10}$) only satellite $\chi_{2}$,
$\Delta\chi^{2}\sim 10^{-17}$. If we perturb ($\sigma=10^{-10}$) only one
parameter of a satellite, we get from $\Delta\chi^{2}\sim 10^{-16}$ to
$\Delta\chi^{2}\sim 10^{-23}$ depending on the perturbation and on the
parameter we decide to perturb (besides, it obviously depends on the
unperturbed configuration we start from).

Figure 6: The signals of generation 1 on the left, of generation 1 and 2 on
the right.

The response of the function $\chi^{2}$ to a perturbation, in this case, is
essentially meaningless, because we are using homogeneous units that depend on
the central mass. However, in preparation for the Earth case, we need to know
the precision required to compute the function to determine the parameters. In
the case of a BH, for example, we see that if we do not guarantee at least a
precision to $10^{-20}$ in the value of $\chi^{2}$, we will eventually be
unable to say very much. For example, in this case, we need a precision in
measuring the signals $\tau_{A}$ at least of the order $10^{-10}$.

##  6 Test the model: an Earth case

In the case of Earth, we model a situation similar to the actual GPS. All
results are understood to be in SI units.

We consider 3 satellites $\chi_{i}$, each of which carries a proper clock.
$\chi_{0}$ has orbital parameters $r_{-}=2.3289\cdot 10^{7}\>m$ and
$r_{+}=2.4089\cdot 10^{7}\>m$. At perihelion of branch 0 we have
$\theta_{0}=\hbox{$\pi\over 6$}$, $t_{0}=-1500\>s$ and $\tau_{0}=-1500\>s$.

The coordinate period of $\chi_{0}$ is $2T_{0}=36284.71748955\>s$, while its
proper period (if we want to be able to see the difference) is
$2\tau_{0}=36284.71747936\>s$. Comparison of these results shows that the
clock slows down due to the gravitational field. It precesses by
$\delta=3.5301\cdot 10^{-9}$ radians per orbit.

The satellite $\chi_{1}$ has orbital parameters $r_{-}=2.4289\cdot 10^{7}m$
and $r_{+}=2.5089\cdot 10^{7}m$. For this satellite we have
$\theta_{1}=-\hbox{$\pi\over 6$}$, $t_{1}=-1600\>s$ and $\tau_{1}=-1600\>s$.
The coordinate period of $\chi_{1}$ is $2T_{1}=38606.36422605\>s$, while its
proper period is $2\tau_{1}=38606.36421565\>s$. It precesses by
$\delta=3.3871\cdot 10^{-9}$ radians per orbit.

The satellite $\chi_{2}$ has orbital parameters $r_{-}=2.5289\cdot 10^{7}m$
and $r_{+}=2.6089\cdot 10^{7}m$. Further we have $\theta_{2}=0$,
$t_{2}=-1400\>s$ and $\tau_{2}=-1400\>s$. The coordinate period of $\chi_{2}$
is $2T_{1}=40975.51797465\>s$, while its proper period is
$2\tau_{2}=40975.51796404\>s$. It precesses by $\delta=3.2551\cdot 10^{-9}$
radians per orbit. In this case, we have (all in seconds):

Generation 2 $\displaystyle\Delta\tau_{3}=0.177670545858865720368199359727407$
(27) $\displaystyle\Delta\tau_{4}=0.006161739629485806051906548644238$
$\displaystyle\Delta\tau_{5}=0.057840923371459509637356702512431$
$\displaystyle\Delta\tau_{6}=0.125990984860198315939016097030856$ Generation 3
$\displaystyle\Delta\tau_{7}=0.064002699645967673093551408614620$
$\displaystyle\Delta\tau_{8}=0.177670272908770755554022606444972$
$\displaystyle\Delta\tau_{9}=0.183831819371646273448205999862825$
$\displaystyle\Delta\tau_{10}=0.177670391080652676094777905921607$
$\displaystyle\Delta\tau_{11}=0.125990791303975115694999501565061$
$\displaystyle\Delta\tau_{12}=0.183832280120454500258023229898599$
$\displaystyle\Delta\tau_{13}=0.125990875102861017455161133288578$
$\displaystyle\Delta\tau_{14}=0.064002666904713921376083935866949$ Generation
4 $\displaystyle\Delta\tau_{15}=0.241673189747825261701984791678735$
$\displaystyle\Delta\tau_{16}=0.183831659226008694553003310344026$
$\displaystyle\Delta\tau_{17}=0.183831536956154698351275644437245$
$\displaystyle\Delta\tau_{18}=0.355340782158789752780184396885066$
$\displaystyle\Delta\tau_{19}=0.183832125100949090748893101070209$
$\displaystyle\Delta\tau_{20}=0.121843663566318426499010983039103$
$\displaystyle\Delta\tau_{21}=0.303660984816986754999415382305216$
$\displaystyle\Delta\tau_{22}=0.132152610754713273885450011424911$
$\displaystyle\Delta\tau_{23}=0.235511427015118169710635646895444$
$\displaystyle\Delta\tau_{24}=0.303661311202605558019309110571608$
$\displaystyle\Delta\tau_{25}=0.070164447082984281655998663824688$
$\displaystyle\Delta\tau_{26}=0.183831878295795158020429178386318$
$\displaystyle\Delta\tau_{27}=0.251981582606611857425290985243845$
$\displaystyle\Delta\tau_{28}=0.183832119974415537678739505974830$
$\displaystyle\Delta\tau_{29}=0.183831997704255089293768828935526$
$\displaystyle\Delta\tau_{30}=0.189993553438716616270862526046800$
$\displaystyle\dots$

We stress that these are the results of one simulation for more or less
generic parameters. To solve the positioning problem we need to sample
parameter space with a significant number of simulations. The results are
discussed below.

We begin by repeating the perturbation testing. We perturb the unknows by a
Gaussian noise signal ($\sigma=10^{-10}\>m$) and examine how the function
$\chi^{2}$ responds by $\Delta\chi^{2}\sim 10^{-20}\>s^{2}$.

If we perturb only the base event ($\sigma=10^{-10}$) on $\chi_{0}$, we get
$\Delta\chi^{2}\sim 10^{-33}\>s^{2}$. If we perturb ($\sigma=10^{-10}$) only
satellite $\chi_{0}$, $\Delta\chi^{2}\sim 10^{-22}\>s^{2}$. If we perturb
($\sigma=10^{-10}$) only satellite $\chi_{1}$, $\Delta\chi^{2}\sim
10^{-22}\>s^{2}$. If we perturb ($\sigma=10^{-10}$) only satellite $\chi_{2}$,
$\Delta\chi^{2}\sim 10^{-24}\>s^{2}$.

If we perturb ($\sigma=10^{-10}$) only one parameter of a satellite, we get
from $\Delta\chi^{2}\sim 10^{-21}\>s^{2}$ to $\Delta\chi^{2}\sim
10^{-37}\>s^{2}$ depending on the perturbation and on the parameter we decide
to perturb and the configuration we start from.

This needs extra care since $\Delta\chi^{2}$ is quite close to the precision
at which we are computing everything. So we went back and recomputed
everything with a precision $10^{-80}$ instead of $10^{-40}$. We were able to
confirm that these figures are really there.

We see that, with a clock sensitive to $10^{-20}\>s$, we can measure
$\chi^{2}$ with a precision near to $10^{-40}$. This corresponds to a change
of parameters of approximately $10^{-20}$.

Indeed, this is the best we can hope for, not necessarily what we can manage
to attain. However, we note that, regarding the positions, we are talking
about $1\>mm$ for a satellite $2\cdot 10^{7}\>m$ away!

Figure 7: The signals of $2^{nd}$ and $3^{rd}$ generations

We summarize by remarking that we have been able to reduce the whole issue of
orbital parameters determination to a minimization. It turns out that we have
several very good algorithms to minimize a function $\chi^{2}$ even when we do
not know much regarding its differentiability. After all, that is what machine
learning is all about.

##  7 Conclusions and perspectives

In summary, we show that in dimension 2+1 we can cope with problems we did not
have in 1+1. We are able to analytically describe the exact bounded orbits to
model the satellites and light rays to model signals exchanged between them.
We are able to solve the control problem in quite a general way, by
determining the signal emitted from a satellite to a target event, possibly on
another satellite. Finally, the simulation phase is robust enough to actually
minimize the function, from which one can determine the orbital parameters.

We would like to provide some preliminary evidence that, in principle, one can
solve the minimizing problem for the function
$\chi^{2}(r_{0},\theta_{i},t_{i},r^{i}_{\pm})$. There are many algorithms
available to do this computation; they mainly differ in efficiency and the
regularity required for the function. In particular, there are several
algorithms which rely only on function evaluations. We remark that we are
using a parameter space which is 11 dimensional. Since we plan to perform the
minimization continuously, we would like to use the previous solution found as
a guess for the solution of the following iteration.

We expect that the performance can be improved by tailoring the algorithm to
the specific situation we are considering. However, at this stage, we are not
concerned with efficiency. We took the standard minimizer of Maple (in the
Optimization package) with a little juggling on its configuration. Then we
model rPS in the Earth situation, we set the mathematical precision to
$10^{-80}$, and the orbital parameters as described above. With this setup we
simulate signals so that we know a solution quite precisely.

We proceed by perturbing the solution ($\sigma=10^{-5}$) to simulate a
previous solution near it. We want to show that the minimizer can start from
this initial perturbation and converge to the known solution. We ask for a
precision of $10^{-20}$ (for the function $\chi^{2}$), which corresponds to a
precision of about $10^{-10}$ for the parameters. The Maple minimizer
returned, as expected, to the solution, which is quite close to the known
solution, namely we observe that each parameter starts from the perturbed
value and ends near the known solution. The results of the calculation are
given in the following table of the absolute accuracy of perturbed values and
the final minimum found with respect to the calculated solution.

$\Delta$perturbed $\Delta$minimum $r_{0}$ $0.57987\cdot 10^{-5}$ $1.0\cdot
10^{-10}$ $r^{0}_{-}$ $1.0037\cdot 10^{-5}$ $1.0\cdot 10^{-10}$ $r^{0}_{+}$
$1.8841\cdot 10^{-5}$ $1.0\cdot 10^{-10}$ $\theta_{1}$ $2.8649\cdot 10^{-7}$
$-0.4096\cdot 10^{-10}$ $t_{1}$ $0.23341\cdot 10^{-5}$ $1.0\cdot 10^{-10}$
$r^{1}_{-}$ $1.662\cdot 10^{-5}$ $1.0\cdot 10^{-10}$ $r^{1}_{+}$ $-1.2022\cdot
10^{-5}$ $1.0\cdot 10^{-10}$ $\theta_{2}$ $-2.3878\cdot 10^{-5}$ $1.0\cdot
10^{-10}$ $t_{2}$ $0.66986\cdot 10^{-5}$ $1.0\cdot 10^{-10}$ $r^{2}_{-}$
$0.93391\cdot 10^{-5}$ $1.0\cdot 10^{-10}$ $r^{2}_{+}$ $0.58673\cdot 10^{-5}$
$1.0\cdot 10^{-10}$ In order to compute the minimization the algorithm has to
do about 100 evaluations of the function $\chi^{2}$, each of which requires
that a complete simulation is made. Currently, each simulation takes about
$60\>s$, namely, about $2\>s$ per signal. By the use of more performance-
oriented programming languages (Swift and C++, using GSL libraries) we
obtained better performance by about a factor of 100, albeit at the price of
having less control on precision and accuracy and of doing the computations
numerically. Here we preferred to go for accuracy and precision, neglecting
performance, since we are discussing viability in principle. Maple’s
simulation can be found at [27].

We stress that even without juggling much with its configuration, the
algorithm converges (though quite slowly) to the solution. The main point here
is that the minimizer requires about 100 evaluations to run. Thus the main
challenge was to obtain a robust simulation procedure, which allows us to
evaluate the function $\chi^{2}$ for all initial conditions in a neighbourhood
of the solution. The results indicate that the function $\chi^{2}$ is well-
defined and continuous almost everywhere near a solution.

We are able to argue that the precision of the clocks, which sets a limit to
our ability to measure $\chi^{2}$, implies a limit in the precision with which
we can determine, in principle, the orbital parameters. In the Earth
configuration, we can reasonably imagine an order of magnitude in the range of
$10^{-3}\>m$.

Under the best conditions, we can envisage an optical clock with a precision
of approximately $10^{-20}\>s$ which allows us the measure $\chi^{2}$ to order
$10^{-40}\>s^{2}$. This precision leads to orbital parameters of order
$10^{-20}$. For the error in position, we have $(\Delta r)^{2}=\Delta
r^{2}+r^{2}\Delta\theta^{2}$ which implies an error of approximately $\Delta
r\sim 10^{-13}\>m$. These results are adequate for a completely theoretical
estimate, since currently this scenario is not yet possible. It can be
considered as an approximate maximal precision estimate.

In this regard we remark that many other sources of error affect the final
result, first of all, atmospheric turbulence. This is why we have designed the
process of determining orbital parameters in order that it is separated from
the determination of user positioning, namely, it happens entirely out of the
atmosphere.

Further investigation needs to be devoted to the domain in which our procedure
is working, as well as for a more realistic evaluation of errors, for which
one needs details about the satellites that are unavailable at the moment.

At this point of our research we have learned that a relativistic theory
provides us with the tools to measure the motion of objects around a source so
that we can bootstrap laws of dynamics without relying on Newtonian physics or
special relativity. The information obtained may be used to measure the
gravitational field (at least in a discretized version of it) by fitting
perturbations of Schwarzschild, and to discuss the observability of modified
gravitational theories in terms of explicitly observable quantities.

We are confident that the whole derivation can be carried out for other
solutions of Einstein equations, e.g., the Kerr spacetime, provided the
solutions possess enough first integrals to determine the geodesic
trajectories in terms of Weierstrass equations, which corresponds to having a
complete integral of the Hamilton-Jacobi equation.

## Appendix A: drawing the light conoid

Using equations of light rays (16), we can fix the vertex event
$p_{0}=(r_{0},\theta_{0},t_{0})$ and parameterize the past light conoid by
$(r,r_{m})$ (or $(r,K^{2})$ depending on the surface region).

As an example, we can consider scattered, ingoing, counterclockwise, light
rays through $p_{0}$ which are parameterized as

$r=r\qquad\qquad\theta=\theta(r,r_{m})\qquad\qquad t=t(r,r_{m})$ (28)

where the functions $\theta(r,r_{m})$ and $t(r,r_{m})$ are given by the
solution in (18) (in this case with $e=-1$ and $\sigma=1$). Accordingly, we
can draw the corresponding region of the light conoid, for clockwise rays, for
counterclockwise and clockwise outgoing rays (see Fig. 8). When we consider
outgoing rays, we can trace back on the surface the curve on which they reach
their minimal $r=r_{m}$ from the BH, and draw, for each outgoing ray, its
ingoing branch, which fills another region on the light conoid. Finally, we
can consider in-falling (incoming and outgoing) rays, which give us two more
regions of the light conoid using the solutions (17). We observe that since we
have analytic solutions, we are drawing the exact analytic light conoid with
no approximations.

It is interesting to draw the past light conoid up to some time $\delta t$
antecedent to its vertex $p_{0}$. In this case, for small enough $\delta t$,
the surface appears to be a conoid, as expected. If we then increase $\delta
t$, multisheets appear. By definition of the light conoid, each of its points
is connected to the vertex by a light ray. Thus, if we consider a particle $P$
with a worldline which intersects the light conoid with vertex $p_{0}$ more
than once, we have a light ray starting at each intersection and reaching
$p_{0}$. This situation shows that the vertex sees multiple images of the
particle $P$ emitted at different times although received at the same time.

a)
b)
c)

Figure 8: The three figures show the light conoid, precisely the past conoid
of the event $p_{0}=(r=5,\theta=\hbox{$\pi\over 4$},t=13)$. The red regions
are the in-falling rays, the blue regions the scattered rays in the branch
containing $p_{0}$, the green parts are the opposite branches of the scattered
rays. Lengths are measured in units $\alpha$, times so that $c=1$. The central
cylinder is the BH horizon. a) shows the interval $t\in[8,15]$. Here the
surface is a cone with a vertex in $p_{0}$ deformed by the BH bending of light
and synchronization. The red line shows a ray through $p_{0}$, which reaches
the minimal spatial distance at the red dot and then goes out on the other
branch. b) shows the interval $t\in[-5,15]$. We see that the generators of the
light conoid self-intersect after going around the BH. c) shows the interval
$t\in[-25,15]$ with two sheet of the light conoid.

## Appendix B: Hamiltonian control problem

The control problem can be extended to a general Hamiltonian system. The
problem consists of the determination of a solution through given initial and
final positions by finding the initial (and final) momenta.

Hamiltonian systems are defined on a symplectic manifold, e.g. the cotangent
bundle $T^{\ast}Q$ of the configuration space $Q$ which is parameterized by
the position coordinates $q$. The dynamics is described by a Hamiltonian
$H(q,p)$ which can be determined by the Legendre transform of the Lagrangian.
In the case of the Lagrangian (3) the momenta are defined by

$p_{r}=-{\dot{r}\over\sqrt{A(c^{2}A^{2}-\dot{r}^{2}-Ar^{2}\dot{\theta}^{2})}}\qquad
p_{\theta}=-{Ar^{2}\dot{\theta}\over\sqrt{A(c^{2}A^{2}-\dot{r}^{2}-Ar^{2}\dot{\theta}^{2})}}$
(29)

and the Hamiltonian by

$H(r,\theta,p_{r},p_{\theta})={c\sqrt{A}\over
r}\sqrt{r^{2}+Ar^{2}p_{r}^{2}+p_{\theta}^{2}}$ (30)

Since the Hamiltonian is independent of $t$ and $\theta$ is cyclic, one can
look for a complete integral of the Hamilton-Jacobi (HJ) equation in the form

$S(t,r,\theta;E,J)=-{c^{2}}Et+J\theta+W(r;E,J)$ (31)

It follows that the HJ equation reduces to

${dW\over dr}=\pm{\sqrt{({c^{2}}E^{2}-A)r^{2}-J^{2}A}\over rA}$ (32)

Thus the solution $W$ (and thence a complete integral $S$) is given by

$W(r;E,J)=\pm\int{\sqrt{({c^{2}}E^{2}-A)r^{2}-J^{2}A}\over rA}dr$ (33)

where the value of the integration constant is immaterial.

The complete integral $S(t,r,\theta;E,J)$ is a generating function for the
canonical flow of transformations

$\Phi_{S}:(t,r,\theta,p_{r},p_{\theta})\mapsto(t,E,J,P_{E},P_{J})$ (34)

where the conjugate momenta $(P_{E},P_{J})$ are also first integrals.

We note that, for any Hamiltonian system, one has a complete integral of the
Hamilton-Jacobi equations. Indeed, the evolution of the system is always a
canonical flow, since the complete integral is nothing but the generating
function of the evolution flow. In fact, one can consider the action
$A=\int_{t_{0}}^{t_{1}}L\circ\sigma dt$ computed along solutions $\sigma$ with
a fixed initial and final positions, write it as a function of initial and
final positions $S(t_{0},q_{0},t_{1},q_{1})$, which is called the Hamilton’s
principal function, and show directly that it satisfies the Hamilton-Jacobi
equation. We note that the distance-function defined in Section 2 is the
square of Hamilton’s principal function [23].

This proves existence. However, to compute Hamilton’s principal function in
this way one needs first to know the general solution of the system. However,
if such a solution is unavailable, we still may be able to find a complete
integral, e.g., by separation of variables, and then show that it satisfies
the HJ equation. As mentioned above the complete integral generates the
evolution flow, hence the general solution of the system. Thus, given a
complete integral of the HJ equation, we can construct the Hamilton’s
principal function, which, in fact, is a kind of classical propagator, since
it solves the control problem of finding solutions that go from
$(t_{0},q_{0})$ to $(t_{1},q_{1})$.

We can compose two such flows

$\Phi_{F}:(t,r,\theta,p_{r},p_{\theta})\mapsto(t,E,J,P_{E},P_{J})\mapsto(t_{0},r_{0},\theta_{0},p^{0}_{r},p^{0}_{\theta})$
(35)

associated with the generating function

$F(t,r,\theta;t_{0},r_{0},\theta_{0})=S(t,r,\theta;E,J)-S(t_{0},r_{0},\theta_{0};E,J)=:S-S_{0}$
(36)

This is called the evolution generator.

Since we know that $(P_{E},P_{J})$ are conserved along solutions one has

$\begin{cases}{\partial F\over\partial E}={\partial S\over\partial
E}-{\partial S_{0}\over\partial E}=P_{E}-P_{E}=0\cr{\partial F\over\partial
J}={\partial S\over\partial J}-{\partial S_{0}\over\partial
J}=P_{J}-P_{J}=0\cr\end{cases}$ (37)

Then, in principle, one can solve this system for
$E=E(t,r,\theta;t_{0},r_{0},\theta_{0})$ and
$J=J(t,r,\theta;t_{0},r_{0},\theta_{0})$, and replace them back into $S-S_{0}$
to obtain $F$ as a function of $(t,r,\theta;t_{0},r_{0},\theta_{0})$ only.

The evolution generator $F(t,r,\theta;t_{0},r_{0},\theta_{0})$ obtained in
this way represents a congruence of solutions (i.e. geodesics) which start at
$t_{0}$ from the point $(r_{0},\theta_{0})$ (i.e. from the event
$(t_{0},r_{0},\theta_{0})$) and go to the event $(t,r,\theta)$ which
corresponds to the first integrals
$\left(E(t,r,\theta;t_{0},r_{0},\theta_{0}),J(t,r,\theta;t_{0},r_{0},\theta_{0})\right)$.
In other words, it solves the typical control problem in which we control the
initial and final position and determine initial and final momenta (or
equivalently, the first integrals $(E,J)$) needed to go from the initial to
the final position (see [22]).

Finally, we remark that if $E$ is negative, we have timelike geodesics, while
if $(E,J)$ goes to infinity and the ratio

$K:={J\over E}={r^{2}\dot{\theta}\over
A}=-{crp_{\theta}\over\sqrt{A(p_{r}^{2}r^{2}A+p_{\theta}^{2}+r^{2})}}$ (38)

remains finite, we have light rays.

Hence, we have a procedure to obtain an explicit description of light rays. We
consider the equations (37), we replace $J=K\epsilon^{-1}$ and
$E=\epsilon^{-1}$ and expand in series of $\epsilon$. Since we want solutions
where $E\rightarrow-\infty$, (i.e. $\epsilon\rightarrow 0$), one necessarily
needs the zero-order term in the series to vanish, otherwise, $\epsilon=0$
would not be a solution. This process yields the equations

$\begin{cases}{\partial F\over\partial
E}(t,r,\theta;t_{0},r_{0},\theta_{0};K,\epsilon=0)=0\cr{\partial
F\over\partial
J}(t,r,\theta;t_{0},r_{0},\theta_{0};K,\epsilon=0)=0\cr\end{cases}$ (39)

the solutions of which: $t=t(r;t_{0},r_{0},\theta_{0},K)$,
$\theta=\theta(r;t_{0},r_{0},\theta_{0},K)$, determine a light ray through
$(t_{0},r_{0},\theta_{0})$. Different values of $K$ parameterize different
“types” of light rays (since the value of $K$ is constant along a single light
ray).

Thus, given two events $(t_{0},r_{0},\theta_{0})$ and $(t,r,\theta)$, it is
easy to verify whether there is a light ray joining them: one checks whether
there is a value of $K$ such that equations (39) are satisfied, i.e., if they
have a solution as equations for $K$.

The typical problem we have to solve for computing signals is given a fixed
event $P=(t_{0},r_{0},\theta_{0})$ and a timelike curve
$\gamma:\mathbb{R}\rightarrow M:s\mapsto(t(s),r(s),\theta(s))$, find a point
on the trajectory of $\gamma$ such that it is connected to $P$ by a light ray.
Then we simply write down equations (39) as

$\begin{cases}{\partial F\over\partial
E}(t(s),r(s),\theta(s);t_{0},r_{0},\theta_{0};K,\epsilon=0)=0\cr{\partial
F\over\partial
J}(t(s),r(s),\theta(s);t_{0},r_{0},\theta_{0};K,\epsilon=0)=0\cr\end{cases}$
(40)

and solve them for $(s,K)$. The value of $s$ gives the event
$Q=(t(s),r(s),\theta(s))$ on the trajectory of $\gamma$, while the value of
$K$ selects the light ray through $P$ and $Q$.

## Acknowledgements

This article is based upon work from COST Action (CA15117 CANTATA), supported
by COST (European Cooperation in Science and Technology).

We also acknowledge the contribution of INFN (Iniziativa Specifica QGSKY and
Iniziativa Specifica Euclid), the local research project Metodi Geometrici in
Fisica Matematica e Applicazioni (2022) of Dipartimento di Matematica of
University of Torino (Italy). This paper is also supported by INdAM-GNFM.

L. Fatibene would like to acknowledge the hospitality of the Department of
Applied Mathematics, University of Waterloo, where part of this research was
done. This work was supported in part by a Discovery Grant from the Natural
Sciences and Engineering Research Council of Canada (R.G. McLenaghan).

## References

  * [1] S.Carloni, L.Fatibene, M.Ferraris, R.G. McLenaghan, P.Pinto, Discrete Relativistic Positioning Systems, GRG, 52, (2), (2020)
  * [2] B. Coll, in: Proc. ERE-2000 Meeting on Reference Frames and Gravitomagnetism, eds. J. F. Pascual- Snchez, L. Flora, A. San Miguel, and F. Vicente (World Scientific, Singapore, 2001) 53.
  * [3] M. Blago jevc, J. Garecki, F. W. Hehl, and Y. N. Obukhov, Real null coframes in general relativity and GPS type coordinates, Phys. Rev. D65 (2002) 044018.
  * [4] B.Coll, Relativistic Positioning Systems, in: Proc. Spanish Relativity Meeting ERE-2005, Oviedo (Spain); arXiv:gr-qc/0601110
  * [5] B.Coll, J.J.Ferrando, J.A.Morales, Two-dimensional approach to relativistic positioning systems, Phys. Rev. D73 (2006) 084017; arXiv:gr-qc/0602015
  * [6] B.Coll, J.M.Pozo, Relativistic Positioning Systems: The Emission Coordinates, Class. Quantum Grav.23 (2006) 7395; arXiv:gr-qc/0606044
  * [7] B.Coll, J.J.Ferrando, J.A.Morales, Positioning with stationary emitters in a two-dimensional space-time, Phys. Rev. D74 (2006) 104003; arXiv:gr-qc/0607037
  * [8] A. Tarantola, L. Klimes, J. M. Pozo and B. Coll, Gravimetry, Relativity, and the Global Navigation Satellite Systems, arXiv:gr-qc/0905.3798
  * [9] M. Lachieze-Rey, The Covariance of GPS Coordinates and Frames, Classical and Quantum Gravity 23.10 (2006), p. 3531. doi: 10.1088/0264- 9381/23/10/019.
  * [10] C.F. Gauss, Untersuchungen über Gegenstände der Höheren Geodäsie. Erste Abhandlung, Abhandlungen der Königlichen Gesellschaft der Wissenschaften in Göttingen. Zweiter Band: 3–46.
  * [11] D. Bini et al., Emission versus Fermi Coordinates: Applications to Relativistic Positioning Systems. Classical and Quantum Gravity 25.20 (2008), p. 205011. doi: 10.1088/0264-9381/25/20/205011.
  * [12] C. Rovelli, GPS observables in general relativity, Phys. Rev. D65 (2002) 044017; arXiv:gr- qc/0110003
  * [13] P.Schneider, J.Ehlers, E.E.Falco Graviational Lenses, Springer-Verlag (Berlin, 1975) 514.
  * [14] J. Hadamard Lectures on Cauchy’s Problem in Linear Partial Differential Equations, Yale University Press (1923) 89
  * [15] H.S. Ruse, Taylor?s theorem in the tensor calculus, Proc. London Math. Soc.32 (1931) 87
  * [16] H.S.Ruse, An absolute partial differential calculus, Quart. J. Math. Oxford Ser. 2 (1931) 190
  * [17] H.S.Ruse, A.G.Walker, T.J.Willmore. Harmonic Spaces, Ed. Cremonese Roma, (1961) 13
  * [18] J.L.Synge, A characteristic function in Riemannian space and its applications to the solution of geodesic triangles, Proc. London Math. Soc.32 (1931) 241.
  * [19] K.Yano, Y.Muto, Notes on the derivation of geodesics and the fundamental scalar in a Riemannian space, Proc. Phys.-Math. Soc. Jap. 18 (1936) 142.
  * [20] J.A.Schouten, Ricci-Calculus. An introduction to tensor analysis and its geometrical applications, Springer-Verlag (1954) 382.
  * [21] F.G.Friedlander, The wave equation on a curved space-time, Cambridge University Press (Cambridge, 1975) 17.
  * [22] S.Benenti, Hamiltonian Structures and Generating Families, Universitext (Springer-Verlag New York, 2011).
  * [23] J.L.Synge, Relativity: The General Theory, North-Holland (Amsterdam, 1960).
  * [24] J.Ehlers, F.A.E.Pirani, A.Schild, The Geometry of Free Fall and Light Propagation, in: General Relativity, ed. L.O.‘Raifeartaigh (Clarendon, Oxford, 1972).
  * [25] M. Di Mauro, L. Fatibene, M.Ferraris, M.Francaviglia, Further Extended Theories of Gravitation: Part I, Int. J. Geom. Methods Mod. Phys. Volume: 7, Issue: 5 (2010), pp. 887-898; gr-qc/0911.2841
  * [26] L. Fatibene, M. Ferraris, M. Francaviglia, S. Mercadante, Further Extended Theories of Gravitation: Part II, Int. J. Geom. Meth. Mod. Phys. 7, (2010); arXiv:0911.2842 [gr-qc]
  * [27] Maple simulations can be found at http://www.fatibene.org/gps.html