Kerr_extremal.ipynb

Kernel: SageMath 9.3.rc5

Extremal Kerr black hole

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.

version()

'SageMath version 9.3.rc5, Release Date: 2021-04-30'

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

%display latex

To speed up computations, we ask for running them in parallel on 8 threads:

Parallelism().set(nproc=8)

Spacetime manifold

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

M = Manifold(4, 'M', structure='Lorentzian')
print(M)

4-dimensional Lorentzian manifold M

We then introduce (3+1 version of) the 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:

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

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

KC.coord_range()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}{\tilde{t}} :\ \left( -\infty, +\infty \right) ;\quad r :\ \left( -\infty, +\infty \right) ;\quad {\theta} :\ \left( 0 , \pi \right) ;\quad {\tilde{\varphi}} :\ \left[ 0 , 2 \, \pi \right] \mbox{(periodic)}$$

Metric tensor

The mass parameter $m$ of the extremal Kerr spacetime is declared as a symbolic variable:

m = var('m', domain='real')
assume(m>0)

We get the (yet undefined) spacetime metric:

g = M.metric()

and initialize it by providing its components in the coordinate frame associated with the Kerr coordinates, which is the current manifold's default frame:

rho2 = r^2 + (m*cos(th))^2
g[0,0] = - (1 - 2*m*r/rho2)
g[0,1] = 2*m*r/rho2
g[0,3] = -2*m^2*r*sin(th)^2/rho2
g[1,1] = 1 + 2*m*r/rho2
g[1,3] = -m*(1 + 2*m*r/rho2)*sin(th)^2
g[2,2] = rho2
g[3,3] = (r^2 + m^2 + 2*m^3*r*sin(th)^2/rho2)*sin(th)^2
g.display()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}g = \left( \frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} - 1 \right) \mathrm{d} {\tilde{t}}\otimes \mathrm{d} {\tilde{t}} + \left( \frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} \right) \mathrm{d} {\tilde{t}}\otimes \mathrm{d} r + \left( -\frac{2 \, m^{2} r \sin\left({\theta}\right)^{2}}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} \right) \mathrm{d} {\tilde{t}}\otimes \mathrm{d} {\tilde{\varphi}} + \left( \frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} \right) \mathrm{d} r\otimes \mathrm{d} {\tilde{t}} + \left( \frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} + 1 \right) \mathrm{d} r\otimes \mathrm{d} r -m {\left(\frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} + 1\right)} \sin\left({\theta}\right)^{2} \mathrm{d} r\otimes \mathrm{d} {\tilde{\varphi}} + \left( m^{2} \cos\left({\theta}\right)^{2} + r^{2} \right) \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + \left( -\frac{2 \, m^{2} r \sin\left({\theta}\right)^{2}}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} \right) \mathrm{d} {\tilde{\varphi}}\otimes \mathrm{d} {\tilde{t}} -m {\left(\frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} + 1\right)} \sin\left({\theta}\right)^{2} \mathrm{d} {\tilde{\varphi}}\otimes \mathrm{d} r + {\left(\frac{2 \, m^{3} r \sin\left({\theta}\right)^{2}}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} + m^{2} + r^{2}\right)} \sin\left({\theta}\right)^{2} \mathrm{d} {\tilde{\varphi}}\otimes \mathrm{d} {\tilde{\varphi}}$$

A matrix view of the components with respect to the manifold's default vector frame:

g[:]

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrr} \frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} - 1 & \frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} & 0 & -\frac{2 \, m^{2} r \sin\left({\theta}\right)^{2}}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} \\ \frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} & \frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} + 1 & 0 & -m {\left(\frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} + 1\right)} \sin\left({\theta}\right)^{2} \\ 0 & 0 & m^{2} \cos\left({\theta}\right)^{2} + r^{2} & 0 \\ -\frac{2 \, m^{2} r \sin\left({\theta}\right)^{2}}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} & -m {\left(\frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} + 1\right)} \sin\left({\theta}\right)^{2} & 0 & {\left(\frac{2 \, m^{3} r \sin\left({\theta}\right)^{2}}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} + m^{2} + r^{2}\right)} \sin\left({\theta}\right)^{2} \end{array}\right)$$

The list of the non-vanishing components:

g.display_comp()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\begin{array}{lcl} g_{ \, {\tilde{t}} \, {\tilde{t}} }^{ \phantom{\, {\tilde{t}}}\phantom{\, {\tilde{t}}} } & = & \frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} - 1 \\ g_{ \, {\tilde{t}} \, r }^{ \phantom{\, {\tilde{t}}}\phantom{\, r} } & = & \frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} \\ g_{ \, {\tilde{t}} \, {\tilde{\varphi}} }^{ \phantom{\, {\tilde{t}}}\phantom{\, {\tilde{\varphi}}} } & = & -\frac{2 \, m^{2} r \sin\left({\theta}\right)^{2}}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} \\ g_{ \, r \, {\tilde{t}} }^{ \phantom{\, r}\phantom{\, {\tilde{t}}} } & = & \frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} \\ g_{ \, r \, r }^{ \phantom{\, r}\phantom{\, r} } & = & \frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} + 1 \\ g_{ \, r \, {\tilde{\varphi}} }^{ \phantom{\, r}\phantom{\, {\tilde{\varphi}}} } & = & -m {\left(\frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} + 1\right)} \sin\left({\theta}\right)^{2} \\ g_{ \, {\theta} \, {\theta} }^{ \phantom{\, {\theta}}\phantom{\, {\theta}} } & = & m^{2} \cos\left({\theta}\right)^{2} + r^{2} \\ g_{ \, {\tilde{\varphi}} \, {\tilde{t}} }^{ \phantom{\, {\tilde{\varphi}}}\phantom{\, {\tilde{t}}} } & = & -\frac{2 \, m^{2} r \sin\left({\theta}\right)^{2}}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} \\ g_{ \, {\tilde{\varphi}} \, r }^{ \phantom{\, {\tilde{\varphi}}}\phantom{\, r} } & = & -m {\left(\frac{2 \, m r}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} + 1\right)} \sin\left({\theta}\right)^{2} \\ g_{ \, {\tilde{\varphi}} \, {\tilde{\varphi}} }^{ \phantom{\, {\tilde{\varphi}}}\phantom{\, {\tilde{\varphi}}} } & = & {\left(\frac{2 \, m^{3} r \sin\left({\theta}\right)^{2}}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} + m^{2} + r^{2}\right)} \sin\left({\theta}\right)^{2} \end{array}$$

Einstein equation

Let us check that we are dealing with a solution of the vacuum Einstein equation:

g.ricci().display()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\mathrm{Ric}\left(g\right) = 0$$

Boyer-Lindquist coordinates

Let us introduce on the chart of Boyer-Lindquist coordinates $(t,r,\theta,\varphi)$:

BL.<t,r,th,ph> = M.chart(r"t r th:(0,pi):\theta ph:(0,2*pi):periodic:\varphi") 
print(BL); BL

Chart (M, (t, r, th, ph))

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

BL.coord_range()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}t :\ \left( -\infty, +\infty \right) ;\quad r :\ \left( -\infty, +\infty \right) ;\quad {\theta} :\ \left( 0 , \pi \right) ;\quad {\varphi} :\ \left[ 0 , 2 \, \pi \right] \mbox{(periodic)}$$

KC_to_BL = KC.transition_map(BL, [tt + 2*m^2/(r-m) - 2*m*ln(abs(r-m)/m),
                                  r,
                                  th,
                                  tph + m/(r-m)])
KC_to_BL.display()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left\{\begin{array}{lcl} t & = & -2 \, m \log\left(\frac{{\left| -m + r \right|}}{m}\right) - \frac{2 \, m^{2}}{m - r} + {\tilde{t}} \\ r & = & r \\ {\theta} & = & {\theta} \\ {\varphi} & = & {\tilde{\varphi}} - \frac{m}{m - r} \end{array}\right.$$

KC_to_BL.inverse().display()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left\{\begin{array}{lcl} {\tilde{t}} & = & -\frac{2 \, m^{2} \log\left(m\right) - 2 \, m r \log\left(m\right) - 2 \, m^{2} - {\left(m - r\right)} t - 2 \, {\left(m^{2} - m r\right)} \log\left({\left| -m + r \right|}\right)}{m - r} \\ r & = & r \\ {\theta} & = & {\theta} \\ {\tilde{\varphi}} & = & \frac{m {\varphi} - {\varphi} r + m}{m - r} \end{array}\right.$$

g.display(BL)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}g = \left( -\frac{m^{2} \cos\left({\theta}\right)^{2} - 2 \, m r + r^{2}}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} \right) \mathrm{d} t\otimes \mathrm{d} t + \left( -\frac{2 \, m^{2} r \sin\left({\theta}\right)^{2}}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} \right) \mathrm{d} t\otimes \mathrm{d} {\varphi} + \left( \frac{m^{2} \cos\left({\theta}\right)^{2} + r^{2}}{m^{2} - 2 \, m r + r^{2}} \right) \mathrm{d} r\otimes \mathrm{d} r + \left( m^{2} \cos\left({\theta}\right)^{2} + r^{2} \right) \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + \left( -\frac{2 \, m^{2} r \sin\left({\theta}\right)^{2}}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} \right) \mathrm{d} {\varphi}\otimes \mathrm{d} t + \left( \frac{2 \, m^{3} r \sin\left({\theta}\right)^{4} + {\left(m^{2} r^{2} + r^{4} + {\left(m^{4} + m^{2} r^{2}\right)} \cos\left({\theta}\right)^{2}\right)} \sin\left({\theta}\right)^{2}}{m^{2} \cos\left({\theta}\right)^{2} + r^{2}} \right) \mathrm{d} {\varphi}\otimes \mathrm{d} {\varphi}$$

Plot of the hypersurfaces $t=\mathrm{const}$ in terms of the Kerr coordinates

The plot is performed via the method plot of the Boyer-Lindquist chart:

graph = BL.plot(KC, ambient_coords=(r, tt), fixed_coords={th: pi/2, ph: 0}, 
                ranges={t: (-10, 10), r: (-10, 10)}, steps={t: 2, r: 2}, 
                plot_points=200, style={t: ':', r: '-'}, color={t: 'white', r: 'blue'},
                parameters={m: 1})
Hor = line(((1, -8), (1, 8)), color='black', thickness=3) + \
      text(r'$\mathscr{H}$', (1.5, 4.7), color='black', fontsize=20)
graph += Hor
graph.set_aspect_ratio(1)
graph.save("exk_BL_slicing.pdf", axes_labels=[r"$r/m$", r"$\tilde{t}/m$"], 
           ymin=-5, ymax=5, figsize=7)
graph.show(axes_labels=[r"$r/m$", r"$\tilde{t}/m$"], ymin=-5, ymax=5, figsize=7)

Ingoing principal null geodesics

k = M.vector_field(1, -1, 0, 0, name='k')
k.display()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}k = \frac{\partial}{\partial {\tilde{t}} }-\frac{\partial}{\partial r }$$

Let us check that $k$ is a null vector:

g(k, k).expr()

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

Check of $\nabla_k k = 0$:

nabla = g.connection()

nabla(k).contract(k).display()

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

Expression of $k$ with respect to the Boyer-Lindquist frame:

k.display(BL)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}k = \left( \frac{m^{2} + r^{2}}{m^{2} - 2 \, m r + r^{2}} \right) \frac{\partial}{\partial t } -\frac{\partial}{\partial r } + \left( \frac{m}{m^{2} - 2 \, m r + r^{2}} \right) \frac{\partial}{\partial {\varphi} }$$

Outgoing principal null geodesics

el = M.vector_field(1/2 + m*r/(r^2 + m^2),
                    1/2 - m*r/(r^2 + m^2),
                    0,
                    m/(r^2 + m^2),
                    name='el', latex_name=r'\ell')
