\magnification=1200
%\raggedbottom
\baselineskip=16pt
\parindent=20pt
\parskip=5pt
\def\rv{{\bf r}}
\def\xv{{\bf x}}
\def\vv{{\bf v}}
\def\av{{\bf a}}
\def\Gv{{\bf G}}
\def\dt{\Delta t}
\def\d{\Delta}
\def\half{{1 \over 2}}
\centerline{\bf Simulating Relativistic Orbits About A Black
Hole}
\centerline{Stephen C.\ Bell}
In a previous column,$^1$ we obtained numerical solutions of the
relativistic motion of a particle orbiting a spherical, nonrotating
central mass. In this supplement, we consider the more general case
of a rotating central mass. This case is important because it is
likely that the gravitational collapse that produces a black hole
would yield one with nonzero rotational angular momentum.$^2$
As before, we assume that the central mass (the black hole) is at
the origin of a rectangular perifocal coordinate system.$^1$ We also
assume that the axis of rotation of the central mass is oriented
along the
$z$-axis (see Fig.~1 of Ref.~1). Given a test particle's initial
position and velocity and a time of flight $t$, our goal is to find
its orbit. If the central mass does not rotate, the gravitational
field is spherically symmetrical, the test particle's motion is in
the orbital plane, and the perifocal frame is convenient because the
value of
$z$ can be assumed to be zero. In contrast, if the central mass does
rotate, its gravitational field is asymmetrical, and the test
particle does not remain in the orbital plane. However, we shall
see that the perifocal frame is still convenient.
As before, we assume that a stationary observer can measure the
position and velocity of the test particle, and that there is a
clock attached to the test particle.
We first transform the rectangular coordinates of the test particle
to polar coordinates using the relations $r=\sqrt{x^2+y^2+z^2}$,
$\phi=\tan^{-1}y/x$, and $\theta=\cos^{-1}z/r$. The four-dimensional
spacetime is defined by the square of the proper time differential
$d\tau$:
$$(d\tau)^2=(d\xv)^T\Gv\, d\xv, \eqno(1)$$ where $\xv
=(r,\theta,\phi,t)^T$ is the transpose of the column four vector.
The proper time $\tau$ is the time on the clock attached
to the test particle measured by the
stationary observer, and the coordinate time
$t$ is the time on the stationary clock at the observer's location.
The tensor
$\Gv$ can be expressed as a
$4\times 4$ matrix:
$$\Gv= \left[ {\matrix{g_{rr}&g_{r \theta}&g_{r \phi}&g_{r t}\cr
g_{\theta r}&g_{\theta \theta }&g_{\theta \phi}&g_{\theta t}\cr
g_{\phi r}&g_{\phi \theta}&g_{\phi \phi }&g_{\phi t}\cr
g_{t r}&g_{t \theta}&g_{t\phi }&g_{tt}\cr
}} \right].\eqno(2)$$
For non-Euclidean spacetimes,$^{2,3}$ the matrix elements of $\Gv$
are functions of $\xv$, but for simplicity, we
suppress this functional dependence. We also understand that
each component of $\xv$ is a function of $\tau$.
In the general theory of relativity, the orbital motion of a test
particle is given by the equations of a geodesic defined as
the path of shortest distance in a four-dimensional
spacetime.$^{2,3}$ The proper time equations of motion derived
directly from the geodesic equations can be written in polar
coordinates as
$${d^2r \over d\tau ^2}=g^{rr}\!\left[{1 \over 2}
{{d\xv^T \over d\tau }{\partial \Gv \over \partial r}}
{d\xv \over d\tau }-{dr \over d\tau }\
{{{\partial g_{rr} \over \partial \xv} \cdot {d\xv \over
d\tau} }} \right] \eqno(3)$$
$${d^2\theta \over d\tau ^2}=g^{\theta\theta} \! \left[{1 \over 2}
{{d\xv^T \over d\tau }{\partial \Gv \over \partial
\theta}}
{d\xv \over d\tau }-{d\theta \over d\tau }
{{{\partial g_{\theta\theta} \over \partial \xv} \cdot
{d\xv \over d\tau} }} \right] \eqno(4)$$
$${d^2\phi \over d\tau ^2}=g^{\phi\phi}\!\left[{1 \over 2}
{{d\xv^T \over d\tau }{\partial \Gv \over \partial \phi}}
{d\xv \over d\tau }-{d\phi \over d\tau }
{{{\partial g_{\phi\phi} \over \partial \xv} \cdot {d\xv
\over d\tau} }} \right] \eqno(5)$$
$${d^2t \over d\tau ^2}=g^{tt}\!\left[{1 \over 2}
{{d\xv^T \over d\tau }{\partial \Gv \over \partial t}}
{d\xv \over d\tau }-{dt \over d\tau }
{{{\partial g_{tt} \over \partial \xv} \cdot {d\xv \over d\tau} }}
\right].\eqno(6)$$
The elements $g^{rr}$, $g^{\theta\theta}$, etc. in Eqs.~(3)--(6) are
the matrix elements of the inverse of $\Gv$.
The proper time polar velocities in Eqs.~(3)--(6) are given by
$d\xv/d\tau = (dt/d\tau)(d\xv/dt)$, where the time dilation factor$^{1,3}$
$dt/d\tau$ is given by
$${dt \over d\tau }=\left( {{d\xv^T \over dt}
\Gv{d\xv \over dt}}
\right)^{-1/2}.\eqno(7)$$
To solve Eqs.~(3)--(7) we need to know the form of $\Gv$.
Schwarzschild derived the elements of $\Gv$ for a nonrotating
central mass. Kerr realized that the spacetime surrounding rotating
black holes would be different, and
derived the appropriate form of $\Gv$ for rotating black holes. His
main result is that the mass of a rotating black hole is
insufficient for describing its gravitational field, and that the
angular momentum also plays a role. The matrix elements of $\Gv$
for a rotating black hole are functions of $r$ and $\theta$, its
mass $M$, and its angular momentum $J$, and are given by:$^2$
$$\eqalign{&g_{rr}=-{\Sigma \over {c^2\Delta}},
\quad g_{\theta \theta }=-{\Sigma \over {c^2}},
\quad g_{\phi \phi }=-{1 \over c^2}
\left[ {{{\left( {r^2+a^2} \right)^2-\Delta \, a^2\sin^2 \theta}
\over \Sigma }} \right]\sin ^2\theta \cr
& \quad g_{tt}={{\Delta -a^2\sin ^2\theta} \over \Sigma},
\quad g_{t \phi}=g_{\phi t}={1 \over c^2}
\left[{{{a\sin^2\theta(r^2+a^2-\Delta )} \over \Sigma}}
\right],\cr}
\eqno(8)$$
where
$$\Sigma =r^2+a^2\cos ^2\theta \quad {\rm and}
\quad \Delta =r^2-{2GMr \over c^2}+a^2. \eqno(9)$$
As usual, $c$ is
the speed of light, and $G$ is Newton's gravitational constant. The
characteristic length associated
with the angular momentum is given by $a = J/(Mc)$. All other matrix
elements of $\Gv$ besides those given in Eq.~(8) are zero. Because
the matrix elements in Eq.~(8) are functions of $r$ and $\theta$
only, the terms involving
$\partial \Gv/\partial \phi$ and
$\partial \Gv/\partial t$ in Eqs.~(5) and (6) are zero,
and the derivatives with respect to $\phi$
and $t$ in $\partial g_{rr}/\partial \xv$,
$\partial g_{\theta \theta}/\partial \xv$,
$\partial g_{\phi \phi}/\partial \xv$ and
$\partial g_{tt}/\partial \xv$ are zero. Note that if the central
mass is not rotating, $J = 0$, and Eqs.~(3)--(8) reduce to the
Schwarzschild equations.$^{1,2,3}$
Equations (2)--(8) provide the complete set of equations needed to
simulate the Kerr orbits around rotating black holes. We use a double
precision version of the fourth-order Runge-Kutta subroutine {\tt
RK4} given in {\it Numerical Recipes}.$^4$ This subroutine requires
the subroutine {\tt DERIVS} which computes the proper time
accelerations given by Eqs.~(3)--(6). We suggest that the nonzero
partial derivatives in Eqs.~(3)--(6) be calculated analytically
using a symbolic manipulation program.
Suppose that a spacecraft (the test particle) is placed in orbit about a black
hole of mass $M = 10M_\odot$, where $M_\odot$ is the mass of the Sun.
Another convenient characteristic length, called the {\it
Schwarzschild radius}, is given by$^1$ $r_s = 2GM/c^2$. (For $M =
10M_\odot$,
$r_s = 29.533\, {\rm km}$.) From Eq.~(8), we see that $g_{rr}$ is
undefined for
$\Delta = 0$. The value of $r$ at which $\Delta=0$ is called the {\it Kerr
radius} $r_k$ and is given by $r_k = m +
\sqrt{m^2-a^2}$, where $m=r_s/2$. Because $r_k$ must be real and
positive, we see that
$a\le m$. Because $a$ must be proportional to
$m$, we set $a = pm$, where $0 \le p \le 1$. In the following, we
have set
$p =0.5$, half-way between its minimum and maximum values, which
implies a fairly rapidly rotating black hole. In the following, we
measure lengths in units of $m$, angles in radians, and time in
seconds.
We consider a spacecraft with an apogee point given by
$\xv = (25.0, \pi/2, 0.0, 0.0)^T$. The spacecraft at this point has
a velocity such that the instantaneous plane of its motion makes a
$45^\circ$ angle with the $xy$-plane, that is, the inclination of
the orbit at this point is $45^\circ$. The proper time velocity at
this point is
$d\xv^T\!/d\tau = (dr/d\tau, d\theta/d\tau, d\phi/d\tau,
dt/d\tau)^T$, with $dr/d\tau = 0$, $d\theta/ d\tau = -85.58610183$,
$d\phi/d\tau = 85.58610183$, and $dt/d\tau = 1.0539077376$. These
values for $\xv$ and $d\xv^T\!/d\tau$ correspond to an initial speed
of approximately $0.14\, {\rm c}$ and a Newtonian eccentricity of $e
= 0.5$. At this speed, the spacecraft completes one revolution about
the black hole in approximately $0.013\, {\rm s}$ of proper time. To
see the relativistic effects on the spacecraft's orbit, we consider
the spacecraft's motion for $0.04\, {\rm s}$ of proper time. The
proper time step should be no greater than
$10^{-5}\, {\rm s}$.
The view of the stationary observer for the first $0.04\, {\rm s}$
(proper time) of the spacecraft's motion is shown in Fig.~1. The
total amount of coordinate time that passed on the stationary clock
was $0.04593\, {\rm s}$. The difference is due to time dilation.
Note
the large precession of the major axes of the first and second
orbits caused by the extreme relativistic effects. Due to its
rotation, the black hole's gravitational field is asymmetrical
because
$g_{t\phi}$ and $g_{\phi t}$ are nonzero everywhere. This asymmetry
causes the spacecraft to continually come out of the plane producing
inclination shifts of its orbits. These inclination shifts are
better seen from the perspective of the stationary observer as shown
in Fig.~2, where the spacecraft has been propagated for an additional
time of $0.042\, {\rm s}$. If the spacecraft's motion is considered
for an extended time, its trajectory would form a torus around the
rotating black hole. The spacecraft could then explore a volume of
space surrounding the black hole without initiating any
orbit-changing maneuvers. In contrast, Newtonian mechanics implies
that the spacecraft would need to perform rocket engine burns to
change its inclination.
\noindent {\bf Suggestions for further study}
\noindent (1) Derive closed form solutions for the
partial derivatives needed in Eqs. (3)--(6).
\noindent (2) Suppose that a test particle's initial position is
$(25.0, 0.0, 0.0)$ in rectangular coordinates. Its velocity is
orthogonal to this position vector and makes a $45^\circ$ angle with
the $xy$-plane. The test particle's coordinate time speed is
$0.1\, {\rm c}$. The central mass is a black hole with a mass of
$M = 10M_\odot$, and the dimensionless angular momentum of $a =
0.5\,m$. Predict the orbit of the test particle.
\noindent (3) Suppose that
the spacecraft has completed its scientific mission, and needs to
escape from the black hole to travel to its next destination.
According to classical Newtonian mechanics, the procedure would be
for the spacecraft to do an engine burn at the appropriate
point on its last bound orbit (usually the perigee point). After the
burn, the spacecraft would be injected to a hyperbolic orbit whose
asymptote points in the desired escape direction.$^5$ The escape
orbit's eccentricity should be greater than unity, so that the
spacecraft would no longer be bound to the black hole.
Show that, due to relativistic effects, the burn would actually
not have to immediately produce an eccentricity greater than unity
for the spacecraft to eventually escape. Produce a three-dimensional
plot showing the pathway of the spacecraft before and after the
escape-producing engine burn. (Assume that the engine
burn produces an instantaneous increase in the velocity of the
spacecraft at the perigee point of one the orbits depicted in Fig.~1.)
\noindent (4) Some researchers$^{2,3,6}$ have stated
that no stable circular orbits exist around Schwarzs\-child
(nonrotating) black holes for $3m \le r_c \le 6m$, where
$r_c$ is the radius of the circular orbit. That is, a stationary
observer cannot see a test particle continually circle a
Schwarzschild black hole at a radius in this interval. Use the Kerr
equations with $a=0$ and show that this implication is incorrect.
Hint: choose an initial state vector whose Newtonian eccentricity is
zero. Explain why all continually circular orbits
that lie above the event horizon are stable.
\medskip \noindent {\bf Acknowledgments}
\noindent I thank the editors for their suggestions and comments
during the writing of this manu\-script. A Fortran program that
calculates Kerr orbits around rotating black holes is available for
downloading from CiP.
\medskip
%\vfill\eject\noindent {\bf References}
\noindent {\bf References}
\frenchspacing
\item{1.} S. C. Bell, ``A numerical solution of the relativistic
Kepler problem, {\sl Computers in Physics} {\bf 9}, 281 (1995).
\item{2.} R. M. Wald, {\sl General Relativity}, University of
Chicago Press (1984).
\item{3.} D. F. Lawden, {\sl An Introduction to Tensor Calculus,
Relativity and Cosmology}, third edition, John Wiley \& Sons, New
York (1982).
\item{4.} W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P.
Flannery, {\sl Numerical Recipes}, second edition, Cambridge
University Press (1992).
\item{5.} S. C. Bell, P. P. Rao, and M. A. Ginsburg, ``Monte Carlo
Analysis of the Titan/Transfer Orbit Stage Guidance System for a
Planetary Mission,'' {\sl Journal of Guidance, Control, and
Dynamics} {\bf 18}, 121 (1995).
\item{6.} S. Chandrasekhar, {\sl The Mathematical Theory of Black
Holes}, Oxford University Press, New York (1983).
\smallskip \parindent=0pt
{\bf Figure Captions}
\smallskip
Fig.~1. Three-dimensional plot of the spacecraft's orbital path.
Distances are measured in units of $m = 14.767$ km. The
precession at the beginning of the second orbit is
$276.7^\circ$.
Fig.~2. View of the spacecraft's path showing the shifts
of the inclination of its orbit due to the rotation of the black
hole.
\bye