# Geodesics in Kerr spacetime¶

## Computation with kerrgeodesic_gw¶

This Jupyter notebook requires SageMath (version $\geq$ 8.2), with the package kerrgeodesic_gw (version $\geq$ 0.3). To install the latter, simply run

sage -pip install kerrgeodesic_gw

in a terminal.

In [1]:
version()

Out[1]:
'SageMath version 9.4, Release Date: 2021-08-22'

First, we set up the notebook to use LaTeX-formatted display:

In [2]:
%display latex


and we ask for CPU demanding computations to be performed in parallel on 8 processes:

In [3]:
Parallelism().set(nproc=8)


A Kerr black bole is entirely defined by two parameters $(m, a)$, where $m$ is the black hole mass and $a$ is the black hole angular momentum divided by $m$. In this notebook, we shall set $m=1$ and we denote the angular momentum parameter $a$ by the symbolic variable a, using a0 for a specific numerical value:

In [4]:
a = var('a')
a0 = 0.998


The spacetime object is created as an instance of the class KerrBH:

In [5]:
from kerrgeodesic_gw import KerrBH
M = KerrBH(a)
print(M)

Kerr spacetime M


The object M is endowed with many methods, which can be discovered via the TAB key:

In [6]:
# M.<TAB>


One of them returns the Boyer-Lindquist coordinate $r$ of the event horizon:

In [7]:
rH = M.event_horizon_radius()
rH

Out[7]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\sqrt{-a^{2} + 1} + 1$
In [8]:
rH0 = rH.subs({a: a0})
rH0

Out[8]:
$\newcommand{\Bold}[1]{\mathbf{#1}}1.06321392251712$

Another one returns the chart of Boyer-Lindquist coordinates and allows the user to instanciate the Python variables (t, r, th, ph) to the coordinates $(t,r,\theta,\phi)$:

In [9]:
BL.<t, r, th, ph> = M.boyer_lindquist_coordinates()
BL

Out[9]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\left(M,(t, r, {\theta}, {\phi})\right)$

The metric tensor is naturally returned by the method metric():

In [10]:
g = M.metric()
g.display()

Out[10]:
$\newcommand{\Bold}[1]{\mathbf{#1}}g = \left( -\frac{a^{2} \cos\left({\theta}\right)^{2} + r^{2} - 2 \, r}{a^{2} \cos\left({\theta}\right)^{2} + r^{2}} \right) \mathrm{d} t\otimes \mathrm{d} t + \left( -\frac{2 \, a r \sin\left({\theta}\right)^{2}}{a^{2} \cos\left({\theta}\right)^{2} + r^{2}} \right) \mathrm{d} t\otimes \mathrm{d} {\phi} + \left( \frac{a^{2} \cos\left({\theta}\right)^{2} + r^{2}}{a^{2} + r^{2} - 2 \, r} \right) \mathrm{d} r\otimes \mathrm{d} r + \left( a^{2} \cos\left({\theta}\right)^{2} + r^{2} \right) \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + \left( -\frac{2 \, a r \sin\left({\theta}\right)^{2}}{a^{2} \cos\left({\theta}\right)^{2} + r^{2}} \right) \mathrm{d} {\phi}\otimes \mathrm{d} t + \left( \frac{2 \, a^{2} r \sin\left({\theta}\right)^{4} + {\left(a^{2} r^{2} + r^{4} + {\left(a^{4} + a^{2} r^{2}\right)} \cos\left({\theta}\right)^{2}\right)} \sin\left({\theta}\right)^{2}}{a^{2} \cos\left({\theta}\right)^{2} + r^{2}} \right) \mathrm{d} {\phi}\otimes \mathrm{d} {\phi}$

Again many methods are available on metric objects and one can discover them via the TAB key:

In [11]:
# g.<TAB>


For instance christoffel_symbols_display() computes the Christofell symbols with respect to the default chart (Boyer-Lindquist) and displays them:

In [12]:
g.christoffel_symbols_display()

Out[12]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\begin{array}{lcl} \Gamma_{ \phantom{\, t} \, t \, r }^{ \, t \phantom{\, t} \phantom{\, r} } & = & -\frac{a^{4} - r^{4} - {\left(a^{4} + a^{2} r^{2}\right)} \sin\left({\theta}\right)^{2}}{a^{2} r^{4} + r^{6} - 2 \, r^{5} + {\left(a^{6} + a^{4} r^{2} - 2 \, a^{4} r\right)} \cos\left({\theta}\right)^{4} + 2 \, {\left(a^{4} r^{2} + a^{2} r^{4} - 2 \, a^{2} r^{3}\right)} \cos\left({\theta}\right)^{2}} \\ \Gamma_{ \phantom{\, t} \, t \, {\theta} }^{ \, t \phantom{\, t} \phantom{\, {\theta}} } & = & -\frac{2 \, a^{2} r \cos\left({\theta}\right) \sin\left({\theta}\right)}{a^{4} \cos\left({\theta}\right)^{4} + 2 \, a^{2} r^{2} \cos\left({\theta}\right)^{2} + r^{4}} \\ \Gamma_{ \phantom{\, t} \, r \, {\phi} }^{ \, t \phantom{\, r} \phantom{\, {\phi}} } & = & -\frac{{\left(a^{3} r^{2} + 3 \, a r^{4} - {\left(a^{5} - a^{3} r^{2}\right)} \cos\left({\theta}\right)^{2}\right)} \sin\left({\theta}\right)^{2}}{a^{2} r^{4} + r^{6} - 2 \, r^{5} + {\left(a^{6} + a^{4} r^{2} - 2 \, a^{4} r\right)} \cos\left({\theta}\right)^{4} + 2 \, {\left(a^{4} r^{2} + a^{2} r^{4} - 2 \, a^{2} r^{3}\right)} \cos\left({\theta}\right)^{2}} \\ \Gamma_{ \phantom{\, t} \, {\theta} \, {\phi} }^{ \, t \phantom{\, {\theta}} \phantom{\, {\phi}} } & = & -\frac{2 \, {\left(a^{5} r \cos\left({\theta}\right) \sin\left({\theta}\right)^{5} - {\left(a^{5} r + a^{3} r^{3}\right)} \cos\left({\theta}\right) \sin\left({\theta}\right)^{3}\right)}}{a^{6} \cos\left({\theta}\right)^{6} + 3 \, a^{4} r^{2} \cos\left({\theta}\right)^{4} + 3 \, a^{2} r^{4} \cos\left({\theta}\right)^{2} + r^{6}} \\ \Gamma_{ \phantom{\, r} \, t \, t }^{ \, r \phantom{\, t} \phantom{\, t} } & = & \frac{a^{2} r^{2} + r^{4} - 2 \, r^{3} - {\left(a^{4} + a^{2} r^{2} - 2 \, a^{2} r\right)} \cos\left({\theta}\right)^{2}}{a^{6} \cos\left({\theta}\right)^{6} + 3 \, a^{4} r^{2} \cos\left({\theta}\right)^{4} + 3 \, a^{2} r^{4} \cos\left({\theta}\right)^{2} + r^{6}} \\ \Gamma_{ \phantom{\, r} \, t \, {\phi} }^{ \, r \phantom{\, t} \phantom{\, {\phi}} } & = & -\frac{{\left(a^{3} r^{2} + a r^{4} - 2 \, a r^{3} - {\left(a^{5} + a^{3} r^{2} - 2 \, a^{3} r\right)} \cos\left({\theta}\right)^{2}\right)} \sin\left({\theta}\right)^{2}}{a^{6} \cos\left({\theta}\right)^{6} + 3 \, a^{4} r^{2} \cos\left({\theta}\right)^{4} + 3 \, a^{2} r^{4} \cos\left({\theta}\right)^{2} + r^{6}} \\ \Gamma_{ \phantom{\, r} \, r \, r }^{ \, r \phantom{\, r} \phantom{\, r} } & = & \frac{{\left(a^{2} r - a^{2}\right)} \sin\left({\theta}\right)^{2} + a^{2} - r^{2}}{a^{2} r^{2} + r^{4} - 2 \, r^{3} + {\left(a^{4} + a^{2} r^{2} - 2 \, a^{2} r\right)} \cos\left({\theta}\right)^{2}} \\ \Gamma_{ \phantom{\, r} \, r \, {\theta} }^{ \, r \phantom{\, r} \phantom{\, {\theta}} } & = & -\frac{a^{2} \cos\left({\theta}\right) \sin\left({\theta}\right)}{a^{2} \cos\left({\theta}\right)^{2} + r^{2}} \\ \Gamma_{ \phantom{\, r} \, {\theta} \, {\theta} }^{ \, r \phantom{\, {\theta}} \phantom{\, {\theta}} } & = & -\frac{a^{2} r + r^{3} - 2 \, r^{2}}{a^{2} \cos\left({\theta}\right)^{2} + r^{2}} \\ \Gamma_{ \phantom{\, r} \, {\phi} \, {\phi} }^{ \, r \phantom{\, {\phi}} \phantom{\, {\phi}} } & = & \frac{{\left(a^{4} r^{2} + a^{2} r^{4} - 2 \, a^{2} r^{3} - {\left(a^{6} + a^{4} r^{2} - 2 \, a^{4} r\right)} \cos\left({\theta}\right)^{2}\right)} \sin\left({\theta}\right)^{4} - {\left(a^{2} r^{5} + r^{7} - 2 \, r^{6} + {\left(a^{6} r + a^{4} r^{3} - 2 \, a^{4} r^{2}\right)} \cos\left({\theta}\right)^{4} + 2 \, {\left(a^{4} r^{3} + a^{2} r^{5} - 2 \, a^{2} r^{4}\right)} \cos\left({\theta}\right)^{2}\right)} \sin\left({\theta}\right)^{2}}{a^{6} \cos\left({\theta}\right)^{6} + 3 \, a^{4} r^{2} \cos\left({\theta}\right)^{4} + 3 \, a^{2} r^{4} \cos\left({\theta}\right)^{2} + r^{6}} \\ \Gamma_{ \phantom{\, {\theta}} \, t \, t }^{ \, {\theta} \phantom{\, t} \phantom{\, t} } & = & -\frac{2 \, a^{2} r \cos\left({\theta}\right) \sin\left({\theta}\right)}{a^{6} \cos\left({\theta}\right)^{6} + 3 \, a^{4} r^{2} \cos\left({\theta}\right)^{4} + 3 \, a^{2} r^{4} \cos\left({\theta}\right)^{2} + r^{6}} \\ \Gamma_{ \phantom{\, {\theta}} \, t \, {\phi} }^{ \, {\theta} \phantom{\, t} \phantom{\, {\phi}} } & = & \frac{2 \, {\left(a^{3} r + a r^{3}\right)} \cos\left({\theta}\right) \sin\left({\theta}\right)}{a^{6} \cos\left({\theta}\right)^{6} + 3 \, a^{4} r^{2} \cos\left({\theta}\right)^{4} + 3 \, a^{2} r^{4} \cos\left({\theta}\right)^{2} + r^{6}} \\ \Gamma_{ \phantom{\, {\theta}} \, r \, r }^{ \, {\theta} \phantom{\, r} \phantom{\, r} } & = & \frac{a^{2} \cos\left({\theta}\right) \sin\left({\theta}\right)}{a^{2} r^{2} + r^{4} - 2 \, r^{3} + {\left(a^{4} + a^{2} r^{2} - 2 \, a^{2} r\right)} \cos\left({\theta}\right)^{2}} \\ \Gamma_{ \phantom{\, {\theta}} \, r \, {\theta} }^{ \, {\theta} \phantom{\, r} \phantom{\, {\theta}} } & = & \frac{r}{a^{2} \cos\left({\theta}\right)^{2} + r^{2}} \\ \Gamma_{ \phantom{\, {\theta}} \, {\theta} \, {\theta} }^{ \, {\theta} \phantom{\, {\theta}} \phantom{\, {\theta}} } & = & -\frac{a^{2} \cos\left({\theta}\right) \sin\left({\theta}\right)}{a^{2} \cos\left({\theta}\right)^{2} + r^{2}} \\ \Gamma_{ \phantom{\, {\theta}} \, {\phi} \, {\phi} }^{ \, {\theta} \phantom{\, {\phi}} \phantom{\, {\phi}} } & = & -\frac{{\left({\left(a^{6} + a^{4} r^{2} - 2 \, a^{4} r\right)} \cos\left({\theta}\right)^{5} + 2 \, {\left(a^{4} r^{2} + a^{2} r^{4} - 2 \, a^{2} r^{3}\right)} \cos\left({\theta}\right)^{3} + {\left(a^{2} r^{4} + r^{6} + 2 \, a^{4} r + 4 \, a^{2} r^{3}\right)} \cos\left({\theta}\right)\right)} \sin\left({\theta}\right)}{a^{6} \cos\left({\theta}\right)^{6} + 3 \, a^{4} r^{2} \cos\left({\theta}\right)^{4} + 3 \, a^{2} r^{4} \cos\left({\theta}\right)^{2} + r^{6}} \\ \Gamma_{ \phantom{\, {\phi}} \, t \, r }^{ \, {\phi} \phantom{\, t} \phantom{\, r} } & = & -\frac{a^{3} \cos\left({\theta}\right)^{2} - a r^{2}}{a^{2} r^{4} + r^{6} - 2 \, r^{5} + {\left(a^{6} + a^{4} r^{2} - 2 \, a^{4} r\right)} \cos\left({\theta}\right)^{4} + 2 \, {\left(a^{4} r^{2} + a^{2} r^{4} - 2 \, a^{2} r^{3}\right)} \cos\left({\theta}\right)^{2}} \\ \Gamma_{ \phantom{\, {\phi}} \, t \, {\theta} }^{ \, {\phi} \phantom{\, t} \phantom{\, {\theta}} } & = & -\frac{2 \, a r \cos\left({\theta}\right)}{{\left(a^{4} \cos\left({\theta}\right)^{4} + 2 \, a^{2} r^{2} \cos\left({\theta}\right)^{2} + r^{4}\right)} \sin\left({\theta}\right)} \\ \Gamma_{ \phantom{\, {\phi}} \, r \, {\phi} }^{ \, {\phi} \phantom{\, r} \phantom{\, {\phi}} } & = & \frac{r^{5} + {\left(a^{4} r - a^{4}\right)} \cos\left({\theta}\right)^{4} - a^{2} r^{2} - 2 \, r^{4} + {\left(2 \, a^{2} r^{3} + a^{4} - a^{2} r^{2}\right)} \cos\left({\theta}\right)^{2}}{a^{2} r^{4} + r^{6} - 2 \, r^{5} + {\left(a^{6} + a^{4} r^{2} - 2 \, a^{4} r\right)} \cos\left({\theta}\right)^{4} + 2 \, {\left(a^{4} r^{2} + a^{2} r^{4} - 2 \, a^{2} r^{3}\right)} \cos\left({\theta}\right)^{2}} \\ \Gamma_{ \phantom{\, {\phi}} \, {\theta} \, {\phi} }^{ \, {\phi} \phantom{\, {\theta}} \phantom{\, {\phi}} } & = & \frac{a^{4} \cos\left({\theta}\right)^{5} + 2 \, {\left(a^{2} r^{2} - a^{2} r\right)} \cos\left({\theta}\right)^{3} + {\left(r^{4} + 2 \, a^{2} r\right)} \cos\left({\theta}\right)}{{\left(a^{4} \cos\left({\theta}\right)^{4} + 2 \, a^{2} r^{2} \cos\left({\theta}\right)^{2} + r^{4}\right)} \sin\left({\theta}\right)} \end{array}$

The Riemann curvature tensor is naturally returned by the method riemann()

In [13]:
Riem = g.riemann()
print(Riem)

Tensor field Riem(g) of type (1,3) on the Kerr spacetime M


The $R^t_{rtr}$ and $R^t_{r\theta\phi}$ components:

In [14]:
Riem[0,1,0,1]

Out[14]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\frac{3 \, a^{4} r \cos\left({\theta}\right)^{4} + 3 \, a^{2} r^{3} + 2 \, r^{5} - {\left(9 \, a^{4} r + 7 \, a^{2} r^{3}\right)} \cos\left({\theta}\right)^{2}}{a^{2} r^{6} + r^{8} - 2 \, r^{7} + {\left(a^{8} + a^{6} r^{2} - 2 \, a^{6} r\right)} \cos\left({\theta}\right)^{6} + 3 \, {\left(a^{6} r^{2} + a^{4} r^{4} - 2 \, a^{4} r^{3}\right)} \cos\left({\theta}\right)^{4} + 3 \, {\left(a^{4} r^{4} + a^{2} r^{6} - 2 \, a^{2} r^{5}\right)} \cos\left({\theta}\right)^{2}}$
In [15]:
Riem[0,1,2,3]

Out[15]:
$\newcommand{\Bold}[1]{\mathbf{#1}}-\frac{{\left({\left(a^{7} + a^{5} r^{2} - 2 \, a^{5} r\right)} \cos\left({\theta}\right)^{5} - {\left(3 \, a^{7} + 8 \, a^{5} r^{2} + 5 \, a^{3} r^{4} - 2 \, a^{5} r - 6 \, a^{3} r^{3}\right)} \cos\left({\theta}\right)^{3} + 3 \, {\left(3 \, a^{5} r^{2} + 5 \, a^{3} r^{4} + 2 \, a r^{6} - 2 \, a^{3} r^{3}\right)} \cos\left({\theta}\right)\right)} \sin\left({\theta}\right)}{a^{2} r^{6} + r^{8} - 2 \, r^{7} + {\left(a^{8} + a^{6} r^{2} - 2 \, a^{6} r\right)} \cos\left({\theta}\right)^{6} + 3 \, {\left(a^{6} r^{2} + a^{4} r^{4} - 2 \, a^{4} r^{3}\right)} \cos\left({\theta}\right)^{4} + 3 \, {\left(a^{4} r^{4} + a^{2} r^{6} - 2 \, a^{2} r^{5}\right)} \cos\left({\theta}\right)^{2}}$

Of course, since the Kerr metric is a solution of the vacuum Einstein equation, the Ricci tensor identically vanishes:

In [16]:
g.ricci().display()

Out[16]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\mathrm{Ric}\left(g\right) = 0$

## Bound timelike geodesic¶

Let us choose the initial point $P$ for the geodesic:

In [17]:
P = M.point((0, 6, pi/2, 0), name='P')
print(P)

Point P on the Kerr spacetime M


A geodesic is constructed by providing the range $[\lambda_{\rm min},\lambda_{\rm max}]$ of the affine parameter $\lambda$, the initial point and either

• (i) the Boyer-Lindquist components $(p^t_0, p^r_0, p^\theta_0, p^\phi_0)$ of the initial 4-momentum vector $p_0 = \left. \frac{\mathrm{d}x}{\mathrm{d}\lambda}\right| _{\lambda_{\rm min}}$,
• (ii) the four integral of motions $(\mu, E, L, Q)$
• or (iii) some of the components of $p_0$ along with with some integrals of motion. We shall also specify some numerical value for the Kerr spin parameter $a$.

Examples of (i) and (iii) are provided below. Here, we choose $\lambda\in[0, 300\, m]$, the option (ii) and $a=0.998 \,m$, where $m$ in the black hole mass::

In [18]:
Li = M.geodesic([0, 300], P, mu=1, E=0.883, L=1.982, Q=0.467, a_num=0.998,
name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector:

$\newcommand{\Bold}[1]{\mathbf{#1}}p = 1.29225788954106 \frac{\partial}{\partial t } + 0.00438084990626460 \frac{\partial}{\partial r } + 0.0189826106258554 \frac{\partial}{\partial {\theta} } + 0.0646134478134985 \frac{\partial}{\partial {\phi} }$
The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

In [19]:
print(Li)

Geodesic Li of the Kerr spacetime M


The numerical integration of the geodesic equation is performed via integrate(), by providing the integration step $\delta\lambda$ in units of $m$:

In [20]:
Li.integrate(step=0.005)


We can then plot the geodesic:

In [21]:
Li.plot()

Out[21]:

Actually, many options can be passed to the method plot(). For instance to a get a 3D spacetime diagram:

In [22]:
Li.plot(coordinates='txy')

Out[22]:

or to get the trace of the geodesic in the $(x,y)$ plane:

In [23]:
Li.plot(coordinates='xy', plot_points=2000)

Out[23]:

or else to get the trace in the $(x,z)$ plane:

In [24]:
Li.plot(coordinates='xz')

Out[24]:

As a curve, the geodesic $\mathcal{L}$ is a map from an interval of $\mathbb{R}$ to the spacetime $M$:

In [25]:
Li.display()

Out[25]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\begin{array}{llcl} \mathcal{L}:& \left(0, 300\right) & \longrightarrow & M \end{array}$
In [26]:
Li.domain()

Out[26]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0, 300\right)$
In [27]:
Li.codomain()

Out[27]:
$\newcommand{\Bold}[1]{\mathbf{#1}}M$

It maps values of $\lambda$ to spacetime points:

In [28]:
Li(0)

Out[28]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\mbox{Point on the Kerr spacetime M}$
In [29]:
Li(0).coordinates()  # coordinates in the default chart (BL)

Out[29]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0.0, 6.0, 1.5707963267948966, 0.0\right)$
In [30]:
BL(Li(0))  # equivalent to above

Out[30]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\left(0.0, 6.0, 1.5707963267948966, 0.0\right)$
In [31]:
Li(300).coordinates()

Out[31]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\left(553.4637326813786, 3.703552505462962, 1.6613834863942039, 84.62814710987239\right)$

The initial 4-momentum vector $p_0$ is returned by the method initial_tangent_vector():

In [32]:
p0 = Li.initial_tangent_vector()
print(p0)

Tangent vector p at Point P on the Kerr spacetime M

In [33]:
p0 in M.tangent_space(P)

Out[33]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}$
In [34]:
p0.display()

Out[34]:
$\newcommand{\Bold}[1]{\mathbf{#1}}p = 1.29225788954106 \frac{\partial}{\partial t } + 0.00438084990626460 \frac{\partial}{\partial r } + 0.0189826106258554 \frac{\partial}{\partial {\theta} } + 0.0646134478134985 \frac{\partial}{\partial {\phi} }$
In [35]:
p0[:]

Out[35]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\left[1.29225788954106, 0.00438084990626460, 0.0189826106258554, 0.0646134478134985\right]$

For instance, the components $p^t_0$ and $p^\phi_0$ are recovered by

In [36]:
p0[0], p0[3]

Out[36]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\left(1.29225788954106, 0.0646134478134985\right)$

Let us check that the scalar square of $p_0$ is $-1$, i.e. is consistent with the mass parameter $\mu = 1$ used in the construction of the geodesic:

In [37]:
g = M.metric()
g.at(P)(p0, p0).subs(a=0.998)

Out[37]:
$\newcommand{\Bold}[1]{\mathbf{#1}}-1.00000000000000$

The 4-momentum vector $p$ at any value of the affine parameter $\lambda$ is obtained by the method evaluate_tangent_vector(); for instance, for $\lambda=200\,m$:

In [38]:
p = Li.evaluate_tangent_vector(200)
p.display()

Out[38]:
$\newcommand{\Bold}[1]{\mathbf{#1}}1.316592599498746 \frac{\partial}{\partial t } -0.07370434215844164 \frac{\partial}{\partial r } -0.01091195426423706 \frac{\partial}{\partial {\theta} } + 0.07600209768075264 \frac{\partial}{\partial {\phi} }$
In [39]:
p in M.tangent_space(Li(200))

Out[39]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}$

The particle mass $\mu$ computed at a given value of $\lambda$ is returned by the method evaluate_mu():

In [40]:
Li.evaluate_mu(0)

Out[40]:
$\newcommand{\Bold}[1]{\mathbf{#1}}1.00000000000000$

Of course, it should be conserved along $\mathcal{L}$; actually it is, up to the numerical accuracy::

In [41]:
Li.evaluate_mu(300)

Out[41]:
$\newcommand{\Bold}[1]{\mathbf{#1}}1.0000117978600134$

Similarly, the conserved energy $E$, conserved angular momentum $L$ and Carter constant $Q$ are computed at any value of $\lambda$ by respectively evaluate_E(), evaluate_L() and evaluate_Q():

In [42]:
Li.evaluate_E(0)

Out[42]:
$\newcommand{\Bold}[1]{\mathbf{#1}}0.883000000000000$
In [43]:
Li.evaluate_L(0)

Out[43]:
$\newcommand{\Bold}[1]{\mathbf{#1}}1.98200000000000$
In [44]:
Li.evaluate_Q(0)

Out[44]:
$\newcommand{\Bold}[1]{\mathbf{#1}}0.467000000000000$

Let us check that the values of $\mu$, $E$, $L$ and $Q$ evaluated at $\lambda=300 \, m$ are equal to those at $\lambda=0$ up to the numerical accuracy of the integration scheme:

In [45]:
Li.check_integrals_of_motion(300)

Out[45]:
 quantity value initial value diff. relative diff. $\mu^2$ $1.0000235958592163$ $1.00000000000000$ $0.00002360$ $0.00002360$ $E$ $0.883067996080701$ $0.883000000000000$ $0.00006800$ $0.00007701$ $L$ $1.98248080818931$ $1.98200000000000$ $0.0004808$ $0.0002426$ $Q$ $0.467214137649741$ $0.467000000000000$ $0.0002141$ $0.0004585$

Decreasing the integration step leads to smaller errors:

In [46]:
Li.integrate(step=0.001)
Li.check_integrals_of_motion(300)

Out[46]:
 quantity value initial value diff. relative diff. $\mu^2$ $1.0000047183936422$ $1.00000000000000$ $4.718 \times 10^{-6}$ $4.718 \times 10^{-6}$ $E$ $0.883013604456676$ $0.883000000000000$ $0.00001360$ $0.00001541$ $L$ $1.98209626120918$ $1.98200000000000$ $0.00009626$ $0.00004857$ $Q$ $0.467042771975860$ $0.467000000000000$ $0.00004277$ $0.00009159$

## Various ways to initialize a geodesic¶

Instead of providing the integral of motions, as for Li above, one can initialize a geodesic by providing the Boyer-Lindquist components $(p^t_0, p^r_0, p^\theta_0, p^\phi_0)$ of the initial 4-momentum vector $p_0$. For instance:

In [47]:
print(p0)

Tangent vector p at Point P on the Kerr spacetime M

In [48]:
p0[:]

Out[48]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\left[1.29225788954106, 0.00438084990626460, 0.0189826106258554, 0.0646134478134985\right]$
In [49]:
Li2 = M.geodesic([0, 300], P, pt0=p0[0], pr0=p0[1], pth0=p0[2], pph0=p0[3],
a_num=0.998)
Li2.initial_tangent_vector() == p0

Out[49]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}$

As a check, we recover the same values of $(\mu, E, L, Q)$ as those that were used to initialize Li:

In [50]:
Li2.evaluate_mu(0)

Out[50]:
$\newcommand{\Bold}[1]{\mathbf{#1}}1.00000000000000$
In [51]:
Li2.evaluate_E(0)

Out[51]:
$\newcommand{\Bold}[1]{\mathbf{#1}}0.883000000000000$
In [52]:
Li2.evaluate_L(0)

Out[52]:
$\newcommand{\Bold}[1]{\mathbf{#1}}1.98200000000000$
In [53]:
Li2.evaluate_Q(0)

Out[53]:
$\newcommand{\Bold}[1]{\mathbf{#1}}0.467000000000000$

We may also initialize a geodesic by providing the mass $\mu$ and the three spatial components $(p^r_0, p^\theta_0, p^\phi_0)$ of the initial 4-momentum vector:

In [54]:
Li3 = M.geodesic([0, 300], P, mu=1, pr0=p0[1], pth0=p0[2], pph0=p0[3],
a_num=0.998)


The component $p^t_0$ is then automatically computed:

In [55]:
Li3.initial_tangent_vector()[:]

Out[55]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\left[1.29225788954106, 0.00438084990626460, 0.0189826106258554, 0.0646134478134985\right]$

and we check the identity with the initial vector of Li, up to numerical errors:

In [56]:
(Li3.initial_tangent_vector() - p0)[:]

Out[56]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\left[2.22044604925031 \times 10^{-16}, 0, 0, 0\right]$

Another way to initialize a geodesic is to provide the conserved energy $E$, the conserved angular momentum $L$ and the two components $(p^r_0, p^\theta_0)$ of the initial 4-momentum vector:

In [57]:
Li4 = M.geodesic([0, 300], P, E=0.8830, L=1.982, pr0=p0[1], pth0=p0[2],
a_num=0.998)
Li4.initial_tangent_vector()[:]

Out[57]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\left[1.29225788954106, 0.00438084990626460, 0.0189826106258554, 0.0646134478134985\right]$

Again, we get a geodesic equivalent to Li:

In [58]:
(Li4.initial_tangent_vector() - p0)[:]

Out[58]:
$\newcommand{\Bold}[1]{\mathbf{#1}}\left[0, 0, 0, 0\right]$

## Ingoing null geodesic with negative angular momentum¶

We choose a ingoing null geodesic in the equatorial plane with $L = -6 E < 0$, starting at the point of Boyer-Lindquist coordinates $(t,r,\theta,\phi) = (0, 12, \pi/2, 0)$:

In [59]:
lambda_max = 13.063
Li = M.geodesic([0, lambda_max], M((0,12,pi/2,0)), mu=0, E=1, L=-6, Q=0,
r_increase=False, a_num=a0, name='Li', latex_name=r'\mathcal{L}',
verbose=True)

Initial tangent vector:

$\newcommand{\Bold}[1]{\mathbf{#1}}p = 1.20797381595070 \frac{\partial}{\partial t } -0.901996260873419 \frac{\partial}{\partial r } -0.0399489777089388 \frac{\partial}{\partial {\phi} }$
The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

In [60]:
Li.integrate(step=0.00002)


We plot the trajectory of the geodesic in the equatorial plane with the tangent vector at 6 points:

In [61]:
Li.plot(coordinates='xy', color='green', prange=[0, lambda_max], plot_points=40000,
display_tangent=True, scale=0.002, color_tangent='brown',
plot_points_tangent=6, axes_labels=[r'$x/m$', r'$y/m$'])

Out[61]:

Note the turning point in $\phi$ and the final winding in the direction of the black rotation (counterclockwise in the figure).

In [62]:
Li.check_integrals_of_motion(0.99999*lambda_max)

Out[62]:
 quantity value initial value diff. relative diff. $\mu^2$ $0.000211815000511706$ $0.000000000000000$ $0.0002118$ - $E$ $1.00000186057173$ $1.00000000000000$ $1.861 \times 10^{-6}$ $1.861 \times 10^{-6}$ $L$ $-5.99999596865291$ $-6.00000000000000$ $4.031 \times 10^{-6}$ $-6.719 \times 10^{-7}$ $Q$ $1.31244559310830 \times 10^{-31}$ $0$ $1.312 \times 10^{-31}$ -