el.display()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\ell = \left( \frac{m r}{m^{2} + r^{2}} + \frac{1}{2} \right) \frac{\partial}{\partial {\tilde{t}} } + \left( -\frac{m r}{m^{2} + r^{2}} + \frac{1}{2} \right) \frac{\partial}{\partial r } + \left( \frac{m}{m^{2} + r^{2}} \right) \frac{\partial}{\partial {\tilde{\varphi}} }$$

Let us check that $\ell$ is a null vector:

g(el, el).expr()

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

Expression of $\ell$ with respect to the Boyer-Lindquist frame:

el.display(BL)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\ell = \frac{1}{2} \frac{\partial}{\partial t } + \left( -\frac{m r}{m^{2} + r^{2}} + \frac{1}{2} \right) \frac{\partial}{\partial r } + \frac{m}{2 \, {\left(m^{2} + r^{2}\right)}} \frac{\partial}{\partial {\varphi} }$$

Computation of $\nabla_\ell \ell$:

acc = nabla(el).contract(el)
acc.display()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( -\frac{m^{5} + 2 \, m^{4} r - 2 \, m^{2} r^{3} - m r^{4}}{2 \, {\left(m^{6} + 3 \, m^{4} r^{2} + 3 \, m^{2} r^{4} + r^{6}\right)}} \right) \frac{\partial}{\partial {\tilde{t}} } + \left( -\frac{m^{5} - 2 \, m^{4} r + 2 \, m^{2} r^{3} - m r^{4}}{2 \, {\left(m^{6} + 3 \, m^{4} r^{2} + 3 \, m^{2} r^{4} + r^{6}\right)}} \right) \frac{\partial}{\partial r } + \left( -\frac{m^{4} - m^{2} r^{2}}{m^{6} + 3 \, m^{4} r^{2} + 3 \, m^{2} r^{4} + r^{6}} \right) \frac{\partial}{\partial {\tilde{\varphi}} }$$

We check that $\nabla_\ell \ell = \kappa \ell$:

kappa = acc[0] / el[0]
kappa

$$\renewcommand{\Bold}[1]{\mathbf{#1}}-\frac{m^{3} - m r^{2}}{m^{4} + 2 \, m^{2} r^{2} + r^{4}}$$

kappa.factor()

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

(acc/kappa).display()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( \frac{m^{2} + 2 \, m r + r^{2}}{2 \, {\left(m^{2} + r^{2}\right)}} \right) \frac{\partial}{\partial {\tilde{t}} } + \left( \frac{m^{2} - 2 \, m r + r^{2}}{2 \, {\left(m^{2} + r^{2}\right)}} \right) \frac{\partial}{\partial r } + \left( \frac{m}{m^{2} + r^{2}} \right) \frac{\partial}{\partial {\tilde{\varphi}} }$$

acc == kappa*el

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}$$

Plot of principal null geodesics

lamb = var('lamb', latex_name=r'\lambda')
def inPNG(v0, th0, tph0):
    return M.curve({KC: [lamb + v0, -lamb, th0, tph0]}, param=lamb)

def outPNG_I(u0, th0, tph0):
    return M.curve({KC: [u0 + r - 4*m^2/(r - m) + 4*m*ln(abs(r - m)/m), 
                         r, th0, tph0]}, 
                   param=(r, 1, +oo))

def outPNG_III(u0, th0, tph0):
    return M.curve({KC: [u0 + r - 4*m^2/(r - m) + 4*m*ln(abs(r - m)/m), 
                         r, th0, tph0]}, 
                   param=(r, -oo, 1))

u, v = var('u v')
tt_in(r,v) = v - r
tt_in

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( r, v \right) \ {\mapsto} \ -r + v$$

tt_out(r, u) = u + r - 4/(r - 1) + 4*ln(abs(r - 1))
tt_out

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left( r, u \right) \ {\mapsto} \ r + u - \frac{4}{r - 1} + 4 \, \log\left({\left| r - 1 \right|}\right)$$

rmin = -8; rmax = 8
graph = Graphics()
for u0 in range(-40, 40, 2):
    graph += plot(tt_out(r, u0), (r, rmin, rmax), color='green', 
                  axes_labels=[r"$r/m$", r"$\tilde{t}/m$"])
for v0 in range(-20, 20, 2):
    graph += plot(tt_in(r, v0), (r, rmin, rmax), color='green', linestyle='--')

graph += Hor
graph.set_aspect_ratio(1)
graph.show(ymin=-5, ymax=5, figsize=8)

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

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

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

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

el0 = el.at(p0)
print(el0)

Tangent vector el at Point p_0 on the 4-dimensional Lorentzian manifold M

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("exk_princ_null_geod.pdf", ymin=-5, ymax=5, figsize=8)
graph.show(ymin=-5, ymax=5, figsize=8)

umt(r) = u - tt_out(r, u)
umt(r)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}-r + \frac{4}{r - 1} - 4 \, \log\left({\left| r - 1 \right|}\right)$$

