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)
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))
KC.coord_range()
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()
A matrix view of the components with respect to the manifold's default vector frame:
g[:]
The list of the non-vanishing components:
g.display_comp()
Let us check that we are dealing with a solution of the vacuum Einstein equation:
g.ricci().display()
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))
BL.coord_range()
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()
KC_to_BL.inverse().display()
g.display(BL)
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)
k = M.vector_field(1, -1, 0, 0, name='k')
k.display()
Let us check that $k$ is a null vector:
g(k, k).expr()
Check of $\nabla_k k = 0$:
nabla = g.connection()
nabla(k).contract(k).display()
Expression of $k$ with respect to the Boyer-Lindquist frame:
k.display(BL)
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()
Let us check that $\ell$ is a null vector:
g(el, el).expr()
Expression of $\ell$ with respect to the Boyer-Lindquist frame:
el.display(BL)
Computation of $\nabla_\ell \ell$:
acc = nabla(el).contract(el)
acc.display()
We check that $\nabla_\ell \ell = \kappa \ell$:
kappa = acc[0] / el[0]
kappa
kappa.factor()
(acc/kappa).display()
acc == kappa*el
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
tt_out(r, u) = u + r - 4/(r - 1) + 4*ln(abs(r - 1))
tt_out
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
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)
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)
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)
uc = (tt - r)/m + 4*m/(r - m) - 4*ln(abs((r - m)/m))
uc
vc = (tt + r)/m
vc
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()
bool(_ == V(tt, 2*m))
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
r1 = find_root(s, -4, -3)
r1
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
CC.coord_range()
KC_to_CC = KC.transition_map(CC, [Tf(tt, r), Xf(tt, r), th, tph])
KC_to_CC.display()
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)
BL_to_CC = KC_to_CC * KC_to_BL.inverse()
BL_to_CC.display()
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)