ges.tex

\chapter{Geodesics in Schwarzschild spacetime: generic and timelike cases}
\label{s:ges}
\index{geodesic!in Schwarzschild spacetime}

\minitoc

\section{Introduction}

We have already investigated some geodesics of Schwarzschild spacetime in
Chap.~\ref{s:sch}, namely
the radial null geodesics (Secs.~\ref{s:sch:rad_null_geod} and \ref{s:sch:radial_null_IEF}).
Here, we perform an extensive study.
After having established the main properties of generic causal (timelike or null) geodesics
in Sec.~\ref{s:ges:geod_motion}, we investigate timelike
geodesics in Sec.~\ref{s:ges:timelike}. They
are of great physical importance, since they represent
orbits of planets or stars around the black hole, as well as worldlines of
intrepid observers freely falling into the black hole.
The study of null geodesics, which govern images received by observers,
is deferred to Chap.~\ref{s:gis}.

Before studying this chapter, the reader might want to have a look at
Appendix~\ref{s:geo}, which recaps the properties of geodesics in
manifolds equipped with a metric.
It could also be worth to read again Sec.~\ref{s:fra:worldlines} about
worldlines of particles.

\section{Geodesic motion} \label{s:ges:geod_motion}

Let $\Li$ be a geodesic\footnote{The definition and basic properties of geodesics
are recalled in Appendix~\ref{s:geo}; see also Sec.~\ref{s:fra:geod_motion}.} of Schwarzschild spacetime
$(\M,\w{g})$. We shall assume that $\Li$ is causal, i.e. either timelike or null\footnote{As
shown in Sec.~\ref{s:geo:def}, a geodesic cannot be partly timelike and partly
null.}. It therefore can be considered as the worldline
of some particle $\mathscr{P}$, either massive
($\Li$ timelike) or massless ($\Li$ null). As recalled in Sec.~\ref{s:fra:worldlines},
the worldline of a particle $\mathscr{P}$ is a geodesic if, and only if,
$\mathscr{P}$ is submitted only to gravitation, i.e.
$\mathscr{P}$ is in free fall\index{free fall}.


\subsection{First integrals of motion} \label{s:ges:fiom}

The Schwarzschild spacetime $(\M,\w{g})$ is static and spherically symmetric; the
Killing vector $\w{\xi}$ associated with the staticity (cf. Sec.~\ref{s:sch:static_spher})
and the Killing vector $\w{\eta}$ associated with the rotation symmetry along some
axis, give birth to two conserved quantities along $\Li$:
\begin{prop}[conserved quantities along causal geodesics]
Denoting by $\w{p}$ the 4-momentum of particle
$\mathscr{P}$ that follows the geodesic $\Li$
(cf. Sec.~\ref{s:fra:worldlines}), the scalar products
\begin{subequations}
\label{e:ges:conserved_quantities}
\begin{align}
& \encadre{E := - \w{\xi}\cdot \w{p} = - \w{g}(\w{\xi},\w{p}) } \label{e:ges:conserved_energy} \\
& \encadre{L := \w{\eta}\cdot \w{p} = \w{g}(\w{\eta},\w{p}) } , \label{e:ges:conserved_angu_mom}
\end{align}
\end{subequations}
are constant along $\Li$.
The scalar $E$ is called $\mathscr{P}$'s
\defin{conserved energy}\index{conserved!energy}\index{energy!conserved --}
or \defin{energy at infinity}\index{energy!at infinity},
while $L$ is called $\mathscr{P}$'s \defin{conserved angular momentum}\index{conserved!angular momentum}\index{angular momentum!conserved --}
or \defin{angular momentum at infinity}\index{angular momentum!at infinity}.
\end{prop}
\begin{proof}
The 4-momentum $\w{p}$ is a tangent vector associated with an affine parameter
of $\Li$, i.e. it obeys the geodesic equation (\ref{e:fra:p_geodesic}).
The constancy of $E$ and $L$ follow then from
Property~\ref{p:geo:g_xi_v_const} (Sec.~\ref{s:geo:sym}).
\end{proof}
In coordinates $(t,r,\th,\ph)$ adapted to the spacetime symmetries,
i.e. coordinates such that $\w{\xi} = \wpar_t$ and $\w{\eta}=\wpar_\ph$,
for instance the Schwarzschild-Droste
coordinates or the Eddington-Finkelstein ones, one can rewrite
(\ref{e:ges:conserved_quantities})
in terms of the components $p_t = g_{t\mu} \, p^\mu$ and $p_\ph = g_{\ph\mu} \, p^\mu$
of the 1-form $\uu{p}$ associated to $\w{p}$ by metric duality:
\begin{subequations}
\begin{align}
& E = - p_t \\
& L = p_\ph
\end{align}
\end{subequations}
Indeed, in such a coordinate system, $\xi^\mu =  \delta^\mu_{\ \, t}$
and $\eta^\mu = \delta^\mu_{\ \, \ph}$, so that $E = -g_{\mu\nu} \, \xi^\mu p^\nu = -g_{t\nu} \, p^\nu = -p_t$
and $L = g_{\mu\nu} \, \eta^\mu p^\nu = g_{\ph\nu} \, p^\nu = p_\ph$.


It is worth stressing that $E$ is not a genuine energy, i.e. it is not
an energy measured by some observer. Indeed, the latter is defined by
Eq.~(\ref{e:fra:E_obs}), which resembles Eq.~(\ref{e:ges:conserved_energy})
but differs from it by $\w{\xi}$ not being a unit vector in general:
$\w{\xi}\cdot\w{\xi} \not = -1$. In other words, $\w{\xi}$ cannot
be interpreted as the 4-velocity of some observer, so that the quantity
$E$ defined by (\ref{e:ges:conserved_energy}) cannot be a \emph{physically measured}
particle energy. It is only in the asymptotic region, where $\w{\xi}\cdot\w{\xi} = g_{tt}
\to -1$, that $\w{\xi}$ is eligible as a 4-velocity, hence the name
\emph{energy at infinity}. Note that this name is commonly used, even in the
particle $\mathscr{P}$ never visits the asymptotic region.
Similarly, $L$ is not some (component of a) genuine angular momentum. Only in the
asymptotic region do we have
\be \label{e:ges:ell_asympt}
    L \simeq g_{\ph\ph} \, p^\ph \simeq r^2\sin^2\th \, p^\ph \simeq r^2\sin^2\theta \, P^\ph
    \simeq r\sin\th \, P^{(\ph)} ,
\ee
where $P^{(\ph)}$ is the azimuthal component of the momentum $\w{P}$ of particle $\mathscr{P}$
as measured by an asymptotic inertial observer  $\Obs$ (cf. Sec.~\ref{s:fra:measure}), i.e.
the component of $\w{P}$ along $\w{e}_{(\ph)}$ in the orthonormal basis $(\w{e}_{(r)}, \w{e}_{(\th)}, \w{e}_{(\ph)})$, with $\w{e}_{(\ph)} = (r\sin\th)^{-1} \wpar_\ph$.
In view of (\ref{e:ges:ell_asympt}), we may say that $L$ is the angular momentum
about the symmetry axis $\th=0$ that $\Obs$ would attribute to
particle $\mathscr{P}$ if the latter would move close to him. Equivalently,
in a Cartesian coordinate system defined by $(x,y,z) = (r\sin\th\cos\ph, r\sin\th\sin\ph, r\cos\th)$,
$L$ is the component $L_{\rm tot}^z$ of the total angular momentum
of $\mathscr{P}$ as measured by $\Obs$:
\be \label{e:ges:L_tot_asympt}
    \w{L}_{\rm tot} := \w{r}\times \w{P} = L_{\rm tot}^x \, \wpar_x + L_{\rm tot}^y \, \wpar_y
        + L \, \wpar_z .
\ee

\begin{prop}[positivity of the conserved energy]
\label{p:ges:E_positive_M_I}
The conserved energy $E$ is a positive
quantity as soon as the geodesic $\Li$ has some part in $\M_{\rm I}$, i.e.
some part with $r>2m$:
\be \label{e:ges:E_positive_M_I}
    \Li \cap \M_{\rm I} \neq \varnothing \quad \Longrightarrow \quad E > 0 .
\ee
\end{prop}
\begin{proof}
In $\M_{\rm I}$, the Killing vector $\w{\xi}$ is timelike and future-directed.
The 4-momentum $\w{p}$ is either timelike or null and always future-directed.
By Eq.~(\ref{e:fra:orient_time1}) in Lemma~\ref{p:fra:lem1} (Sec.~\ref{s:fra:time_orientation}), one has then necessarily $\w{\xi}\cdot\w{p} < 0$; hence Eq.~(\ref{e:ges:conserved_energy})
implies $E > 0$ in $\M_{\rm I}$. Since $E$ is constant along $\Li$, it
follows that $E > 0$ everywhere.
\end{proof}
\begin{remark}
If the geodesic $\Li$ is confined to $\M_{\rm II}$, i.e. to the black hole
region (cf. Sec.~\ref{s:sch:BH}),
where $\w{\xi}$ is spacelike (cf. Sec.~\ref{s:sch:SD_domain}),
it is possible to have $E \leq 0$, since the
scalar product of $\w{p}$ with a spacelike vector can take any value.
\end{remark}
\begin{remark}\label{r:ges:L_any_sign}
The Killing vector $\w{\eta}$ being always spacelike,
the scalar product $\w{g}(\w{\eta},\w{p})$ can a priori take any real value, and
thus there is no constraint
on the sign of $L$.
\end{remark}

To be specific, let us describe Schwarzschild spacetime in terms of the
Schwarzschild-Droste coordinates $(t,r,\th,\ph)$ introduced in Sec.~\ref{s:sch:solving_EE}.
Without any loss of generality, we may choose these coordinates so that
at $t=0$, the particle $\mathscr{P}$ is located in the equatorial plane $\th=\pi/2$ and
the spatial projection of the worldline $\Li$ lies in that plane, i.e. $\w{p}$ has
no component along $\wpar_\th$:
\be
    \w{p} \stackrel{t=0}{=} p^t \wpar_{t} + p^r \wpar_r + p^\ph \wpar_\ph .
\ee
Now, for $t>0$, if the geodesic $\Li$ were departing from $\th=\pi/2$, this
would constitute some breaking of spherical symmetry, making a difference
between the ``Northern'' hemisphere and the ``Southern'' one.
Hence\footnote{More rigorously, Eq.~(\ref{e:ges:pth_zero}) can
be derived from the geodesic equation (\ref{e:fra:p_geodesic}): given
the expression of the Christoffel symbols of $\w{g}$ in Schwarzschild-Droste
coordinates (cf. Sec.~\ref{s:sam:Kottler_solution}), Eq.~(\ref{e:fra:p_geodesic})
yields \[\frac{\D p^\th}{\D\lambda} +\frac{2}{r}  p^r p^\th - \sin\th\cos\th \left( p^\ph \right)^2 = 0,\] where $\lambda$ is the affine parameter of $\Li$ associated with $\w{p}$,
so that $p^\th=\D\th/\D\lambda$. Whatever the values of $r(\lambda)$,
$p^r(\lambda)$ and $p^\ph(\lambda)$,
the solution of this ordinary differential equation with the initial conditions $p^\th=0$ and
$\cos\th=0$ is $p^\th=0$ for all values of $\lambda$.} $\Li$
must stay at $\th=\pi/2$, which implies
\be \label{e:ges:pth_zero}
    \encadre{p^\th = 0} .
\ee
We conclude:
\begin{prop}[planar character of geodesics]
A geodesic $\Li$ of Schwarzschild spacetime is necessarily confined to a timelike hypersurface.
Without any loss of generality, we can choose Schwarzschild-Droste coordinates $(t,r,\th,\ph)$
such that this hypersurface is the ``equatorial hyperplane'' $\th=\pi/2$.
Then the component $p^\th$ of the 4-momentum of the particle having $\Li$ as worldline
vanishes identically [Eq.~(\ref{e:ges:pth_zero})].
\end{prop}

Let us denote by $\mu$ the mass of particle $\mathscr{P}$, with possibly
$\mu=0$ if $\mathscr{P}$ is a photon. The scalar square of the 4-momentum $\w{p}$ is
then [cf. Eq.~(\ref{e:fra:def_mass})]
\be \label{e:ges:p2_mu2}
    \w{g}(\w{p},\w{p}) = - \mu^2 .
\ee
So $\mu^2$ can be seen as an integral of motion.

\subsection{Equations of motion and generic properties} \label{s:ges:eq_to_be_solved}

Contemplating Eqs.~(\ref{e:ges:conserved_energy}), (\ref{e:ges:conserved_angu_mom}),
(\ref{e:ges:pth_zero}) and (\ref{e:ges:p2_mu2}), we realize that we have
four first integral of motions. The problem is then completely integrable.
More specifically, let $\lambda$ be the affine parameter along the geodesic
$\Li$ associated with the 4-momentum $\w{p}$ [cf. Eq.~(\ref{e:geo:v_dxdlambda})]:
\be \label{e:ges:def_lambda}
   \encadre{ \w{p} = \frac{\D\w{x}}{\D\lambda} },
\ee
where  $\D\w{x}$ is the infinitesimal displacement along $\Li$ corresponding
to the parameter change $\D\lambda$. Note that $\lambda$ is dimensionless
and necessarily increases towards the future\footnote{Let us recall
that Schwarzschild spacetime is time-oriented, cf. Sec.~\ref{s:sch:time_orientation}.},
since $\w{p}$ is by definition future-oriented (cf. Sec.~\ref{s:fra:worldlines_def}).
In terms of the components with respect to Schwarzschild-Droste coordinates,
this yields
\be \label{e:ges:comp_4_momentum}
    \dot{t} := \frac{\D t}{\D\lambda} = p^t,\qquad
    \dot{r} := \frac{\D r}{\D\lambda} = p^r,\qquad
    \dot{\th} := \frac{\D \th}{\D\lambda} = p^\th,\qquad
    \dot{\ph} := \frac{\D \ph}{\D\lambda} = p^\ph .
\ee
In the present case, where $\th(\lambda)=\pi/2$, we have of course $\dot{\th}=0$,
in agreement with Eq.~(\ref{e:ges:pth_zero}).
Given the components (\ref{e:sch:Schwarz_metric_SD}) of Schwarzschild metric
with respect to the Schwarzschild-Droste coordinates,
Eq.~(\ref{e:ges:conserved_energy}) can be written as
\[
    E = - g_{t\mu} p^\mu = - g_{tt} p^t = - g_{tt} \dot{t}
    = \left(1 - \frac{2m}{r} \right) \dot{t} ,
\]
hence
\be \label{e:ges:dot_t}
    \encadre{\frac{\D t}{\D \lambda} = E \left(1 - \frac{2m}{r} \right) ^{-1} }.
\ee
Similarly, Eq.~(\ref{e:ges:conserved_angu_mom}) becomes
\[
    L = g_{\ph\mu} p^\mu = g_{\ph\ph} p^\ph  = g_{\ph\ph} \dot{\ph}
        = r^2 \sin^2\th \, \dot{\ph} .
\]
Since $\th=\pi/2$, we get
\be \label{e:ges:dot_ph}
    \encadre{\frac{\D\ph}{\D\lambda} = \frac{L}{r^2} } .
\ee
We have already noticed that the sign of $L$ is unconstrained (Remark~\ref{r:ges:L_any_sign} on p.~\pageref{r:ges:L_any_sign}). The above equation
shows that it corresponds to the increase ($L>0$) or decrease $(L<0)$
of $\ph$ along the geodesic $\Li$. In other words, we deduce from Eq.~(\ref{e:ges:dot_ph})
that
\begin{prop}[monotonic behaviour of $\ph$ along geodesics]
Along any timelike or null geodesic of Schwarzschild spacetime,
the azimuthal coordinate $\ph$ is either constant ($L=0$)
or increases (resp. decreases) monotonically ($L>0$) (resp. $L<0$).
\end{prop}

The last unexploited first integral of motion is Eq.~(\ref{e:ges:p2_mu2}); it
yields
\[
   - \left( 1 - \frac{2m}{r} \right) (\dot{t})^2 +
   \left( 1 - \frac{2m}{r} \right) ^{-1}  (\dot{r})^2
   + r^2 (\dot{\th})^2 + r^2 \sin^2\th (\dot{\ph})^2  = - \mu^2 .
\]
Using (\ref{e:ges:dot_t}), (\ref{e:ges:dot_ph}), as well as
$\dot{\th}=0$ and $\th=\pi/2$, we get
\[
    -  E^2 \left(1 - \frac{2m}{r} \right) ^{-1}
    +  \left( 1 - \frac{2m}{r} \right) ^{-1}  (\dot{r})^2
    +  \frac{L^2}{r^2} = - \mu^2 ,
\]
which can be recast as
\be \label{e:ges:dot_r_square}
    \encadre{ \left( \frac{\D r}{\D \lambda} \right) ^2
        - \frac{2 \mu^2 m}{r} + \frac{L^2}{r^2}
         \left( 1 - \frac{2m}{r} \right) = E^2 - \mu^2 } .
\ee
To summarize, the geodesic motion in Schwarzschild spacetime is governed by
Eqs.~(\ref{e:ges:dot_t}), (\ref{e:ges:dot_ph}) and (\ref{e:ges:dot_r_square}),
where $r=r(\lambda)$ and $\mu$, $E$ and $L$ are constants.
This constitutes a system of 3 differential equations for the 3 unknown
functions $t(\lambda)$, $r(\lambda)$ and $\ph(\lambda)$.
We observe
that Eq.~(\ref{e:ges:dot_r_square}) is decoupled from the other two equations.
The task is then to first solve this equation for $r(\lambda)$ and to inject
the solution into Eqs.~(\ref{e:ges:dot_t}) and (\ref{e:ges:dot_ph}), which
can then be integrated separately.

A constraint to keep in mind is that the 4-momentum vector $\w{p}$, whose
components are related to the solution $(t(\lambda), r(\lambda), \ph(\lambda))$
by Eq.~(\ref{e:ges:comp_4_momentum}), has to be a future-directed causal vector.
In $\M_{\rm I}$, as we have seen above, this is guaranteed by choosing $E > 0$
[cf. Eq.~(\ref{e:ges:E_positive_M_I})]. In $\M_{\rm II}$, a future-directed
timelike vector is $-\wpar_r$ (cf. Sec.~\ref{s:sch:time_orientation}).
According to Eq.~(\ref{e:fra:orient_time1}) in Lemma~\ref{p:fra:lem1} (Sec.~\ref{s:fra:time_orientation}),
we have then
$\w{p}$ future-directed iff $-\wpar_r\cdot\w{p} < 0$, i.e. iff
\[
    \left(\frac{2m}{r} - 1  \right) ^{-1} p^r < 0  .
\]
Since $2m/r - 1 > 0$ in $\M_{\rm II}$, this is equivalent to $p^r < 0$, i.e.
to $\D r/\D\lambda < 0$. Hence
\begin{prop}[decreasing of $r$ in the black hole region]
In the black hole region $\M_{\rm II}$, i.e. for $r<2m$,
the solution $r(\lambda)$ of Eq.~(\ref{e:ges:dot_r_square})
must be a strictly decreasing function of $\lambda$.
\end{prop}
Actually, we recover Property~\ref{p:sch:r_decreasing}, which has been
derived for any causal worldline (not necessarily
a geodesic) in Sec.~\ref{s:sch:time_orientation}.

\begin{remark}
We have derived the system of Eqs.~(\ref{e:ges:dot_t}), (\ref{e:ges:dot_ph}) and (\ref{e:ges:dot_r_square}) without invoking explicitly the
famous \emph{geodesic equation}\index{geodesic!equation}, i.e.
Eq.~(\ref{e:geo:eq_geod}) in Appendix~\ref{s:geo}.
This is because we had enough
first integrals of the second-order differential equation
(\ref{e:geo:eq_geod})
to completely reduce it to
a system of first order equations.
\end{remark}

\subsection{Trajectories in the orbital plane} \label{s:ges:trajectories}

If the conserved angular momentum vanishes, $L=0$, the equation of motion
(\ref{e:ges:dot_ph}) implies that $\ph = \mathrm{const} = \ph_0$. The geodesic
$\Li$ is then confined to the 2-dimensional timelike surface
$(\th, \ph) = (\pi/2,\ph_0)$, which is spanned by the coordinates $(t, r)$.
One says that $\Li$ is a \defin{radial geodesic}\index{radial!geodesic}.

In the remainder of this section, we discuss the opposite case, namely we
assume
\be
    L \not = 0 .
\ee
We have stressed above that $\ph$ is then a strictly monotonic function of $\lambda$,
increasing (resp. decreasing) continuously along $\Li$ for $L>0$ (resp. $L<0$).
Consequently,
\begin{prop}[$\ph$ as a parameter along geodesics]
Along any timelike or null geodesic with $L\not =0$,
$\ph$ can be chosen as a parameter, provided one does not restrict its
range to $(0,2\pi)$.
\end{prop}
Contrary to $\lambda$, $\ph$ is not in general an affine parameter of $\Li$.
Indeed, the dependency $r = r(\lambda)$ in Eq.~(\ref{e:ges:dot_ph})
does not correspond to an affine relation between $\ph$ and $\lambda$, except for
$r(\lambda) = \mathrm{const}$ (case of circular orbits).

Let us use Eq.~(\ref{e:ges:dot_ph}) to write
$\D r/ \D\lambda = \D r/\D\ph \times \D\ph/\D\lambda = L/r^2 \D r /\D\ph$
and substitute this expression into the equation of motion
(\ref{e:ges:dot_r_square}). We get
\be
    \frac{1}{r^4} \left( \frac{\D r}{\D\ph} \right)^2 =
     \frac{2m}{r^3} - \frac{1}{r^2} + 2  \left( \frac{\mu}{L} \right)^2 \frac{m}{r}
     + \left( \frac{E}{L} \right)^2 - \left( \frac{\mu}{L} \right)^2 .
\ee
To simplify this equation, it is natural to introduce the dimensionless variable
\be \label{e:ges:def_u}
    \encadre{ u := \frac{m}{r} }
\ee
instead of $r$. We then get
\be \label{e:ges:DuDph_trajectories}
    \encadre{
    \left( \frac{\D u}{\D\ph} \right)^2 = 2 u^3 - u^2 + 2 \left( \frac{m\mu}{L} \right)^2 u
    + \left( \frac{m E}{L} \right)^2 - \left( \frac{m \mu}{L} \right)^2   }.
\ee
This differential equation determines entirely the $(r,\ph)$-part of the geodesic $\Li$,
which we shall call the
\defin{trajectory of $\Li$ in the orbital plane}\index{trajectory!of a geodesic}.
The term ``orbital plane'' is a slight abuse of language for
the 2-dimensional \emph{surface} $(t,\th) = (t_0, \pi/2)$, where $t_0$ is a constant.

In general Eq.~(\ref{e:ges:DuDph_trajectories}) is not solvable in terms of
elementary functions. The exceptions are circular orbits ($u = \mathrm{const}$, the
constant being one of the roots of the cubic polynomial in $u$ in the right-hand side)
and the critical null geodesics, which we shall discuss in
Sec.~\ref{s:gis:crit_geod}.
For the generic case, exact solutions are expressible in terms of
some non-elementary special functions; there
are basically two strategies:
\begin{itemize}
\item The first one is to invoke the
\emph{Weierstrass elliptic function}\index{Weierstrass elliptic function}\index{elliptic!Weierstrass -- function}\footnote{The character $\wp$
is a kind of calligraphic lowercase p, which is standard to denote this function.}
$\wp(z; \omega_1,\omega_2)$, which is a doubly-periodic meromorphic function of
the complex variable $z$, of periods $\omega_1\in\mathbb{C}$ and $\omega_2\in\mathbb{C}$.
Among the many properties of this function, the one relevant here is that
$\wp$ is a solution of the differential equation
\be \label{e:ges:Weierstrass_eq}
    \left( \wp'(z) \right)^2  = 4 \wp(z)^3 - g_2 \wp(z) - g_3,
\ee
where $g_2$ and $g_3$ are two constants entirely determined by the
periods $\omega_1$ and $\omega_2$ of $\wp$. Indeed, via
the change of variables $v := u - 1/6$ and $\tilde\ph := \ph / \sqrt{2}$,
Eq.~(\ref{e:ges:DuDph_trajectories}) is equivalent to
\be \label{e:ges:dvdtph_cubic}
    \left( \frac{\D v}{\D\tilde\ph} \right)^2 =
        4 v^3 - \left[\frac{1}{3} - 4 \left( \frac{m\mu}{L} \right)^2 \right] v
        + 2 \left(\frac{m E}{L} \right) ^2
        - \frac{4}{3} \left( \frac{m\mu}{L} \right) ^2 - \frac{1}{27},
\ee
which is obviously of type (\ref{e:ges:Weierstrass_eq}) (no square in the cubic polynomial).
The solution is thus
\be
    u = \frac{m}{r} = \wp\left( \frac{\ph}{\sqrt{2}} + C; \omega_1, \omega_2 \right) + \frac{1}{6} ,
\ee
where $C\in \mathbb{C}$ is a constant and $\omega_1$ and $\omega_2$ are
determined by $m$, $\mu$, $E$ and $L$ (see e.g. Ref.~\cite{GibboV12} for more details
regarding this method applied to null geodesics).
\item The second approach consists in noticing that the method of separation
of variables can easily be applied to Eq.~(\ref{e:ges:DuDph_trajectories}),
leading to
\be
    \ph = \pm \int_{u_0}^u \frac{\D\bar{u}}{\sqrt{
     2 \bar{u}^3 - \bar{u}^2 + 2 \left( \frac{m\mu}{L} \right)^2 \bar{u}
    + \left( \frac{m E}{L} \right)^2 - \left( \frac{m \mu}{L} \right)^2  } }
    + \ph_0 ,
\ee
where $u_0$ and $\ph_0$ are two constants, and the $\pm$ sign can be $+$
on some parts of the geodesic $\Li$ and $-$ on some other parts of $\Li$.
The integral in right-hand side is expressible in terms of the so-called
\emph{incomplete elliptic integrals
of the first kind}\index{incomplete!elliptic integral}\index{elliptic!integral!incomplete --}.
We shall detail such a technique for null geodesics in Sec.~\ref{s:gis:planar_trajectories}.
Note that this approach leads to $\ph = \ph(r)$, whereas the method involving the
Weierstrass function leads to the ``polar equation'' form: $r = r(\ph)$. It is possible though to
get the polar form by invoking the inverses of elliptic integrals, namely
\emph{Jacobi elliptic functions}\index{Jacobi!elliptic function}\index{elliptic!Jacobi -- function}.
\end{itemize}

Once the solution $r=r(\ph)$ of Eq.~(\ref{e:ges:DuDph_trajectories})
has been obtained, it can be injected into the equation for $t = t(\ph)$
that can be deduced from the equations of motions (\ref{e:ges:dot_t})-(\ref{e:ges:dot_ph}):
\be
    \frac{\D t}{\D\ph} = \frac{E}{L} \frac{r(\ph)^3}{r(\ph) - 2 m} .
\ee
This is an ordinary differential equation for $t=t(\ph)$, the solution of
which amounts to finding a primitive with respect to $\ph$ of the
right-hand side. Unfortunately, this is not an easy task in general, the function
$r(\ph)$ being quite involved, except for circular orbits ($r(\ph)=\mathrm{const}$).

In what follows, we discuss separately the resolution of
the system (\ref{e:ges:dot_t})-(\ref{e:ges:dot_r_square}) or of
Eq.~(\ref{e:ges:DuDph_trajectories})
for timelike geodesics
(Sec.~\ref{s:ges:timelike}) and for null geodesics (Chap.~\ref{s:gis}).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Timelike geodesics} \label{s:ges:timelike}

\subsection{Effective potential} \label{s:ges:eff_pot_timelike}

When the geodesic $\Li$ is timelike, it is natural to use the proper time $\tau$
as an affine parameter along it, instead of the parameter $\lambda$ associated
with the 4-momentum $\w{p}$. Since the tangent vector associated with $\tau$
is the 4-velocity $\w{u}$ (cf. Sec.~\ref{s:fra:massive_part}) and $\w{p}$ and $\w{u}$ are related by
Eq.~(\ref{e:fra:p_m_u}): $\w{p} = \mu \, \w{u}$, we get
$\D\w{x}/\D\lambda = \mu \, \D\w{x}/\D\tau$, from which we infer the relation
between $\tau$ and $\lambda$:
\be
    \tau = \mu \lambda ,
\ee
up to some additive constant. This is
of course a special case of the generic relation (\ref{e:geo:affine_transf})
between two affine parameters of the same geodesic.
Equation~(\ref{e:ges:dot_r_square}) becomes then
\be \label{e:ges:1d_motion_timelike}
    \encadre{ \frac{1}{2} \left( \frac{\D r}{\D \tau} \right) ^2
        + V_{\ell}(r) = \frac{\veps^2 - 1}{2} } ,
\ee
where
\be \label{e:ges:V_eff_timelike}
     \encadre{V_{\ell}(r) := - \frac{m}{r} + \frac{\ell^2}{2 r^2}
     \left( 1 - \frac{2 m}{r} \right)}
\ee
and $\veps$ and $\ell$ are respectively the
\defin{specific conserved energy}\index{specific!conserved!energy}
and \defin{specific conserved angular momentum}\index{specific!conserved!angular momentum}
of particle $\mathscr{P}$:
\be \label{e:ges:def_eps_ell}
  \encadre{ \veps := \frac{E}{\mu} = - \w{\xi}\cdot \w{u}}\qquad\mbox{and}\qquad
  \encadre{ \ell :=  \frac{L}{\mu} = \w{\eta}\cdot\w{u}} ,
\ee
where $\w{u}$ is the 4-velocity of $\mathscr{P}$ and the second equalities
result from definitions (\ref{e:ges:conserved_quantities})
and the relation $\w{p} = \mu \, \w{u}$ [Eq.~(\ref{e:fra:p_m_u})].
Note that $\veps$ is dimensionless (in units $c=1$) and that it shares the
same positiveness property (\ref{e:ges:E_positive_M_I}) as $E$:
\begin{prop}[positivity of the specific conserved energy]
The specific conserved energy $\veps$ is positive as soon as the
timelike geodesic $\Li$ has some part
in $\M_{\rm I}$, i.e. some part with $r>2m$:
\be \label{e:ges:eps_positive_M_I}
    \Li \cap \M_{\rm I} \not= \varnothing \quad \Longrightarrow \quad \veps > 0 .
\ee
\end{prop}
On the contrary, $\ell$ can be either positive, zero or negative, depending
on the variation of $\ph$ along $\Li$, as was already noticed above for $L$.

We note that Eq.~(\ref{e:ges:1d_motion_timelike}) has the shape of the
first integral of the
1-dimensional motion of a non-relativist particle in the potential
$V_{\ell}$ (called hereafter the \defin{effective potential}\index{effective!potential!timelike geodesic}), the term $1/2 \, (\D r/\D\tau)^2$ being interpreted as
the kinetic energy per unit mass, $V_{\ell}(r)$ as the potential
energy per unit mass and the constant right-hand side $(\veps^2-1)/2$ as the total
mechanical energy per unit mass.
\begin{remark} \label{r:ges:V_eff_Newt}
The effective potential (\ref{e:ges:V_eff_timelike}) differs from its
non-relativistic (Newtonian) counterpart only by the factor $1-2m/r$ instead of
$1$. This difference plays an important role for small values of $r$,
leading to some orbital instability, as we shall see in Sec.~\ref{s:ges:circular_orbits}.
\end{remark}

In $\M_{\rm I}$, where the Killing vector $\w{\xi}$ is timelike,
we may introduce the \defin{static observer}\index{static!observer!in Schwarzschild spacetime}\index{observer!static --}
$\Obs$, whose 4-velocity $\w{u}_{\Obs}$ is collinear to $\w{\xi}$:
\be \label{e:ges:u_Obs_xi}
    \w{u}_{\Obs} = \left( 1 - \frac{2m}{r} \right) ^{-1/2} \, \w{\xi} ,
\ee
the proportionality coefficient ensuring that $\w{u}_\Obs\cdot\w{u}_\Obs=-1$
given that $\w{\xi}\cdot\w{\xi} = g_{tt} = - (1-2m/r)$.
We have then, from (\ref{e:ges:def_eps_ell}),
\be
    \veps = - \left( 1 - \frac{2m}{r} \right) ^{1/2} \w{u}_{\Obs}\cdot\w{u}
        = \Gamma \left( 1 - \frac{2m}{r} \right) ^{1/2} ,
\ee
where $\Gamma = - \w{u}_{\Obs}\cdot\w{u}$ is the Lorentz factor of $\mathscr{P}$
with respect to $\Obs$ (cf. Sec.~\ref{s:fra:measure}; in particular Eq.~(\ref{e:fra:def_Lorentz_factor})). We may express $\Gamma$ in terms of the norm $v$ of the
velocity of $\mathscr{P}$ with respect to $\Obs$, according to Eq.~(\ref{e:fra:Gam_V2}):
$\Gamma = (1-v^2)^{-1/2}$ and get
\be
    \veps = \left( 1 - v^2 \right) ^{-1/2} \left( 1 - \frac{2m}{r} \right) ^{1/2}  .
\ee
In the region $r\gg m$, we may perform a first order expansion, assuming
that  $\mathscr{P}$ moves at nonrelativistic velocity with respect to $\Obs$
($v\ll 1$), thereby obtaining:
\be \label{e:ges:veps_far}
    \encadre{ \veps - 1 \simeq \frac{1}{2} v^2 - \frac{m}{r} } \qquad
    \left(r\gg m\quad\mbox{and}\quad v\ll 1\right).
\ee
We recognize in the right-hand side the \emph{Newtonian mechanical energy per
unit mass}\index{Newtonian!mechanical energy}\index{mechanical!energy} of particle $\mathscr{P}$ with respect to observer $\Obs$, who can
then be considered as an inertial observer, $v^2/2$ being the kinetic energy per unit mass
and $-m/r$ the gravitational potential energy per unit mass.

\begin{figure}
\centerline{\includegraphics[width=0.8\textwidth]{ges_eff_pot.pdf}}
\caption[]{\label{f:ges:eff_pot} \footnotesize
Effective potential $V_{\ell}(r)$ governing the $r$-part of the
motion along a timelike geodesic in
Schwarzschild spacetime via Eq.~(\ref{e:ges:1d_motion_timelike}).
The vertical dashed line marks $r=2m$, i.e. the
location of the event horizon.
The numerical value $\ell/m=3.4641$ is that of the critical
specific angular momentum (\ref{e:ges:ell_crit}).}
\end{figure}

\begin{figure}
\centerline{\includegraphics[width=0.8\textwidth]{ges_eff_pot_zoom.pdf}}
\caption[]{\label{f:ges:eff_pot_zoom} \footnotesize
Same as Fig.~\ref{f:ges:eff_pot}, but with a zoom in along the $y$-axis
and a zoom out along the $x$-axis. The dots mark the mimima of
$V_{\ell}$, locating stable circular orbits.}
\end{figure}

The profile  of $V_{\ell}(r)$ for selected values of $\ell$ is
plotted in Figs.~\ref{f:ges:eff_pot} and \ref{f:ges:eff_pot_zoom} .
Its extrema are given by
$\D V_{\ell}/\D r=0$, which is equivalent to
\[
    m r^2 - \ell^2 r + 3 \ell^2 m = 0 .
\]
This quadratic equation admits real roots iff $|\ell| \geq \ell_{\rm crit}$,
with
\be \label{e:ges:ell_crit}
   \encadre{ \ell_{\rm crit} = 2\sqrt{3}\, m \simeq 3.464102 \, m}.
\ee
For $|\ell| \geq \ell_{\rm crit}$, the two roots are
\be
    r_{\rm max} = \frac{\ell}{2m} \left( \ell -
    \sqrt{\ell^2 - \ell_{\rm crit}^2} \right)
    \qquad\mbox{and}\qquad
    r_{\rm min} = \frac{\ell}{2m} \left( \ell +
    \sqrt{\ell^2 - \ell_{\rm crit}^2} \right) ,
\ee
corresponding respectively to a maximum of $V_{\ell}$ and a minimum
of $V_{\ell}$, hence the indices ``max'' and ``min''. Note that
$r_{\rm max} \leq r_{\rm min}$.
In the marginal case $|\ell| = \ell_{\rm crit}$, the two roots
coincide and correspond to an inflection point of $V_{\ell}$ (the circled
dot in Fig.~\ref{f:ges:eff_pot_zoom}).

For $|\ell| < \ell_{\rm crit}$, there is no extremum and
$V_{\ell}$ is a strictly increasing function of $r$.


To get a full solution in terms of the Schwarzschild-Droste coordinates,
once Eq.~(\ref{e:ges:1d_motion_timelike}) is solved for $r(\tau)$,
one has still to solve Eqs.~(\ref{e:ges:dot_t}) and (\ref{e:ges:dot_ph}),
which can be rewritten in terms of the proper time $\tau$ as
\be \label{e:ges:Dt_Dtau}
    \frac{\D t}{\D \tau} = \veps \left(1 - \frac{2m}{r(\tau)} \right) ^{-1} ,
\ee
\be \label{e:ges:Dph_Dtau}
    \frac{\D\ph}{\D\tau} = \frac{\ell}{r(\tau)^2}  .
\ee


\subsection{Radial free fall} \label{s:ges:radial_free_fall}

\subsubsection{Generic case}

The radial geodesics correspond to a vanishing conserved angular momentum:
$\ell = 0$. Indeed, setting $\ell=0$ in Eq.~(\ref{e:ges:dot_ph}) yields
$\ph = \mathrm{const}$, which defines a purely radial trajectory in the
plane $\th=\pi/2$.
The effective potential (\ref{e:ges:V_eff_timelike}) reduces then
to $V_{\ell}(r) = - m/r$, so that the equation of radial motion
(\ref{e:ges:1d_motion_timelike}) becomes
\be \label{e:ges:radial_motion}
    \frac{1}{2} \left( \frac{\D r}{\D \tau} \right) ^2
        - \frac{m}{r} = \frac{\veps^2 - 1}{2} .
\ee
This equation is identical to that governing radial free fall in the gravitational field generated by a mass $m$ in Newtonian gravity.
The solution is well known and
depends on the sign of the ``mechanical energy'' in the right-hand
side, i.e. of the position of $\veps$ with respect to $1$:
\begin{itemize}
\item if $\veps>1$, the solution
is given in parameterized form (parameter $\eta$) by
\be \label{e:ges:sol_E_pos}
    \left\{ \begin{array}{l}
    \displaystyle\tau = \frac{m}{(\veps^2 - 1)^{3/2}} \left( \sinh\eta - \eta \right)
        + \tau_0 \\[2ex]
    \displaystyle r = \frac{m}{\veps^2 - 1} \left( \cosh\eta - 1 \right),
    \end{array} \right.
\ee
\item if $\veps=1$, the solution is
\be \label{e:ges:sol_E_zero}
    r(\tau) =  \left( \frac{9 m}{2} (\tau -\tau_0)^2 \right) ^{1/3} ,
\ee
\item if $\veps<1$, the solution
is given in parameterized form (parameter $\eta$) by
\be \label{e:ges:sol_E_neg}
    \left\{ \begin{array}{l}
    \displaystyle\tau =  \frac{m}{(1 - \veps^2) ^{3/2}} \left( \eta + \sin\eta \right)
    + \tau_0  \\[2ex]
    \displaystyle r = \frac{m}{1 - \veps^2} \left( 1 + \cos\eta \right),
    \end{array} \right.
\ee
\end{itemize}
In the above formulas, $\tau_0$ is a constant; for $\veps\geq \mu$,
$\tau_0$ is the value of $\tau$ for which $r\to 0$, while
for $\veps<1$, it is the value of $\tau$ at which $r$ takes its maximal value.

\subsubsection{Radial free fall from rest}

Let us focus on the radial free fall from rest, starting at some
position $r=r_0$ at $\tau=0$. Starting from rest
means $\D r/\D\tau = 0$ at $\tau=0$. The equation of radial motion
(\ref{e:ges:radial_motion}) leads then to $-m/r_0 = (\veps^2 - 1)/2$,
or equivalently
\be \label{e:ges:beps2}
         \veps^2 = 1 - \frac{2m}{r_0} .
\ee
The right-hand side of this equation must be non-negative. This implies
$r_0\geq 2m$. We recover the fact that one cannot be momentarily at rest
(in terms of $r$) if $r_0<2m$, for $r$ has to decrease along any causal
geodesic in the black hole region $\M_{\rm II}$
(cf. Sec.~\ref{s:ges:eq_to_be_solved}).

Equation~(\ref{e:ges:beps2}) implies $\veps < 1$, i.e. $E < \mu$. The solution
is thus given by Eq.~(\ref{e:ges:sol_E_neg}); expressing $1-\veps^2$
in it via (\ref{e:ges:beps2}), we get
\be \label{e:ges:sol_radial_infall}
    \encadre{
    \left\{ \begin{array}{l}
    \displaystyle\tau = \sqrt{\frac{r_0^3}{8 m}}  \left( \eta + \sin\eta \right) \\[3ex]
    \displaystyle r = \frac{r_0}{2} \left( 1 + \cos\eta \right)
    \end{array} \right. }
    \qquad 0 \leq \eta \leq \pi ,
\ee
where the range of $\eta$ is such that $r=r_0$ for $\tau=0$ ($\eta=0$) and
$r$ decays to $0$ when $\eta\to \pi$. The function $r(\tau)$ resulting
from (\ref{e:ges:sol_radial_infall}) is depicted in
Fig.~\ref{f:ges:radial_infall_tau}.

\begin{figure}
\centerline{\includegraphics[width=0.8\textwidth]{ges_radial_infall_tau.pdf}}
\caption[]{\label{f:ges:radial_infall_tau} \footnotesize
Coordinate $r$ as a function of the proper time $\tau$
for the radial free fall from rest, for various initial values $r_0$ of $r$.}
\end{figure}

The solution for $t=t(\tau)$ is obtained by combining $\D t/\D\tau$ as expressed
by (\ref{e:ges:Dt_Dtau}) and $\D\tau/\D\eta$ deduced from (\ref{e:ges:sol_radial_infall}):
\[
    \frac{\D\tau}{\D\eta} = \sqrt{\frac{r_0^3}{8 m}}  \left( 1 + \cos\eta \right)
        = \sqrt{\frac{r_0}{2m}} \, r .
\]
We get
\[
    \frac{\D t}{\D\eta} =  \frac{\D t}{\D\tau} \frac{\D\tau}{\D\eta}
        = \veps \sqrt{\frac{r_0}{2m}} \, r \left(1 - \frac{2m}{r} \right) ^{-1}
        = \sqrt{\frac{r_0}{2m} - 1} \; \frac{r^2}{r - 2m} ,
\]
where we have used (\ref{e:ges:beps2}) and $\veps > 0$ [Eq.~(\ref{e:ges:eps_positive_M_I})]
to write
$\veps = \sqrt{1-2m/r_0}$.
Substituting $r$ from Eq.~(\ref{e:ges:sol_radial_infall}), we get
\[
    \frac{\D t}{\D\eta} =   \frac{r_0}{2} \sqrt{\frac{r_0}{2m} - 1} \;
        \frac{(1+\cos\eta)^2}{1+\cos\eta - 4m/r_0} .
\]
This equation can be integrated to (cf. the \textsf{SageMath} computation in Sec.~\ref{s:sam:ges_radial_free_fall})
\be \label{e:ges:sol_t_radial_infall}
   \encadre{ t = 2m \left\{ \sqrt{\frac{r_0}{2m} - 1} \left[ \eta + \frac{r_0}{4m}
    (\eta + \sin\eta) \right]
    + \ln \left| \frac{\sqrt{\frac{r_0}{2m} - 1} + \tan\frac{\eta}{2}}{\sqrt{\frac{r_0}{2m} - 1} - \tan\frac{\eta}{2}} \right| \right\} } ,
\ee
where we have assumed $t=0$ at $\tau=0$ ($\eta=0$).

\begin{figure}
\centerline{\includegraphics[height=0.5\textheight]{ges_radial_infall_t.pdf}\qquad\qquad
\includegraphics[height=0.5\textheight]{ges_radial_infall_IEF.pdf}}
\caption[]{\label{f:ges:radial_infall} \footnotesize
Radial free fall from rest, viewed in Schwarzschild-Droste coordinates $(t,r)$
(left) and in the ingoing Eddington-Finkelstein coordinates $(\ti,r)$ (right),
for various values $r_0$ of the coordinate $r$ at $\tau=0$. The grey area is the
black hole region $\M_{\rm II}$.}
\end{figure}


The solution of the radial free fall starting from rest at $r=r_0$ is
thus given in parametric form by Eqs.~(\ref{e:ges:sol_radial_infall}) and
(\ref{e:ges:sol_t_radial_infall}) and is represented in the left panel
of Fig.~\ref{f:ges:radial_infall}. It has been obtained in the Schwarzschild-Droste coordinates
$(t,r,\th,\ph)$, which are singular at the event horizon $\Hor$. So, one might wonder
if such a solution can describe the full infall, with the crossing of $\Hor$.
In particular, we notice that the differential equation for $t(\tau)$,
Eq.~(\ref{e:ges:Dt_Dtau}), is singular at $r(\tau)=2m$, i.e. on $\Hor$. The solution
$t(\eta)$, as given by Eq.~(\ref{e:ges:sol_t_radial_infall}), is singular at
$\eta=\eta_{\rm h}$, where
\be \label{e:ges:def_eta_star}
    \eta_{\rm h} := 2 \arctan \sqrt{\frac{r_0}{2m} - 1}
\ee
is precisely the value of $\eta$ yielding $r=2m$ in Eq.~(\ref{e:ges:sol_radial_infall})
[to see it, rewrite the second part of Eq.~(\ref{e:ges:sol_radial_infall}) as
$r=r_0\cos^2(\eta/2) = r_0/(1+\tan^2(\eta/2))$].
This singularity of $t(\eta)$ appears also clearly on Fig.~\ref{f:ges:radial_infall} (left panel).
On the other hand, the equation
for $r$, Eq.~(\ref{e:ges:radial_motion}), does not exhibit any pathology
at $r=2m$, nor its solution (\ref{e:ges:sol_t_radial_infall}).
Actually, had we started from the ingoing Eddington-Finkelstein (IEF) coordinates
$(\ti,r,\th,\ph)$, instead of the Schwarzschild-Droste ones, we would have
found\footnote{{\emph Exercise:} do it!}
exactly the same solution for $r(\tau)$ (which is not surprising since
$r$, considered as a scalar field on $\M$, is perfectly regular at $\Hor$).
The solution for $\ti(\tau)$ can be deduced from that for $t(\tau)$ by
the coordinate transformation law (\ref{e:sch:ti_t_r}). Noticing that and
$r/(2m) = \cos^2(\eta/2)/\cos^2(\eta_{\rm h}/2)$, we get
\bea
    \frac{r}{2m} - 1 & = & \frac{\cos^2(\eta/2)}{\cos^2(\eta_{\rm h}/2)} - 1
        = \cos^2(\eta/2) \left( \frac{1}{\cos^2(\eta_{\rm h}/2)}
        - \frac{1}{\cos^2(\eta/2)} \right) \nonumber \\
       & = &\cos^2(\eta/2) \left(\tan^2(\eta_{\rm h}/2) - \tan^2(\eta/2) \right) .
\eea
Using this identity, as well as (\ref{e:ges:def_eta_star}) to express
$\sqrt{r_0/(2m)-1}$ in the logarithm term of Eq.~(\ref{e:ges:sol_t_radial_infall}), the transformation law (\ref{e:sch:ti_t_r}) yields
\bea
    \ti & = & 2m \left\{ \sqrt{\frac{r_0}{2m} - 1} \left[ \eta + \frac{r_0}{4m}
    (\eta + \sin\eta) \right]
    + \ln \left| \frac{\tan\frac{\eta_{\rm h}}{2} + \tan\frac{\eta}{2}}{\tan\frac{\eta_{\rm h}}{2}- \tan\frac{\eta}{2}} \cos^2\frac{\eta}{2} \left( \tan^2\frac{\eta_{\rm h}}{2}
        - \tan^2\frac{\eta}{2} \right) \right| \right\} \nonumber \\
    & = & 2m \left\{ \sqrt{\frac{r_0}{2m} - 1} \left[ \eta + \frac{r_0}{4m}
    (\eta + \sin\eta) \right]
    + \ln \left| \cos^2\frac{\eta}{2}
    \left( \tan\frac{\eta_{\rm h}}{2} + \tan\frac{\eta}{2} \right) ^2 \right| \right\} \nonumber \\
    & = & 2m \left\{ \sqrt{\frac{r_0}{2m} - 1} \left[ \eta + \frac{r_0}{4m}
    (\eta + \sin\eta) \right]
    +  2 \ln \left( \cos\frac{\eta}{2} \tan\frac{\eta_{\rm h}}{2} + \sin \frac{\eta}{2} \right) \right\} . \nonumber
\eea
From this expression, we have $\ti= 4 m \ln \tan(\eta_{\rm h}/2)$ for $\eta=0$.
Now, we can change the origin of the IEF coordinate $\ti$ to ensure $\ti=0$ for $\eta=0$,
i.e. $\tau=0$.
We get then
\be \label{e:ges:sol_radial_infall_ti}
    \encadre{ \ti = 2m \left\{ \sqrt{\frac{r_0}{2m} - 1} \left[ \eta + \frac{r_0}{4m}
    (\eta + \sin\eta) \right]
    +  2 \ln \left[ \cos\frac{\eta}{2}  + \left(\frac{r_0}{2m} - 1\right)^{-1/2} \sin \frac{\eta}{2} \right] \right\} } .
\ee
This expression is perfectly regular for all values of $\eta$ in $[0,\pi]$, reflecting
the fact that the ingoing Eddington-Finkelstein coordinates cover all $\M$
in a regular way. The radial free fall solution in terms of $(\ti, r)$
is represented in the right panel of Fig.~\ref{f:ges:radial_infall}. We note the
smooth crossing of the event horizon $\Hor$.

\begin{figure}
\centerline{\includegraphics[height=0.4\textheight]{ges_infall_time.pdf}}
\caption[]{\label{f:ges:infall_time} \footnotesize
Elapsed proper time to reach the event horizon ($\tau_{\rm h}$, dashed curve)
and the central singularity ($\tau_{\rm f}$, solid curve), as a function
of the initial value of $r$ for a radial free fall from rest.}
\end{figure}

\begin{figure}
\centerline{\includegraphics[height=0.4\textheight]{ges_time_inside.pdf}}
\caption[]{\label{f:ges:time_inside} \footnotesize
Proper time spent inside the black hole region as a function
of the initial value of $r$ for a radial free fall from rest.
Note that $r_0=2m$ does not correspond to any asymptote but to the
finite value $\Delta\tau_{\rm in} = \pi\, m$ with a vertical tangent.
On the other side, there is an horizontal asymptote
$\Delta\tau_{\rm in} \to 4m/3$ for $r_0\to +\infty$.}
\end{figure}

In view of Eq.~(\ref{e:ges:sol_radial_infall}), we may say that the radial
infall starts at $\eta=0$, for which $\tau=0$ and $r=r_0$, and terminates
at $\eta=\pi$, for which $r=0$, which means that the particle hits the curvature
singularity (cf. Sec.~\ref{s:sch:singularities}).
The final value of the particle's proper time is obtained by
setting $\eta=\pi$ in Eq.~(\ref{e:ges:sol_radial_infall}):
\be
    \encadre{\tau_{\rm f} = \frac{\pi}{2} \sqrt{\frac{r_0^3}{2 m}} } .
\ee
Similarly, the final value of $\ti$ is obtained by setting $\eta=\pi$ in (\ref{e:ges:sol_radial_infall_ti}):
\be
    \encadre{\ti_{\rm f} = 2 m \left[ \pi \sqrt{\frac{r_0}{2m} - 1 }
        \left( \frac{r_0}{4m} + 1 \right)
        -  \ln  \left( \frac{r_0}{2m} - 1 \right) \right] } .
\ee
As noticed above, the event horizon $\Hor$ is crossed at $\eta=\eta_{\rm h}$;
via (\ref{e:ges:sol_radial_infall}) and (\ref{e:ges:def_eta_star}), this corresponds to the following value
of the proper time:
\be
    \encadre{\tau_{\rm h} = \sqrt{\frac{r_0^3}{2m}}
        \left[ \arctan \sqrt{\frac{r_0}{2m} - 1} +
        \sqrt{\frac{2m}{r_0} \left( 1 - \frac{2m}{r_0} \right) } \right] } ,
\ee
while (\ref{e:ges:sol_radial_infall_ti}) leads to the following value
of the IEF coordinate $\ti$:
\be
    \encadre{\ti_{\rm h} = 2m \left[ 2\left(1+\frac{r_0}{4m}\right)
        \sqrt{\frac{r_0}{2m} - 1} \arctan \sqrt{\frac{r_0}{2m} - 1}
        + \frac{r_0}{2m} - 1 - \ln\frac{r_0}{2m} \right] } .
\ee
The variation of $\tau_{\rm h}$ and $\tau_{\rm f}$ with $r_0$ are depicted in
Fig.~\ref{f:ges:infall_time} and numerical values for $r_0=6m$ and
standard astrophysical black holes are provided in
Table~\ref{t:ges:time_free_fall}.

\begin{table}
\centerline{
\begin{tabular}{|c|c|c|c|c|}
\hline
$m$ & $r_{\rm S} = 2m$ & $\tau_{\rm h}$ & $\tau_{\rm f}$ & $\Delta\tau_{\rm in}$ \\[0.5ex]
\hline\hline
$15\, M_\odot$ (Cyg X-1) & $44.3 {\rm\; km}$ & $1.10 {\rm\; ms}$ & $1.21 {\rm\; ms}$ & $0.11 {\rm\; ms}$ \\[0.5ex]
\hline
$4.3\; 10^6 \, M_\odot$ (Sgr A*) & $12.7\; 10^6 {\rm\; km} = 0.085 {\rm\; au}$ & $5 {\rm\; min\; } 14 {\rm\; s}$ & $5 {\rm\; min\; } 46 {\rm\; s}$ & $32 {\rm\; s}$ \\[0.5ex]
\hline
$6\; 10^9 \, M_\odot$ (M87*) & $118 {\rm\; au}$ & $5.07 {\rm\; days}$ & $5.58 {\rm\; days}$ & $12 {\rm\; h\;} 17 {\rm\; min}$ \\[0.5ex]
\hline
\end{tabular}
}
\caption[]{\label{t:ges:time_free_fall} \footnotesize
Proper time to reach the event horizon ($\tau_{\rm h}$) and
the central curvature singularity ($\tau_{\rm f}$), as well as elapsed proper
time inside the black hole region ($\Delta\tau_{\rm in}$), when freely falling
from rest at $r_0 = 6m$. The numerical values are given for various
black hole masses $m$, corresponding to astrophysical objects:
the stellar black hole Cygnus X-1\index{Cygnus X-1} \cite{Orosz_al11,Gou_al14}, the
supermassive black hole at the center of our galaxy
(Sagittarius A*)\index{Sagittarius A*}  \cite{GenzeEG10,Johan16} and the supermassive black hole M87* in the nucleus of the galaxy Messier 87\index{M87*} \cite{Gebha_al11,EHT19a}.}
\end{table}

The proper time spent inside the black hole is
\be
    \encadre{\Delta\tau_{\rm in}  = \tau_{\rm f} - \tau_{\rm h}
    = \sqrt{\frac{r_0^3}{2m}}
        \left[ \frac{\pi}{2} - \arctan \sqrt{\frac{r_0}{2m} - 1} -
        \sqrt{\frac{2m}{r_0} \left( 1 - \frac{2m}{r_0} \right) } \right] } .
\ee
It varies between $\pi m$ ($r_0 \to 2m$) and $4m/3$ ($r_0 \to +\infty$)
(cf. Fig.~\ref{f:ges:time_inside} and Sec.~\ref{s:sam:ges_radial_free_fall}
for the computation of $\lim_{r_0\to +\infty} \Delta\tau_{\rm in}$).
Numerical values for astrophysical black holes are provided in
Table~\ref{t:ges:time_free_fall}.

%%%%%%%%%%%%%%%%%

\subsection{Circular orbits}\index{circular!orbit}\index{orbit!circular --}
\label{s:ges:circular_orbits}

Circular orbits are defined as timelike geodesics with $r = \mathrm{const}$.
We have then $\D r/\D\tau = 0$ and $\D^2 r/\D\tau^2 = 0$, so that
Eq.~(\ref{e:ges:1d_motion_timelike}) implies
\begin{subequations}
\begin{align}
& V_{\ell}(r) = \frac{\veps^2 - 1}{2} \label{e:ges:Veff_circ} \\
& \frac{\D V_{\ell}}{\D r} = 0 . \label{e:ges:dVeffdr_zero}
\end{align}
\end{subequations}
Given the expression (\ref{e:ges:V_eff_timelike}) of $V_{\ell}$,
Eq.~(\ref{e:ges:dVeffdr_zero}) is equivalent to
\be \label{e:ges:eq_r_ell_circ}
    m r^2 - \ell^2 r + 3 \ell^2 m = 0 .
\ee
As already noticed in Sec.~\ref{s:ges:eff_pot_timelike}, this
quadratic equation in $r$ admits two real roots iff $|\ell| \geq \ell_{\rm crit}$,
with $\ell_{\rm crit} = 2\sqrt{3}\, m$ [Eq.~(\ref{e:ges:ell_crit})],
which are
\be \label{e:ges:r_circ_ell}
    \encadre{ r_{\rm circ}^\pm(\ell) = \frac{\ell}{2m} \left( \ell \pm
    \sqrt{\ell^2 - \ell_{\rm crit}^2} \right) }.
\ee
$r_{\rm circ}^+(\ell)$ corresponds to a minimum of the effective potential
$V_{\ell}$ and thus
to a stable orbit (see the dots in Fig.~\ref{f:ges:eff_pot_zoom}),
while $r_{\rm circ}^-(\ell)$ corresponds to a
maximum of $V_{\ell}$ and thus to an unstable orbit.
When $\ell$ varies from $\ell_{\rm crit}$ to $+\infty$,
$r_{\rm circ}^+(\ell)$ increases from $6 m$ to $+\infty$, while
$r_{\rm circ}^-(\ell)$ decreases from $6 m$ to $3 m$ (cf. Fig.~\ref{f:ges:ell_circ_orbit}). We conclude that
\begin{prop}[circular orbits]
Circular orbits in Schwarzschild spacetime exist for all values of $r>3m$.
Those with $r<6m$ are unstable and those with $r>6m$ are stable. The marginal
case $r=6m$ is called the
\defin{innermost stable circular orbit}\index{stable!circular orbit}\index{innermost!stable circular orbit}\index{circular!orbit!innermost stable --}, often abridged as \defin{ISCO}\index{ISCO}.
\end{prop}

\begin{remark}
In the Newtonian spherical gravitational field generated by a point mass $m$,
there is no unstable orbit, and thus no ISCO. The existence of unstable orbits
in the relativistic case can be understood by the extra term in
the effective potential $V_{\ell}(r)$ (cf. Remark~\ref{r:ges:V_eff_Newt} on p.~\pageref{r:ges:V_eff_Newt}), which
adds the attractive part $-\ell^2 m/r^3$ to the two
terms constituting the Newtonian potential: $-m/r$ (attractive) and
$\ell^2/(2 r^2)$ (repulsive). The latter is responsible for the infinite ``centrifugal
barrier''\index{centrifugal!barrier} at small $r$ in the Newtonian problem, leading
always to a minimum of $V_\ell(r)$ and thus to a stable circular orbit.
In the relativistic case, for $r$ small enough, the attractive term, which is $O(r^{-3})$, dominates over
the centrifugal one, which is only $O(r^{-2})$.
Equivalently, we may say that the ``centrifugal barrier'' is weakened by the
factor $1-2m/r$ (cf. the expression (\ref{e:ges:V_eff_timelike}) of $V_\ell(r)$)
and ceases to exist for small values of $|\ell|$ (i.e.
$|\ell|<\ell_{\rm crit}$).
\end{remark}

\begin{figure}
\centerline{\includegraphics[height=0.4\textheight]{ges_ell_circ_orbit.pdf}}
\caption[]{\label{f:ges:ell_circ_orbit} \footnotesize
Specific conserved angular momentum $\ell = L/\mu$
on circular orbits as a function of the
orbit circumferential radius $r$. The dashed part of the curve
corresponds to unstable orbits ($r=r_{\rm circ}^-(\ell)$, as given by Eq.~(\ref{e:ges:r_circ_ell})), while
the solid part corresponds to stable orbits ($r=r_{\rm circ}^+(\ell)$).
The minimal value of $\ell$ is $\ell_{\rm crit} = 2\sqrt{3}\, m \simeq 3.46\, m$.}
\end{figure}

From Eq.~(\ref{e:ges:eq_r_ell_circ}), we can easily express $\ell$
as a function of $r$ on a circular orbit:
\be \label{e:ges:ell_r_circ}
    \encadre{ |\ell| = r \sqrt{\frac{m}{r - 3m}} } .
\ee
This function is represented in Fig.~\ref{f:ges:ell_circ_orbit} (for $\ell > 0$).

If we substitute (\ref{e:ges:ell_r_circ}) for $\ell$ in
the expression (\ref{e:ges:V_eff_timelike}) of $V_{\ell}$
and use Eq.~(\ref{e:ges:Veff_circ}), we obtain the value of the
specific conserved energy along a circular orbit, in terms of $r$:
\be \label{e:ges:eps_r_circ}
    \encadre{ \veps = \frac{r - 2m}{\sqrt{r(r-3m)}} } .
\ee
This function is represented in Fig.~\ref{f:ges:ener_circ_orbit}.
The minimal value of $\veps$ is achieved for $r=6m$, i.e. at the ISCO:
\be
    \encadre{ \mathrm{min}\,  \veps = \frac{2\sqrt{2}}{3} \simeq 0.9428} .
\ee
From Fig.~\ref{f:ges:ener_circ_orbit}, we notice that
\be \label{e:ges:r_marginally_bound}
    r > 4 m \iff \veps < 1 \iff E < \mu .
\ee
This corresponds to \defin{bound orbits}\index{bound!orbit}\index{orbit!bound --},
i.e. to geodesics that, if slightly perturbed, cannot reach the asymptotically
flat region $r\gg 2m$, since $E \geq \mu$ there. Indeed, when
$r\to +\infty$, the Killing vector $\w{\xi}$ can be interpreted
as the 4-velocity of some asymptotically inertial observer (at rest with
respect to the black hole) and $E$ is the particle energy measured by
that observer; the famous Einstein relation (\ref{e:fra:E_Gam_m}) is then
$E = \Gamma \mu$, where $\Gamma$ is the Lorentz factor
of the particle with respect to the observer.
Since $\Gamma \geq 1$ [Eq.~(\ref{e:fra:Gam_V2})], we have obviously\footnote{Similarly,
the radial-motion solutions (\ref{e:ges:sol_E_pos})-(\ref{e:ges:sol_E_zero}), which
allow for $r\to +\infty$, have $E \geq\mu$, while the solution (\ref{e:ges:sol_E_neg}),
which is relevant for a free fall from rest, has $E < \mu$.} $E \geq \mu$.
For this reason, the circular orbit at $r=4m$ is called the
\defin{marginally bound circular orbit}\index{marginally!bound!circular orbit}\index{orbit!marginally bound --}\index{marginally!bound!geodesic}. Note that the marginally bound
circular orbit is unstable, since it has $r<6m$.

\begin{figure}
\centerline{\includegraphics[height=0.4\textheight]{ges_ener_circ_orbit.pdf}}
\caption[]{\label{f:ges:ener_circ_orbit} \footnotesize
Specific conserved energy $\veps = E/\mu$
on circular orbits as a function of the
orbit circumferential radius $r$.
The dashed part of the curve
corresponds to unstable orbits,
while the solid part corresponds to stable ones.
The horizontal red line $\veps=1$ marks the limit of bound orbits.}
\end{figure}

The track of circular orbits in the $(\ell,\veps)$ plane is depicted in
Fig.~\ref{f:ges:circ_eps_ell}. The ISCO, which is a minimum for both $\veps$
and $\ell$, appears as a cusp point.

\begin{figure}
\centerline{\includegraphics[height=0.4\textheight]{ges_circ_eps_ell.pdf}}
\caption[]{\label{f:ges:circ_eps_ell} \footnotesize
Circular orbits in the $(\ell,\veps)$ plane.
The solid (resp. dashed) curve
corresponds to stable (resp. unstable) orbits. The ISCO is located at the cusp point.}
\end{figure}

The \defin{angular velocity}\index{angular!velocity} of a circular orbit $\Li$ is defined by
\be \label{e:ges:def_Omega}
    \Omega := \left. \frac{\D\ph}{\D t} \right| _\Li = \frac{u^\ph}{u^t} ,
\ee
where $u^\ph = \D\ph/\D\tau$ and $u^t = \D t/\D\tau$ are the only
nonzero components w.r.t. Schwarzschild-Droste coordinates
of the 4-velocity $\w{u}$ along the worldline
$\Li$. It follows from (\ref{e:ges:def_Omega}) that
$\Omega$ enters into the linear combination of the two Killing
vectors $\w{\xi}$ and $\w{\eta}$ expressing the 4-velocity on a circular
orbit according to
\be
    \w{u} = u^t\left( \w{\xi} + \Omega \w{\eta} \right) .
\ee
We have the following nice physical interpretation:
\begin{prop}[angular velocity measured at infinity]
The quantity $\Omega$ defined by Eq.~(\ref{e:ges:def_Omega}) is nothing but
the angular velocity of the orbiting particle $\mathscr{P}$ monitored by a infinitely
distant static observer $\Obs$.
\end{prop}
\begin{proof}
Suppose that $\Obs$ is located at fixed coordinates
$(r,\th,\ph)=(r_\Obs,\pi/2,0)$ with $r_\Obs\gg m$ and that
$\mathscr{P}$ emits a photon at the event
$(t_1, r, \pi/2, 0)$ along a \emph{radial} null geodesic. This photon is received
by $\Obs$ at $t=t'_1$. After one orbit, at the event $(t_2, r, \pi/2, 2\pi)$,
$\mathscr{P}$ emits a second photon in the radial direction, which is received
at $t=t'_2$ by $\Obs$. According to the definition (\ref{e:ges:def_Omega})
of $\Omega$, we have
\[
    2\pi = \Omega(t_2 - t_1) .
\]
On the other hand, since $r_\Obs\gg m$, the proper time of $\Obs$ is $t$, so that the
angular velocity measured by $\Obs$ is
\[
    \Omega_\Obs = \frac{2\pi}{t'_2 - t'_1} .
\]
Now, since $t$ is the coordinate associated to
the spacetime invariance by time translation (stationarity), we have
necessarily $t'_2 - t_2 = t'_1 - t_1$
(the increment in coordinate $t$ for the second signal is the same as for the first one),
so that $t'_2 - t'_1 = t_2 - t_1$. Accordingly, the above two equations
combine to $\Omega_\Obs = \Omega$.
\end{proof}

By combining Eqs.~(\ref{e:ges:Dt_Dtau}) and (\ref{e:ges:Dph_Dtau}), we get
\[
    \Omega = \frac{1}{r^2} \left( 1 -\frac{2m}{r} \right) \frac{\ell}{\veps} .
\]
Substituting expression (\ref{e:ges:ell_r_circ}) for $\ell$ and
expression (\ref{e:ges:eps_r_circ}) for $\veps$, we obtain
\be \label{e:ges:Omega_m_r}
    \encadre{ \Omega = \pm \sqrt{\frac{m}{r^3}} } ,
\ee
with the $+$  (resp. $-$) sign for $\ell > 0$ (resp. $\ell < 0$).

\begin{table}
\centerline{
\begin{tabular}{|c|c|c|c|c|}
\hline
$m$ & $r_{\rm ISCO} = 6m$ & $\dfrac{\Omega_{\rm ISCO}}{2\pi}$ & $T_{\rm ISCO}$  & $T_{\mathscr{P},\rm ISCO}$ \\[1ex]
\hline\hline
$15\, M_\odot$ (Cyg X-1) & $133 {\rm\; km}$ & $147 {\rm\; Hz}$ &
$6.80 {\rm\; ms}$ &  $4.81 {\rm\; ms}$ \\[0.5ex]
\hline
$4.3\; 10^6 \, M_\odot$ (Sgr A*) & ${\displaystyle 38.1\; 10^6 {\rm\; km}\atop \displaystyle (0.255 {\rm\; au})}$ & $5.11\; 10^{-4} {\rm\; Hz}$
& $32 {\rm\; min}\; 37 {\rm\; s}$ &
 $23 {\rm\; min}\; 4 {\rm\; s}$  \\[0.5ex]
\hline
$6\; 10^9 \, M_\odot$ (M87*) & $355 {\rm\; au}$ & $3.66\; 10^{-7} {\rm\; Hz}$
& $31 {\rm\; d}\; 15 {\rm\; h}$ & $22 {\rm\; d}\; 9 {\rm\; h}$ \\[0.5ex]
\hline
\end{tabular}
}
\caption[]{\label{t:ges:freq_ISCO} \footnotesize
Values of various quantities at the ISCO for masses $m$ of some astrophysical black holes
(see Table~\ref{t:ges:time_free_fall} for details):
areal radius $r$, orbital frequency $\Omega_{\rm ISCO}/(2\pi)$,
orbital period seen from infinity $T_{\rm ISCO}$ and orbital period measured
by the orbiting observer/particle $T_{\mathscr{P},\rm ISCO}$ (proper time).}
\end{table}


\begin{remark}
This formula is identical to that of Newtonian gravity (Kepler's third law for circular
orbits) for all values of $r$. This is a mere coincidence, valid only for
Schwarzschild-Droste coordinates. Only for $r\gg m$, i.e. in the weak-field limit, this
agreement is physically meaningful; it can be then used to interpret the
parameter $m$ as the \emph{gravitational mass}\index{mass!gravitational --}\index{mass!parameter of Schwarzschild solution} of Schwarzschild spacetime,
as mentioned in Sec.~\ref{s:sch:mass_parameter}.
\end{remark}

\begin{remark}
$\Omega$ is not the orbital angular frequency experienced by the particle/observer
$\mathscr{P}$ on the circular orbit $\Li$, because the proper time of $\mathscr{P}$ is
$\tau$ and not $t$. The actual orbital frequency measured by $\mathscr{P}$
is
\[
    \Omega_{\mathscr{P}} = \frac{\D t}{\D \tau} \Omega = u^t \Omega ,
\]
with $u^t = \D t/\D\tau$ obtained from (\ref{e:ges:Dt_Dtau}) and
(\ref{e:ges:eps_r_circ}): $u^t = \sqrt{r/(r-3m)}$. Hence
\be
    \Omega_{\mathscr{P}} = \sqrt{\frac{r}{r-3m}}\, \Omega  =
       \pm \frac{1}{r}\sqrt{\frac{m}{r-3m}} .
\ee
Note that $|\Omega_{\mathscr{P}}|> |\Omega|$; in particular, at the ISCO ($r=6m$),
$\Omega_{\mathscr{P}} = \sqrt{2} \Omega$. The orbital period\index{orbital!period} measured by $\mathscr{P}$
is $T_{\mathscr{P}} = 2\pi/|\Omega_{\mathscr{P}}|$. Some ISCO values of $T_{\mathscr{P}}$
for astrophysical black holes are provided in Table~\ref{t:ges:freq_ISCO}.
\end{remark}


At the ISCO, $r=6m$ and formula (\ref{e:ges:Omega_m_r}) yields (for orbits with positive $\ell$)
\be
    \encadre{ \Omega_{\rm ISCO} = \frac{1}{6\sqrt{6}\, m} } .
\ee
Numerical values of $\Omega_{\rm ISCO}$ (actually the frequency $\Omega_{\rm ISCO}/(2\pi)$,
which is more relevant from an observational point of view) are provided in
Table~\ref{t:ges:freq_ISCO}.

\subsection{Other orbits}

Let us relax the assumption $r=\mathrm{const}$ and consider generic orbits $\Li$
obeying
\be
    | \ell | > \ell_{\rm crit} \qquad\mbox{and}\qquad
      0 < \veps < 1 .
\ee
The first condition ensures that the effective potential $V_{\ell}(r)$
takes the shape of a well in the region $r>2m$ (cf. Fig.~\ref{f:ges:eff_pot_bound})
and the second one that the particle $\mathscr{P}$ is trapped in this well.
Indeed, $0 < \veps < 1$  makes the right-hand of Eq.~(\ref{e:ges:1d_motion_timelike}) negative, so
that the region $r\to +\infty$, where $V_{\ell}(r)\to 0$,
cannot be reached. We have also argued in Sec.~\ref{s:ges:circular_orbits}
that $\veps < 1$ is forbidden in the region $r\to +\infty$ on physical
grounds [cf. the discussion below Eq.~(\ref{e:ges:r_marginally_bound})].

\begin{figure}
\centerline{\includegraphics[height=0.4\textheight]{ges_eff_pot_bound.pdf}}
\caption[]{\label{f:ges:eff_pot_bound} \footnotesize
Effective potential $V_{\ell}(r)$ for $\ell = 4.2 m$ (one of the values displayed
in Figs.~\ref{f:ges:eff_pot}  and \ref{f:ges:eff_pot_zoom}).
The horizontal red line marks $V_{\ell}(r) = (\veps^2-1)/2$
with $\veps = 0.973$, leading to $r_{\rm per} = 9.058\, m$ and $r_{\rm apo} = 25.634\, m$.
The corresponding orbit is shown in Fig.~\ref{f:ges:orbit_e973_l42}.
}
\end{figure}

\begin{figure}
\centerline{\includegraphics[height=0.4\textheight]{ges_orbit_e973_l42.pdf}}
\caption[]{\label{f:ges:orbit_e973_l42} \footnotesize
Timelike geodesic with $\veps = 0.973$ and $\ell=4.2m$ (same values as in
Fig.~\ref{f:ges:eff_pot_bound}), plotted in terms of the coordinates
$(x,y) := (r\cos\ph, r\sin\ph)$. The dotted circles correspond to $r=r_{\rm per}$
(periastron) and $r=r_{\rm apo}$ (apoastron). The grey disk indicates the
black hole region $r<2m$. \textsl{[Figure produced with the notebook \ref{s:sam:ges_orbits}]}
}
\end{figure}

\begin{figure}
\centerline{\includegraphics[width=0.48\textwidth]{ges_orbit_e967_l42.pdf}\quad
\includegraphics[width=0.48\textwidth]{ges_orbit_e990_l42.pdf}}
\caption[]{\label{f:ges:orbit_e967_990_l42} \footnotesize
Timelike geodesics with the same value of $\ell$ as in Fig.~\ref{f:ges:orbit_e973_l42} ($\ell=4.2m$), but for different values of $\veps$: $\veps = 0.967$ (left) and $\veps=0.990$ (right).
Note that the left and right figures have different scales. \textsl{[Figure produced with the notebook \ref{s:sam:ges_orbits}]}
}
\end{figure}


In the potential well, the $r$-coordinate along $\Li$ varies between two
extrema: a minimum $r_{\rm per}$, for \defin{periastron}\index{periastron} (or \defin{pericenter}\index{pericenter} or \defin{periapsis}\index{periapsis}), and
a maximum
$r_{\rm apo}$, for \defin{apoastron}\index{apoastron}
(or \defin{apocenter}\index{apocenter} or \defin{apoapsis}\index{apoapsis}) (cf. Figs.~\ref{f:ges:eff_pot_bound} and \ref{f:ges:orbit_e973_l42}). Being extrema of $r(\tau)$,
the values of $r_{\rm per}$ and $r_{\rm apo}$ are obtained by setting
$\D r /\D\tau =0$ in Eq.~(\ref{e:ges:1d_motion_timelike}), which leads to
\be
    \left( 1 - \frac{2m}{r} \right) \left( 1 + \frac{\ell^2}{r^2} \right)
        = \veps^2 .
\ee
This is a cubic equation in $r^{-1}$, which has three real positive roots, corresponding
to the three intersections of the curve $V_{\ell}(r)$ with the horizontal
line at $(\veps^2-1)/2$ in Fig.~\ref{f:ges:eff_pot_bound}. However, the smaller
root has to be disregarded as a periastron since it would lead to a motion with
$V_{\ell}(r) > (\veps^2-1)/2$, which is forbidden by Eq.~(\ref{e:ges:1d_motion_timelike}).

We get, from Eqs.~(\ref{e:ges:1d_motion_timelike})-(\ref{e:ges:V_eff_timelike}),
\be
    \frac{\D r}{\D\tau} = \pm
    \sqrt{ \veps^2 - \left(1 - \frac{2m}{r}\right) \left( 1 + \frac{\ell^2}{r^2} \right) } .
\ee

The equation governing the trajectory of $\Li$ in the orbital plane is
Eq.~(\ref{e:ges:DuDph_trajectories}), which can be recast as
\be  \label{e:ges:dudph_timelike}
    \frac{\D u}{\D\ph}  = \pm \sqrt{ 2 u^3 - u^2 + 2 \left( \frac{m}{\ell} \right)^2 u
    + \left( \frac{m \veps}{\ell} \right)^2 - \left( \frac{m}{\ell} \right)^2  } .
\ee
Let us recall that $u := m/r$ and that methods for solving this
differential equation have been briefly discussed in Sec.~\ref{s:ges:trajectories}.

\begin{remark}
Far from the black hole, i.e. in the region $r\gg m$, one can easily recover
the Newtonian orbits from Eq.~(\ref{e:ges:dudph_timelike}). Indeed, according to
Eq.~(\ref{e:ges:veps_far}), $\veps = 1 + \veps_0$, where $\veps_0 = v^2/2 -m/r$  is
the Newtonian mechanical energy per unit mass. It obeys $|\veps_0|\ll 1$, so that $\veps^2 \simeq 1 + 2\veps_0$.
Moreover, for $r\gg m$, $u\ll 1$ and we can neglect the $u^3$ term in front
of the $u^2$ one in Eq.~(\ref{e:ges:dudph_timelike}). Hence Eq.~(\ref{e:ges:dudph_timelike})
reduces to
\be \label{e:ges:dudph_Newt_prov}
     \frac{\D u}{\D\ph} \simeq \pm \sqrt{ 2 \left( \frac{m}{\ell} \right)^2 \veps_0 - u^2
     + 2 \left( \frac{m}{\ell} \right)^2 u } .
\ee
Let us introduce the constants
\be \label{e:ges:def_p_e}
    p := \frac{\ell^2}{m} \qquad\mbox{and}\qquad
    e := \sqrt{1 + 2 \frac{\veps_0 \ell^2}{m^2}} .
\ee
Then Eq.~(\ref{e:ges:dudph_Newt_prov}) can be rewritten as
\[
    \frac{\D \ph}{\D u} = \pm \frac{\frac{p}{m e}}{\sqrt{1 - \left(\frac{\frac{p}{m}u-1}{e}\right) ^2}} ,
\]
which is readily integrated into
\[
    \ph = \pm \arccos \left(\frac{\frac{p}{m}u-1}{e}\right) + \ph_0 ,
\]
where $\ph_0$ is a constant. We have then $\frac{p}{m} u = 1 + e\cos(\ph-\ph_0)$, or
equivalently,
\be \label{e:ges:sol_Newt_ellipse}
    r = \frac{p}{1 + e\cos(\ph-\ph_0)}
\ee
Assuming a bound orbit, we have $\veps < 1$, which implies $\veps_0 < 0$
and, via Eq.~(\ref{e:ges:def_p_e}), $e<1$. We recognize then in (\ref{e:ges:sol_Newt_ellipse}) the equation
of an ellipse of eccentricity $e$ and semi-latus rectum $p$. Hence Keplerian orbits\index{Keplerian!orbit}\index{orbit!Keplerian --} are recovered for $r\gg m$, as they should.
\end{remark}

Generic bound orbits differ from the Keplerian ellipses by the fact that
the variation of $\ph$ between two successive \defin{periastron passages}\index{periastron!passage},
i.e. two events along the worldline of $\mathscr{P}$ for which $r=r_{\rm per}$,
is strictly larger than $2\pi$. This phenomenon is called \defin{periastron advance}\index{periastron!advance}
and causes the orbits to be not closed, as illustrated in Figs.~\ref{f:ges:orbit_e973_l42} and
\ref{f:ges:orbit_e967_990_l42}.


\begin{hist} \label{h:ges:geod}
The equations of geodesic motion (\ref{e:ges:dot_t})-(\ref{e:ges:dot_r_square}),
as well as Eq.~(\ref{e:ges:DuDph_trajectories}) for the trajectories in the
orbital plane, have been first given by Karl Schwarzschild\index[pers]{Schwarzschild, K.}
himself in January 1916 in the very same article \cite{Schwa1916} in which he presented
his famous solution\footnote{Equation~(\ref{e:ges:DuDph_trajectories}) for the
trajectories is Eq.~(18) in the article \cite{Schwa1916}, the
link between Schwarzschild's notations and ours being $R=r$, $x=u/m$, $c=L$,
$1=E$, $h=\mu^2$.}. Schwarzschild discussed only the weak field limit of these
equations, to recover the Mercury's perihelion advance computed by
Einstein\index[pers]{Einstein, A.}.

It is quite remarkable that the general solution to the geodesic motion
in Schwarzschild spacetime\footnote{More precisely: in the $\M_{\rm I}$ region
of Schwarzschild spacetime.} has been given as early as May 1916 by
Johannes Droste\index[pers]{Droste, J.}, in the same article \cite{Drost1917} in which
he derived the Schwarzschild solution, independently of Karl Schwarzschild
(cf. historical note on p.~\pageref{h:sch:Schwarzschild_sol}).
Droste derived the equations of geodesic motion in Schwarzschild-Droste coordinates
$(t, r, \th,\ph)$, using $t$ as the parameter
along the geodesics, as well as the equation governing
the trajectories in the orbital plane\footnote{The link between Droste's notations
and ours is $\alpha=2m$, $x=2u$, $z = 2u - 1/3$, $A = \mu^2/E^2$ and $B= L/E$.}.
He gave the solutions for the trajectories
in terms of the Weierstrass elliptic function\index{Weierstrass elliptic function}\index{elliptic!Weierstrass -- function} $\wp$ (cf. Sec.~\ref{s:ges:trajectories}).
He classified the solutions in terms of the roots of the cubic polynomial
in $v$ that appears in the right-hand side of Eq.~(\ref{e:ges:dvdtph_cubic}).
Droste noticed that it takes an infinite amount of coordinate time $t$ for a particle
to reach $r=2m$ (cf. left panel of Fig.~\ref{f:ges:radial_infall}) and
he concluded incorrectly that ``a moving particle outside the sphere $r=2m$ can
never pass that sphere''. He missed that this is only a coordinate effect,
reflecting the pathology of Schwarzschild-Droste coordinates at the horizon.
One shall keep in mind that in 1916, general relativity was just in its infancy
and disentangling coordinate artefacts from physical effects was not so obvious,
especially regarding time.
One can be amused by the fact that the first analysis of the static
black hole of general relativity made the black hole appear, not as an object
from which no particle may escape, but as an object into which no particle
may penetrate...

Finally, it is worth mentioning two early detailed studies of geodesic motion
of massive particles in Schwarzschild spacetime: one by
by Carlo De Jans\index[pers]{De Jans, C.} in 1923 \cite{DeJan1923} and other one by
Yusuke Hagihara\index[pers]{Hagihara, Y.} in 1931 \cite{Hagih1931}.
For a detailed account about the history of geodesic motion in
Schwarzschild spacetime see Ref.~\cite{Eisen87}.
\end{hist}