Plot of principal null geodesics in Kerr spacetime

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

The computations make use of tools developed through the SageManifolds project.

In [1]:
version()
Out[1]:
'SageMath version 9.3.rc2, Release Date: 2021-04-06'

First we set up the notebook to display mathematical objects using LaTeX rendering:

In [2]:
%display latex

Spacetime manifold

We declare the Kerr spacetime as a 4-dimensional Lorentzian manifold $M$:

In [3]:
M = Manifold(4, 'M', structure='Lorentzian')
print(M)
4-dimensional Lorentzian manifold M

and introduce the (3+1 version of) Kerr coordinates $(\tilde{t},r,\theta,\tilde{\varphi})$ as a chart KC on $M$, via the method chart(). The argument of the latter is a string (delimited by r"..." because of the backslash symbols) expressing the coordinates names, their ranges (the default is $(-\infty,+\infty)$) and their LaTeX symbols:

In [4]:
KC.<tt,r,th,tph> = M.chart(r"tt:\tilde{t} r th:(0,pi):\theta tph:(0,2*pi):periodic:\tilde{\varphi}") 
print(KC); KC
Chart (M, (tt, r, th, tph))
Out[4]:
In [5]:
a = 0.9
m = 1 
rp = m + sqrt(m^2-a^2)
rm = m - sqrt(m^2-a^2)
(rp,rm)
Out[5]:

Plot of the principal null geodesics in the $(\tilde{t}, r)$ plane

In [6]:
tt_in(r, v) = v - r
tt_in
Out[6]:
In [7]:
tt_out(r, u) = u + r + 2*m/sqrt(m^2 - a^2)*(rp*ln(abs((r - rp)/(2*m))) - rm*ln(abs((r - rm)/(2*m))))
tt_out
Out[7]:
In [8]:
rmin = -8; rmax = 8
graph = Graphics()
for u0 in range(-20, 20, 2):
    graph += plot(tt_out(r, u0), (r, rmin, rmax), color='green', ticks=2,
                  axes_labels=[r"$r/m$", r"$\tilde{t}/m$"])
In [9]:
for v0 in range(-20, 20, 2):
    graph += plot(tt_in(r, v0), (r, rmin, rmax), color='green', linestyle='--')
In [10]:
H = line(((rp, -8), (rp, 8)), color='black', thickness=3) + \
    text(r'$\mathscr{H}$', (rp+0.5, 4.7), color='black', fontsize=20)
Hin = line(((rm, -8), (rm, 8)), color='peru', thickness=3) + \
      text(r'$\mathscr{H}_{\rm in}$', (rm-0.6, 4.7), color='peru', fontsize=20)
graph += H + Hin
show(graph, aspect_ratio=1, ymin=-5, ymax=5, figsize=8)

Adding the vectors $k$ and $\ell$ to the plot

In [11]:
k = M.vector_field(1, -1, 0, 0, name='k')
k.display()
Out[11]:
In [12]:
el = M.vector_field(1/2 + m*r/(r^2 + a^2),
                    1/2 - m*r/(r^2 + a^2),
                    0,
                    a/(r^2 + a^2),
                    name='el', latex_name=r'\ell')
el.display()
Out[12]:

We add the vectors $k$ and $\ell$ at the intersection of the $v=6m$ ingoing geodesic with the $u=-6m$ outgoing one:

In [13]:
u0, v0 = -6, 6
r0 = RDF(find_root(tt_in(r, v0) == tt_out(r, u0), 2, 6))
tt0 = tt_in(r0, v0)
tt0, r0
Out[13]:
In [14]:
p0 = M((tt0, r0, pi/2, 0), name='p_0')
k0 = k.at(p0)
print(k0)
Tangent vector k at Point p_0 on the 4-dimensional Lorentzian manifold M
In [15]:
el0 = el.at(p0)
print(el0)
Tangent vector el at Point p_0 on the 4-dimensional Lorentzian manifold M
In [16]:
graph +=  k0.plot(chart=KC, ambient_coords=(r, tt), color='green', 
                  scale=1.5, fontsize=18, label_offset=0.3)  \
          + el0.plot(chart=KC, ambient_coords=(r, tt), color='green', 
                     parameters={m: 1}, scale=1.5, fontsize=18,
                     label_offset=0.25)
graph.save("ker_princ_null_geod.pdf", aspect_ratio=1, ymin=-5, ymax=5,
           figsize=8)
show(graph, aspect_ratio=1, ymin=-5, ymax=5, figsize=8)

Plot of $u - \tilde{t}$ and $\tilde{\tilde{\varphi}} - \tilde{\varphi}$ as functions of $r$

In [17]:
U(r) = u - tt_out(r, u)
U(r)
Out[17]:
In [18]:
graph = plot(U(r), (r, -8, 8), color='red', thickness=2, 
             axes_labels=[r"$r/m$", r"$(u - \tilde{t})/m$"], frame=True, 
             gridlines=True)
graph += line(((rp, -10), (rp, 10)), color='black', thickness=1.5, linestyle='--') + \
         line(((rm, -10), (rm, 10)), color='peru', thickness=1.5, linestyle='--')
graph.save('ker_u_r.pdf', aspect_ratio=1, ymin=-8, ymax=8)
show(graph, aspect_ratio=1, ymin=-8, ymax=8)
In [19]:
ttphi(r) = - a/sqrt(m^2 - a^2)*ln(abs((r - rp)/(r - rm)))
ttphi(r)
Out[19]:
In [20]:
graph = plot(ttphi(r), (r, -8, 8), color='red', thickness=2, 
             axes_labels=[r"$r/m$", r"$\tilde{\tilde{\varphi}} - \tilde{\varphi}$"], 
             frame=True, gridlines=True)
graph += line(((rp, -10), (rp, 10)), color='black', thickness=1.5, linestyle='--') + \
         line(((rm, -10), (rm, 10)), color='peru', thickness=1.5, linestyle='--')
graph.save('ker_ttphi_r.pdf', aspect_ratio=1, ymin=-8, ymax=8)
show(graph, aspect_ratio=1, ymin=-8, ymax=8)
In [ ]: