Spherical 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.

In [1]:
version()
Out[1]:
'SageMath version 9.2, Release Date: 2020-10-24'

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.95

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 Boyer-Lindquist coordinate $r$ of the event horizon:

In [6]:
rH = M.event_horizon_radius()
rH
Out[6]:
In [7]:
rH0 = rH.subs({a: a0})
show(LatexExpr(r'r_+ = '), rH0)
In [8]:
show(LatexExpr(r'r_- = '),
     M.inner_horizon_radius().subs({a: a0}))

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)$:

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

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

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

Spherical photon orbits

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

In [11]:
r = var('r')
lsph(a, r) = (r^2*(3 - r) - a^2*(r + 1))/(a*(r -1))
lsph
Out[11]:
In [12]:
qsph(a, r) = r^3 / (a^2*(r - 1)^2) * (4*a^2 - r*(r - 3)^2)
qsph
Out[12]:

$\theta$-turning points:

In [13]:
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
Out[13]:
In [14]:
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
Out[14]:

Spherical photon orbit at $r_0 = 3 m$ ($q = q_{\rm max} = 27 m^2$)

In [15]:
r0 = 3.
E = 1
L = lsph(a0, r0)
Q = qsph(a0, r0)
L, Q
Out[15]:
In [16]:
theta0(a0, L, Q)
Out[16]:
In [17]:
P = M.point((0, r0, 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.
In [18]:
lmax = 100 # lambda_max
In [19]:
Li = M.geodesic([0, lmax], P, mu=0, E=E, L=L, Q=Q, a_num=a0,
                name='Li', latex_name=r'\mathcal{L}', verbose=True)
Initial tangent vector: 
The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].
In [20]:
v0 = Li.initial_tangent_vector()
Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3],
                a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True)
Initial tangent vector: 
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$:

In [21]:
Li.integrate(step=0.001, method='dopri5')
Li.check_integrals_of_motion(0.999*lmax)
Out[21]:
quantity value initial value diff. relative diff.
-
In [22]:
print("Final point: ")
show(BL[:], "=", BL(Li(0.999*lmax)))
Final point: 
In [23]:
lplot = 0.32*lmax
print("max lambda (plot): ", lplot)
Li.plot(prange=(0, lplot), plot_points=2000, thickness=3, color='green',
        display_tangent=True, scale=0.15, width_tangent=2, color_tangent='green', 
        plot_points_tangent=20, horizon_color='lightgrey') \
     + line([(0,0,-3), (0,0,3)], color='black', thickness=2)
max lambda (plot):  32.0000000000000
Out[23]:
In [24]:
graph = Li.plot(coordinates='xy', prange=(0, lplot), plot_points=2000, 
                thickness=3, color='green', display_tangent=True, scale=0.15, 
                width_tangent=2, color_tangent='green', plot_points_tangent=20, 
                horizon_color='lightgrey', axes_labels=[r'$x/m$', r'$y/m$'])
graph.save("gik_spher_3d_r_30_xy.pdf")
graph
Out[24]:
In [25]:
Li.plot(coordinates='xz', prange=(0, lplot), plot_points=2000, 
        thickness=3, color='green', display_tangent=True, scale=0.15, 
        width_tangent=2, color_tangent='green', plot_points_tangent=20, 
        horizon_color='lightgrey', axes_labels=[r'$x/m$', r'$z/m$'])
Out[25]:

Prograde spherical photon orbit at $r_0 = 1.6m$

In [26]:
r0 = 1.6
L = lsph(a0, r0)
Q = qsph(a0, r0)
L, Q
Out[26]:
In [27]:
theta0(a0, L, Q)
Out[27]:
In [28]:
P = M.point((0, r0, pi/2, 0), name='P')
lmax = 70
Li = M.geodesic([0, lmax], P, mu=0, E=E, L=L, Q=Q, a_num=a0,
                name='Li', latex_name=r'\mathcal{L}', verbose=True)
Initial tangent vector: 
The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].
In [29]:
v0 = Li.initial_tangent_vector()
Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3],
                a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True)
Initial tangent vector: 
The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].
In [30]:
Li.integrate(step=0.0004, method='dopri5')
Li.check_integrals_of_motion(0.999*lmax)
Out[30]:
quantity value initial value diff. relative diff.
In [31]:
print("Final point: ")
show(BL[:], "=", BL(Li(0.999*lmax)))
Final point: 
In [32]:
lplot = 0.11*lmax
print("max lambda (plot): ", lplot)
Li.plot(prange=(0, lplot), plot_points=2000, thickness=3, color='green',
        display_tangent=True, scale=0.03, width_tangent=1, color_tangent='green', 
        plot_points_tangent=20, horizon_color='lightgrey') \
    + line([(0,0,-1.5), (0,0,1.5)], color='black', thickness=2)
max lambda (plot):  7.70000000000000
Out[32]:
In [33]:
Li.plot(prange=(0, lmax), plot_points=2000, thickness=3, color='green',
        display_tangent=True, scale=0.03, width_tangent=1, color_tangent='green', 
        plot_points_tangent=40, horizon_color='lightgrey') \
    + line([(0,0,-1.5), (0,0,1.5)], color='black', thickness=2)
Out[33]:

Spherical photon orbit at $r_0=2.8m$

In [34]:
r0 = 2.8
L = lsph(a0, r0)
Q = qsph(a0, r0)
L, Q
Out[34]:
In [35]:
theta0(a0, L, Q)
Out[35]:
In [36]:
P = M.point((0, r0, pi/2, 0), name='P')
lmax = 100
Li = M.geodesic([0, lmax], P, mu=0, E=E, L=L, Q=Q, a_num=a0,
                name='Li', latex_name=r'\mathcal{L}', verbose=True)
Initial tangent vector: 
The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].
In [37]:
v0 = Li.initial_tangent_vector()
Li = M.geodesic([0, lmax], P, pt0=v0[0], pr0=0, pth0=v0[2], pph0=v0[3],
                a_num=a0, name='Li', latex_name=r'\mathcal{L}', verbose=True)
Initial tangent vector: 
The curve was correctly set.
Parameters appearing in the differential system defining the curve are [a].
In [38]:
Li.integrate(step=0.001, method='dopri5')
Li.check_integrals_of_motion(0.999*lmax)
Out[38]:
quantity value initial value diff. relative diff.
In [39]:
print("Final point: ")
show(BL[:], "=", BL(Li(0.999*lmax)))
Final point: 
In [40]:
lplot = 0.38*lmax
print("max lambda (plot): ", lplot)
Li.plot(prange=(0, lplot), plot_points=2000, thickness=3, color='green',
        display_tangent=True, scale=0.12, width_tangent=2, color_tangent='green', 
        plot_points_tangent=20, horizon_color='lightgrey')
max lambda (plot):  38.0000000000000
Out[40]: