Kerr_null_geod_plots.ipynb

Kernel: SageMath 9.2

Plots of null geodesics in Kerr spacetime

Computation with kerrgeodesic_gw

This Jupyter/SageMath notebook is relative to the lectures Geometry and physics of black holes.

It requires SageMath (version $\geq$ 8.2), with the package kerrgeodesic_gw (version $\geq$ 0.3.2). To install the latter, simply run

sage -pip install kerrgeodesic_gw

in a terminal.

version()

'SageMath version 9.2, Release Date: 2020-10-24'

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

%display latex

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

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:

a = var('a')
a0 = 0.95

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

from kerrgeodesic_gw import KerrBH
M = KerrBH(a)
print(M)

Kerr spacetime M

The Boyer-Lindquist coordinate $r$ of the event horizon:

rH = M.event_horizon_radius()
rH

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\sqrt{-a^{2} + 1} + 1$$

rH0 = rH.subs({a: a0})
show(LatexExpr(r'r_+ = '), rH0)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_+ = 1.31224989991992$$

show(LatexExpr(r'r_- = '),
     M.inner_horizon_radius().subs({a: a0}))

$$\renewcommand{\Bold}[1]{\mathbf{#1}}r_- = 0.687750100080080$$

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

BL.<t, r, th, ph> = M.boyer_lindquist_coordinates()
BL

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

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

g = M.metric()
g.display()

$$\renewcommand{\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}$$

Functions $\ell(r_0)$ and $q(r_0)$ for spherical photon orbits:

r = var('r')
lsph(a, r) = (r^2*(3 - r) - a^2*(r + 1))/(a*(r -1))
lsph

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( a, r \right) \ {\mapsto} \ -\frac{a^{2} {\left(r + 1\right)} + {\left(r - 3\right)} r^{2}}{a {\left(r - 1\right)}}$$

qsph(a, r) = r^3 / (a^2*(r - 1)^2) * (4*a^2 - r*(r - 3)^2)
qsph

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( a, r \right) \ {\mapsto} \ -\frac{{\left({\left(r - 3\right)}^{2} r - 4 \, a^{2}\right)} r^{3}}{a^{2} {\left(r - 1\right)}^{2}}$$

$\theta$-turning points:

theta0(a, l, q) = arccos(sqrt(1/2*(1 - (l^2+q)/a^2 + sqrt((1 - (l^2+q)/a^2)^2 + 4*q/a^2))))
theta0

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( a, l, q \right) \ {\mapsto} \ \arccos\left(\sqrt{\frac{1}{2} \, \sqrt{{\left(\frac{l^{2} + q}{a^{2}} - 1\right)}^{2} + \frac{4 \, q}{a^{2}}} - \frac{l^{2} + q}{2 \, a^{2}} + \frac{1}{2}}\right)$$

theta1(a, l, q) = arccos(sqrt(1/2*(1 - (l^2+q)/a^2 - sqrt((1 - (l^2+q)/a^2)^2 + 4*q/a^2))))
theta1

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( a, l, q \right) \ {\mapsto} \ \arccos\left(\sqrt{-\frac{1}{2} \, \sqrt{{\left(\frac{l^{2} + q}{a^{2}} - 1\right)}^{2} + \frac{4 \, q}{a^{2}}} - \frac{l^{2} + q}{2 \, a^{2}} + \frac{1}{2}}\right)$$

Plot of $\mathcal{R}(r)$ for various values of $(\ell, q)$

l = var('l', latex_name=r'\ell')
q = var('q')
R(a, l, q, r) = r^4 + (a^2 - l^2 - q)*r^2 + 2*(q + (l - a)^2)*r - a^2*q
R

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( a, {\ell}, q, r \right) \ {\mapsto} \ r^{4} - a^{2} q + {\left(a^{2} - {\ell}^{2} - q\right)} r^{2} + 2 \, {\left({\left(a - {\ell}\right)}^{2} + q\right)} r$$

xmin = 1.2
xmax = 3.7
ymin = -10
ymax = 40

L = 0.5
Q = 16
g0 = plot(R(a0, L, Q, r), (r, xmin, xmax), color='red', thickness=2,         
          axes_labels=[r'$r/m$', r'$\mathcal{R}(r)/m^4$'],
          frame=True, gridlines=True)
g0 += line([(rH0, ymin), (rH0, ymax)], color='black', thickness=2)
g0 += line([(xmin, -0.4), (1.1*xmax, -0.4)], color='lightgreen', alpha=0.5, thickness=3)
g0 += polygon([(0.9*xmin, ymin), (rH0, ymin), (rH0, ymax), (0.9*xmin, ymax)], 
              color='lightgrey')
g0 += text(r'$\ell = 0.5\, m$', (1.8, 35), fontsize=16, color='red',
           background_color='white')
g0 += text(r'$q = 16\, m^2$', (1.8, 28.5), fontsize=16, color='red',
           background_color='white')
g0.set_axes_range(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
g0.set_aspect_ratio(0.04)
g0.save('gik_R_in_M1_1.pdf', figsize=5)
show(g0, figsize=5)

r0 = 2.2
L0 = lsph(a0, r0)
Q0 = qsph(a0, r0)
L0, Q0

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(0.863157894736842, 18.0416251154201\right)$$

g1 = plot(R(a0, L0, Q0, r), (r, rH0, xmax), color='red', thickness=2,         
          axes_labels=[r'$r/m$', None],
          frame=True, gridlines=True)
g1 += line([(rH0, ymin), (rH0, ymax)], color='black', thickness=2)
g1 += line([(xmin, -0.4), (1.1*xmax, -0.4)], color='lightgreen', thickness=3)
g1 += point((r0, 0), size=60, color='red', zorder=100)
g1 += text(r'$r_0$', (r0, -2.8), color='red', fontsize=16)
g1 += text(r'$\ell = {:.3}\, m$'.format(float(L0)), (1.9, 35), fontsize=16, 
           color='red', background_color='white')
g1 += text(r'$q = {:.4}\, m^2$'.format(float(Q0)), (1.92, 28.5), fontsize=16, 
           color='red', background_color='white')
g1 += polygon([(0.9*xmin, ymin), (rH0, ymin), (rH0, ymax), (0.9*xmin, ymax)], 
              color='lightgrey')
g1.set_aspect_ratio(0.04)
g1.set_axes_range(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
g1.save('gik_R_in_M1_2.pdf', figsize=5)
show(g1, figsize=5)

L = 1
Q = 20
r1 = find_root(R(a0, L, Q, r), rH0, r0)
r2 = find_root(R(a0, L, Q, r), r0, 4)
g2 = plot(R(a0, L, Q, r), (r, xmin, xmax), color='red', thickness=2,         
          axes_labels=[r'$r/m$', None],
          frame=True, gridlines=True)
g2 += line([(rH0, ymin), (rH0, ymax)], color='black', thickness=2)
g2 += line([(xmin, -0.4), (r1, -0.4)], color='lightgreen', thickness=3)
g2 += line([(r2, -0.4), (1.1*xmax, -0.4)], color='lightgreen', thickness=3)
g2 += point((r1, 0), size=60, color='red', zorder=100)
g2 += text(r'$r_1$', (r1, -2.8), color='red', fontsize=16)
g2 += point((r2, 0), size=60, color='red', zorder=100)
g2 += text(r'$r_2$', (1.02*r2, -2.8), color='red', fontsize=16)
g2 += text(r'$\ell = m$', (1.7, 35), fontsize=16, color='red',
           background_color='white')
g2 += text(r'$q = 20\, m^2$', (1.85, 28.5), fontsize=16, color='red',
           background_color='white')
g2 += polygon([(0.9*xmin, ymin), (rH0, ymin), (rH0, ymax), (0.9*xmin, ymax)], 
              color='lightgrey')
g2.set_aspect_ratio(0.04)
g2.set_axes_range(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
g2.save('gik_R_in_M1_3.pdf', figsize=5)
show(g2, figsize=5)

L = 2*rH0/a0
Q = 8
r2 = find_root(R(a0, L, Q, r), r0, 8)
L, Q, r2

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(2.76263136825246, 8, 2.722847233324403\right)$$

g3 = plot(R(a0, L, Q, r), (r, xmin, xmax), color='red', thickness=2,         
          axes_labels=[r'$r/m$', None],
          frame=True, gridlines=True)
g3 += line([(rH0, ymin), (rH0, ymax)], color='black', thickness=2)
g3 += line([(xmin, -0.4), (rH0, -0.4)], color='lightgreen', thickness=3)
g3 += line([(r2, -0.4), (1.1*xmax, -0.4)], color='lightgreen', thickness=3)
g3 += point((r2, 0), size=60, color='red', zorder=100)
g3 += text(r'$r_0$', (1.02*r2, -2.8), color='red', fontsize=16)
g3 += text(r'$\ell = 2mr_+/a$', (2, 35), fontsize=16, 
           color='red', background_color='white')
g3 += text(r'$q = 8 \, m^2$', (1.8, 28), fontsize=16, 
           color='red', background_color='white')
g3 += polygon([(0.9*xmin, ymin), (rH0, ymin), (rH0, ymax), (0.9*xmin, ymax)], 
              color='lightgrey')
g3.set_aspect_ratio(0.04)
g3.set_axes_range(xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax)
g3.save('gik_R_in_M1_4.pdf', figsize=5)
show(g3, figsize=5)

Critical null geodesic

Parameters of a spherical photon orbit at $r_0 = 2.2m$:

r0 = 2.2
L0 = lsph(a0, r0)
Q0 = qsph(a0, r0)
L0, Q0

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(0.863157894736842, 18.0416251154201\right)$$

Chosen parameters for the geodesic:

L = L0
Q = Q0
L, Q

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(0.863157894736842, 18.0416251154201\right)$$

theta0(a0, L, Q)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}0.195946218453965$$

r_init = 40
P = M.point((0, r_init, 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.

lmax = 80 # lambda_max

Li = M.geodesic([0, lmax], P, mu=0, E=1, L=L, Q=Q, r_increase=False,
                a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 1.05260305969646 \frac{\partial}{\partial t } -0.994675862738240 \frac{\partial}{\partial r } + 0.00265471463828243 \frac{\partial}{\partial {\theta} } + 0.000570385018106026 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

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

Li.integrate(step=0.0005, method='dopri5')
#Li.integrate(step=0.001)
Li.check_integrals_of_motion(0.999*lmax)

print("Final point: ")
show(BL[:], "=", BL(Li(0.999*lmax)))

Final point: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(t, r, {\theta}, {\phi}\right) \verb|=| \left(206.66703023467483, 2.2000863112656326, 1.9511762444244498, 53.094380043176905\right)$$

lplot = lmax
print("max lambda (plot): ", lplot)
graph1 = Li.plot(prange=(30, lplot), plot_points=2000, thickness=4, color='green',
                 display_tangent=True, scale=0.1, width_tangent=10, color_tangent='green', 
                 plot_points_tangent=12, horizon_color='lightgrey') \
         + line([(0,0,-3), (0,0,3)], color='black', thickness=2)
graph1

max lambda (plot):  80

Escaping null geodesic close to the critical one

L = 1.0001*L0
Q = Q0
L, Q

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(0.863244210526316, 18.0416251154201\right)$$

theta0(a0, L, Q)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}0.195965348586193$$

r_init = 40
P = M.point((0, r_init, pi/2, 0), name='P')
lmax = 80
Li = M.geodesic([0, lmax], P, mu=0, E=1, L=L, Q=Q, r_increase=False,
                a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 1.05260305700070 \frac{\partial}{\partial t } -0.994675815686174 \frac{\partial}{\partial r } + 0.00265471463828243 \frac{\partial}{\partial {\theta} } + 0.000570438933462204 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

Li.integrate(step=0.001, method='dopri5')
Li.check_integrals_of_motion(0.999*lmax)

print("Final point: ")
show(BL[:], "=", BL(Li(0.999*lmax)))

Final point: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(t, r, {\theta}, {\phi}\right) \verb|=| \left(142.73493888601888, 22.257873149906203, 0.49847494769804584, 23.336351331857166\right)$$

lplot = 0.8*lmax
print("max lambda (plot): ", lplot)
graph2 = Li.plot(prange=(30, lplot), plot_points=2000, thickness=4, color='green',
                display_tangent=True, scale=0.2, width_tangent=10, color_tangent='green', 
                plot_points_tangent=12, horizon_color='lightgrey') \
         + line([(0,0,-3), (0,0,3)], color='black', thickness=2)
graph2

max lambda (plot):  64.0000000000000

Trapped null geodesic close to the critical one

L = 0.9999*L0
Q = Q0
L, Q

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(0.863071578947368, 18.0416251154201\right)$$

theta0(a0, L, Q)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}0.195927088185978$$

r_init = 40
P = M.point((0, r_init, pi/2, 0), name='P')
lmax = 56.81
Li = M.geodesic([0, lmax], P, mu=0, E=1, L=L, Q=Q, r_increase=False,
                a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True)

Initial tangent vector: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}p = 1.05260306239223 \frac{\partial}{\partial t } -0.994675909785857 \frac{\partial}{\partial r } + 0.00265471463828243 \frac{\partial}{\partial {\theta} } + 0.000570331102749847 \frac{\partial}{\partial {\phi} }$$

The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].

Li.integrate(step=0.0005, method='dopri5')
Li.check_integrals_of_motion(0.999*lmax)

/home/eric/sage/9.2/local/lib/python3.8/site-packages/scipy/integrate/_ode.py:1179: UserWarning: dopri5: larger nsteps is needed
  warnings.warn('{:s}: {:s}'.format(self.__class__.__name__,

print("Final point: ")
show(BL[:], "=", BL(Li(0.999*lmax)))

Final point: 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(t, r, {\theta}, {\phi}\right) \verb|=| \left(128.57773967099268, 1.3375219985463933, 1.267137358002562, 28.267159612791097\right)$$

lplot = lmax
print("max lambda (plot): ", lplot)
graph3 = Li.plot(prange=(30, lplot), plot_points=2000, thickness=4, color='orange',
                 plot_horizon=False) 
graph2 + graph3

max lambda (plot):  56.8100000000000