graph = plot(umt(r), (r, -8, 8), color='red', thickness=2, 
             axes_labels=[r"$r/m$", r"$(u - \tilde{t})/m$"], frame=True, 
             gridlines=True)
graph += line(((1, -10), (1, 10)), color='black', thickness=1.5, linestyle='--') 
show(graph, aspect_ratio=1, ymin=-8, ymax=8)

ttphi(r) = 2 / (r - 1)
ttphi(r)

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{2}{r - 1}$$

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(((1, -10), (1, 10)), color='black', thickness=1.5, linestyle='--')
show(graph, aspect_ratio=1, ymin=-8, ymax=8)

Carter-Penrose diagram

uc = (tt - r)/m + 4*m/(r - m) - 4*ln(abs((r - m)/m))
uc 

$$\renewcommand{\Bold}[1]{\mathbf{#1}}-\frac{4 \, m}{m - r} - \frac{r - {\tilde{t}}}{m} - 4 \, \log\left({\left| \frac{m - r}{m} \right|}\right)$$

vc = (tt + r)/m
vc

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\frac{r + {\tilde{t}}}{m}$$

Functions $U(\tilde{t}, r)$ and $V(\tilde{t}, r)$

U(tt, r) = atan(uc/2) + pi*unit_step(m - r)
V(tt, r) = atan(vc/2)

#U(tt, r) = tanh(uc/4) + 2*unit_step(m - r)
#V(tt, r) = tanh(vc/4)

#U(tt, r) = atan(asinh(uc))
#V(tt, r) = atan(asinh(vc))

#U(tt, r) = atan(uc^3)
#V(tt, r) = atan(vc^3)

#U(tt, r) = atan(real_nth_root(uc, 3))
#V(tt, r) = atan(real_nth_root(vc, 3))

plot(U(0, r).subs({m: 1}), (r, -15, 10), axes_labels=[r'$r/m$', r'$U, V$'],
     legend_label=r'$U\  (\tilde{t}=0)$') \
+ plot(U(m, r).subs({m: 1}), (r, -15, 10), linestyle='--', 
       legend_label=r'$U\  (\tilde{t}=m)$') \
+ plot(U(-m, r).subs({m: 1}), (r, -15, 10), linestyle=':',
       legend_label=r'$U\  (\tilde{t}=-m)$') \
+ plot(V(0, r).subs({m: 1}), (r, -15, 10), color='red',
       legend_label=r'$V\  (\tilde{t}=0)$') \
+ plot(V(m, r).subs({m: 1}), (r, -15, 10), color='red', linestyle='--',
       legend_label=r'$V\  (\tilde{t}=m)$') \
+ plot(V(-m, r).subs({m: 1}), (r, -15, 10), color='red', linestyle=':',
       legend_label=r'$V\  (\tilde{t}=-m)$')

Zoom around $r=m$:

plot(U(0, r).subs({m: 1}), (r, 0.8, 1.2), axes_labels=[r'$r/m$', r'$U$']) 

$U = V$ for $r= 2m$:

U(tt, 2*m).simplify_full()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\arctan\left(\frac{2 \, m + {\tilde{t}}}{2 \, m}\right)$$

bool(_ == V(tt, 2*m))

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\mathrm{True}$$

Functions $T(\tilde{t}, r)$ and $X(\tilde{t}, r)$

Tf(tt, r) = U(tt, r) + V(tt, r)
Xf(tt, r) = V(tt, r) - U(tt, r)

plot(Tf(0, r).subs({m: 1}), (r, -15, 10), axes_labels=[r'$r/m$', r'$T, X$'],
     legend_label=r'$T\  (\tilde{t}=0)$') \
+ plot(Tf(m, r).subs({m: 1}), (r, -15, 10), linestyle='--', 
       legend_label=r'$T\  (\tilde{t}=m)$') \
+ plot(Tf(-m, r).subs({m: 1}), (r, -15, 10), linestyle=':',
       legend_label=r'$T\  (\tilde{t}=-m)$') \
+ plot(Xf(0, r).subs({m: 1}), (r, -15, 10), color='red',
       legend_label=r'$X\  (\tilde{t}=0)$') \
+ plot(Xf(m, r).subs({m: 1}), (r, -15, 10), color='red', linestyle='--',
       legend_label=r'$X\  (\tilde{t}=m)$') \
+ plot(Xf(-m, r).subs({m: 1}), (r, -15, 10), color='red', linestyle=':',
       legend_label=r'$X\  (\tilde{t}=-m)$')

Zoom around $r=m$:

plot(Tf(0, r).subs({m: 1}), (r, 0.9, 1.1), axes_labels=[r'$r/m$', r'$T$'])

plot(diff(Tf(0, r), r).subs({m: 1}), (r, 0.9, 1.1), 
     axes_labels=[r'$r/m$', r'$\frac{\partial T}{\partial r}$']) 

Same value of $X$ for all $\tilde{t}$ for $r = r_1$:

s = -2*r + 4/(r -1) - 4*ln(1 -r)
s

$$\renewcommand{\Bold}[1]{\mathbf{#1}}-2 \, r + \frac{4}{r - 1} - 4 \, \log\left(-r + 1\right)$$

r1 = find_root(s, -4, -3)
r1

$$\renewcommand{\Bold}[1]{\mathbf{#1}}-3.4273344950702875$$

Compactified coordinates

CC.<T,X,th,tph> = M.chart(r"T:(-pi,2*pi) X:(-2*pi,pi) th:(0,pi):\theta tph:(0,2*pi):periodic:\tilde{\varphi}")
CC

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left(M,(T, X, {\theta}, {\tilde{\varphi}})\right)$$

CC.coord_range()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}T :\ \left( -\pi , 2 \, \pi \right) ;\quad X :\ \left( -2 \, \pi , \pi \right) ;\quad {\theta} :\ \left( 0 , \pi \right) ;\quad {\tilde{\varphi}} :\ \left[ 0 , 2 \, \pi \right] \mbox{(periodic)}$$

KC_to_CC = KC.transition_map(CC, [Tf(tt, r), Xf(tt, r), th, tph])
KC_to_CC.display()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left\{\begin{array}{lcl} T & = & \pi \mathrm{u}\left(m - r\right) + \arctan\left(-\frac{2 \, m}{m - r} - \frac{r - {\tilde{t}}}{2 \, m} - 2 \, \log\left({\left| \frac{m - r}{m} \right|}\right)\right) + \arctan\left(\frac{r + {\tilde{t}}}{2 \, m}\right) \\ X & = & -\pi \mathrm{u}\left(m - r\right) - \arctan\left(-\frac{2 \, m}{m - r} - \frac{r - {\tilde{t}}}{2 \, m} - 2 \, \log\left({\left| \frac{m - r}{m} \right|}\right)\right) + \arctan\left(\frac{r + {\tilde{t}}}{2 \, m}\right) \\ {\theta} & = & {\theta} \\ {\tilde{\varphi}} & = & {\tilde{\varphi}} \end{array}\right.$$

graph0 = polygon([(0, pi), (-pi, 2*pi), (-2*pi, pi), (-pi, 0)], 
                color='cornsilk', edgecolor='black') \
         + polygon([(pi, 0), (0, pi), (-pi, 0), (0, -pi)], 
                   color='white', edgecolor='black')
graph0

Hor = line([(-pi,0), (0, pi)], color='black', thickness=3) \
      + text(r'$\mathscr{H}$', (-0.6, 3), color='black', fontsize=20)

def plot_const_r(r0, color='red', linestyle=':', thickness=1, plot_points=400):
    return KC.plot(CC, ambient_coords=(X,T), fixed_coords={th: pi/3, tph: 0, r: r0},
                   ranges={tt: (-100, 100)}, color=color, style=linestyle,
                   thickness=thickness, plot_points=plot_points, parameters={m: 1})

def plot_const_tt(tt0, color='darkgrey', linestyle='-', thickness=1, plot_points=100):
    resu = KC.plot(CC, ambient_coords=(X,T), fixed_coords={th: pi/3, tph: 0, tt: tt0},
                   ranges={r: (-100, -10)}, color=color, style=linestyle,
                   thickness=thickness, plot_points=plot_points, parameters={m: 1}) \
           + KC.plot(CC, ambient_coords=(X,T), fixed_coords={th: pi/3, tph: 0, tt: tt0},
                     ranges={r: (-10, 10)}, color=color, style=linestyle,
                     thickness=thickness, plot_points=plot_points, parameters={m: 1}) \
           + KC.plot(CC, ambient_coords=(X,T), fixed_coords={th: pi/3, tph: 0, tt: tt0},
                     ranges={r: (10, 100)}, color=color, style=linestyle,
                     thickness=thickness, plot_points=plot_points, parameters={m: 1})
    return resu

Plot of $r=\mathrm{const}$ curves:

graph_r = Graphics()

rmin, rmax = -10, 0
dr = 2
nr = (rmax - rmin)/dr
for i in range(int(nr)):
    ri = rmin + dr*i
    graph_r += plot_const_r(ri) 

rmin, rmax = 0, 0.8
dr = 0.2
nr = (rmax - rmin)/dr + 1
for i in range(int(nr)):
    ri = rmin + dr*i
    graph_r += plot_const_r(ri) 

rmin, rmax = 1.2, 3
dr = 0.2
nr = (rmax - rmin)/dr + 1
for i in range(int(nr)):
    ri = rmin + dr*i
    graph_r += plot_const_r(ri) 

rmin, rmax = 3, 13
dr = 2
nr = (rmax - rmin)/dr + 1
for i in range(int(nr)):
    ri = rmin + dr*i
    graph_r += plot_const_r(ri) 
graph_r += plot_const_r(0, color='maroon', linestyle='--', thickness=2)

graph1 = graph0 + plot_const_tt(0, color='darkgrey', thickness=2)

tmin, tmax = -20, 20
dt = 2
nt = (tmax - tmin)/dt
for i in range(nt):
    ti = tmin + dt*i
    graph1 += plot_const_tt(ti) 

graph1 += graph_r + Hor
show(graph1, figsize=10)

Adding some principal null geodesics to the plot:

graph_PNG = Graphics()

for L in [inPNG(0, pi/3, 0), inPNG(-4, pi/3, 0), inPNG(4, pi/3, 0)]:
    L.expr(chart2=CC)
    graph_PNG += L.plot(CC, ambient_coords=(X, T), color='green', style='--', 
                        max_range=100, plot_points=5, parameters={m: 1})
    
for L in [outPNG_I(0, pi/3, 0), outPNG_I(-4, pi/3, 0), outPNG_I(4, pi/3, 0)]:
    L.expr(chart2=CC)
    graph_PNG += L.plot(CC, ambient_coords=(X, T), color='green', 
                        prange=(1.001, 100), plot_points=5, parameters={m: 1})

for L in [outPNG_III(0, pi/3, 0), outPNG_III(-4, pi/3, 0), outPNG_III(4, pi/3, 0)]:
    L.expr(chart2=CC)
    graph_PNG += L.plot(CC, ambient_coords=(X, T), color='green', 
                        prange=(-100, 0.999), plot_points=5, parameters={m: 1})

graph1 += graph_PNG
show(graph1, figsize=10, axes=False, frame=True)

We save the plot in png format, in order to add easily labels on it with Inkscape:

graph1.save('exk_CPdiag_Kerr-raw.svg', figsize=10, axes=False, frame=True)

Carter-Penrose diagram with Boyer-Lindquist time slices

BL_to_CC = KC_to_CC * KC_to_BL.inverse()
BL_to_CC.display()

$$\renewcommand{\Bold}[1]{\mathbf{#1}}\left\{\begin{array}{lcl} T & = & \pi \mathrm{u}\left(m - r\right) + \arctan\left(\frac{2 \, m^{2} \log\left(m\right) - 2 \, m^{2} - {\left(2 \, m \log\left(m\right) + m\right)} r + r^{2} + {\left(m - r\right)} t - 2 \, {\left(m^{2} - m r\right)} \log\left({\left| -m + r \right|}\right)}{2 \, {\left(m^{2} - m r\right)}}\right) + \arctan\left(-\frac{2 \, m^{2} \log\left(m\right) - 2 \, m^{2} - {\left(2 \, m \log\left(m\right) + m\right)} r + r^{2} - {\left(m - r\right)} t - 2 \, {\left(m^{2} - m r\right)} \log\left({\left| -m + r \right|}\right)}{2 \, {\left(m^{2} - m r\right)}}\right) \\ X & = & -\pi \mathrm{u}\left(m - r\right) - \arctan\left(\frac{2 \, m^{2} \log\left(m\right) - 2 \, m^{2} - {\left(2 \, m \log\left(m\right) + m\right)} r + r^{2} + {\left(m - r\right)} t - 2 \, {\left(m^{2} - m r\right)} \log\left({\left| -m + r \right|}\right)}{2 \, {\left(m^{2} - m r\right)}}\right) + \arctan\left(-\frac{2 \, m^{2} \log\left(m\right) - 2 \, m^{2} - {\left(2 \, m \log\left(m\right) + m\right)} r + r^{2} - {\left(m - r\right)} t - 2 \, {\left(m^{2} - m r\right)} \log\left({\left| -m + r \right|}\right)}{2 \, {\left(m^{2} - m r\right)}}\right) \\ {\theta} & = & {\theta} \\ {\tilde{\varphi}} & = & \frac{m {\varphi} - {\varphi} r + m}{m - r} \end{array}\right.$$

def plot_const_t(t0, color='blue', linestyle='-', thickness=1, plot_points=50):
    resu = BL.plot(CC, ambient_coords=(X,T), fixed_coords={th: pi/3, ph: 0, t: t0},
                   ranges={r: (-100, -10)}, color=color, style=linestyle,
                   thickness=thickness, plot_points=plot_points, parameters={m: 1}) \
           + BL.plot(CC, ambient_coords=(X,T), fixed_coords={th: pi/3, ph: 0, t: t0},
                     ranges={r: (-10, 0)}, color=color, style=linestyle,
                     thickness=thickness, plot_points=plot_points, parameters={m: 1}) \
           + BL.plot(CC, ambient_coords=(X,T), fixed_coords={th: pi/3, ph: 0, t: t0},
                     ranges={r: (0, 0.999)}, color=color, style=linestyle,
                     thickness=thickness, plot_points=plot_points, parameters={m: 1}) \
           + BL.plot(CC, ambient_coords=(X,T), fixed_coords={th: pi/3, ph: 0, t: t0},
                     ranges={r: (0.999, 3)}, color=color, style=linestyle,
                     thickness=thickness, plot_points=plot_points, parameters={m: 1}) \
           + BL.plot(CC, ambient_coords=(X,T), fixed_coords={th: pi/3, ph: 0, t: t0},
                     ranges={r: (3, 100)}, color=color, style=linestyle,
                     thickness=thickness, plot_points=plot_points, parameters={m: 1})
    return resu

tmin, tmax = -20, 20
dt = 2
nt = (tmax - tmin)/dt
graph2 = graph0
for i in range(nt):
    ti = tmin + dt*i
    graph2 += plot_const_t(ti) 

graph2 += graph_r + Hor

graph2.save('exk_CPdiag_BL-raw.svg', figsize=10, axes=False, frame=True)
show(graph2, figsize=10, axes=False, frame=True)