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}\) -