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.5.beta2, Release Date: 2021-09-26'
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()
M_I = M.open_subset('M_I', latex_name=r'M_{\rm I}', coord_def={KC: r>m})
KC.restrict(M_I).coord_range()
M_III = M.open_subset('M_III', latex_name=r'M_{\rm III}', coord_def={KC: r<m})
KC.restrict(M_III).coord_range()
Let us introduce on the chart of Boyer-Lindquist coordinates $(t,r,\theta,\varphi)$ on $M_{\rm I}$:
BL.<t,r,th,ph> = M_I.chart(r"t r:(m,+oo) th:(0,pi):\theta ph:(0,2*pi):periodic:\varphi")
print(BL); BL
Chart (M_I, (t, r, th, ph))
BL.coord_range()
KC_to_BL = KC.restrict(M_I).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)
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 that $k$ is a geodesic vector field, i.e. obeys $\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((r + m)^2/(2*(r^2 + m^2)),
(r - m)^2/(2*(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*el
OKC.<to,r,th,oph> = M_I.chart(r"to:\tilde{\tilde{t}} r:(m,+oo) th:(0,pi):\theta oph:(0,2*pi):periodic:\tilde{\tilde{\varphi}}")
OKC.coord_range()
BL_to_OKC = BL.transition_map(OKC, [t + 2*m^2/(r-m) - 2*m*ln(abs(r-m)/m),
r, th, ph + m/(r-m)])
BL_to_OKC.display()
BL_to_OKC.inverse().display()
KC_to_OKC = BL_to_OKC * KC_to_BL.restrict(M_I)
KC_to_OKC.display()
KC_to_OKC.inverse().display()
M_I.set_default_chart(OKC)
M_I.set_default_frame(OKC.frame())
gI = g.restrict(M_I)
gI.display()
gI[1,3]
gI[1,3] == m*(1 + 2*m*r/rho2)*sin(th)^2
gI[3,3]
g[3,3] == (r^2 + m^2 + 2*m^3*r*sin(th)^2/rho2)*sin(th)^2
ol = M_I.vector_field({OKC.frame(): (1, 1, 0, 0)}, name='ol',
latex_name=r"\ell'")
ol.display()
ol.display(KC.restrict(M_I).frame())
ol.display(BL.frame())
g(ol, ol).expr()
nabla.coef(OKC.frame())
nabla(ol).contract(ol).display()
elI = el.restrict(M_I)
elI.display()
Check of the relation $\ell' = 2 \frac{r^2 + m^2}{(r - m)^2} \, \ell$:
ol == 2*(r^2 + m^2)/(r - m)^2 * elI
kI = k.restrict(M_I)
kI.display()
ok = (r - m)^2/(2*(r^2 + m^2)) * kI
ok.set_name('ok', latex_name=r"k'")
ok.display()
g(k, el).expr()
g(ok, ol).expr()
g(k, ol).expr().factor()
g(ok, el).expr().factor()
acc_ok = nabla(ok).contract(ok)
acc_ok.display()
kappa_ok = acc_ok[0] / ok[0]
kappa_ok.factor()
We check that $\nabla_{k'} k' = \kappa_{k'} k'$:
acc_ok == kappa_ok * ok
CC.<T,X,th,tph> = M.chart(r"T X th:(0,pi):\theta tph:(0,2*pi):periodic:\tilde{\varphi}")
CC
uc = (tt - r)/m + 4*m/(r - m) - 4*ln(abs((r - m)/m))
vc = (tt + r)/m
KC_to_CC = KC.transition_map(CC, [atan(uc/2) + atan(vc/2) + pi*unit_step(m - r),
atan(vc/2) - atan(uc/2) - pi*unit_step(m - r),
th,
tph])
KC_to_CC.display()
OKC_to_CC = KC_to_CC.restrict(M_I) * KC_to_OKC.inverse()
OKC_to_CC.display()
forget(r>m)
Mp = Manifold(4, "M'", structure='Lorentzian')
OKCp.<to,r,th,oph> = Mp.chart(r"to:\tilde{\tilde{t}} r th:(0,pi):\theta oph:(0,2*pi):periodic:\tilde{\tilde{\varphi}}")
OKCp
OKCp.coord_range()
CCp.<T,X,th,oph> = Mp.chart(r"T X th:(0,pi):\theta oph:(0,2*pi):periodic:\tilde{\tilde{\varphi}}")
CCp
uc = (to - r)/m
vc = (to + r)/m - 4*m/(r - m) + 4*ln(abs((r - m)/m))
OKC_to_CCp = OKCp.transition_map(CCp, [atan(uc/2) + atan(vc/2) - pi*unit_step(m - r),
atan(vc/2) - atan(uc/2) - pi*unit_step(m - r),
th,
oph])
OKC_to_CCp.display()
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(u0, th0, oph0):
return Mp.curve({OKCp: [u0 + r, r, th0, oph0]}, param=r)
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))
def inPNG_IIIp(v0, th0, oph0):
return Mp.curve({OKCp: [v0 + lamb - 4*m^2/(lamb + m) - 4*m*ln(abs(lamb + m)/m),
-lamb, th0, oph0]},
param=(lamb, -1, +oo))
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') \
+ polygon([(0, -pi), (-pi, 0), (-2*pi, -pi), (-pi, -2*pi)],
color='cornsilk', edgecolor='black')
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=4, parameters={m: 1})
for L in [outPNG(0, pi/3, 0), outPNG(-4, pi/3, 0), outPNG(4, pi/3, 0)]:
L.expr(chart2=CCp)
graph_PNG += L.plot(CCp, ambient_coords=(X, T), color='green',
max_range=100, plot_points=4, 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=4, parameters={m: 1})
for L in [inPNG_IIIp(0, pi/3, 0), inPNG_IIIp(-4, pi/3, 0), inPNG_IIIp(4, pi/3, 0)]:
L.expr(chart2=CCp)
graph_PNG += L.plot(CCp, ambient_coords=(X, T), color='green', style='--',
prange=(-0.999, 100), plot_points=4, parameters={m: 1})
graph = graph0 + graph_PNG
graph
def plot_const_r(r0, color='red', linestyle=':', thickness=1, plot_points=300):
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
def plot_const_r_p(r0, color='red', linestyle=':', thickness=1, plot_points=300):
return OKCp.plot(CCp, ambient_coords=(X,T), fixed_coords={th: pi/3, oph: 0, r: r0},
ranges={to: (-100, 100)}, color=color, style=linestyle,
thickness=thickness, plot_points=plot_points, parameters={m: 1})
def plot_const_to(to0, color='violet', linestyle='-', thickness=1, plot_points=100):
resu = OKCp.plot(CCp, ambient_coords=(X,T), fixed_coords={th: pi/3, oph: 0, to: to0},
ranges={r: (-100, -10)}, color=color, style=linestyle,
thickness=thickness, plot_points=plot_points, parameters={m: 1}) \
+ OKCp.plot(CCp, ambient_coords=(X,T), fixed_coords={th: pi/3, oph: 0, to: to0},
ranges={r: (-10, 10)}, color=color, style=linestyle,
thickness=thickness, plot_points=plot_points, parameters={m: 1}) \
+ OKCp.plot(CCp, ambient_coords=(X,T), fixed_coords={th: pi/3, oph: 0, to: to0},
ranges={r: (10, 100)}, color=color, style=linestyle,
thickness=thickness, plot_points=plot_points, parameters={m: 1})
return resu
for r0 in [-10, -8, -6, -4, -2, 0.5, 1.5, 2, 2.5, 3, 5, 7, 9]:
graph += plot_const_r(r0)
graph += plot_const_r(0, color='maroon', linestyle='--', thickness=2)
for r0 in [-10, -8, -6, -4, -2, 0.5]:
graph += plot_const_r_p(r0)
graph += plot_const_r_p(0, color='maroon', linestyle='--', thickness=2)
show(graph, figsize=10, axes=False, frame=True)
tmin, tmax, dt = -10, 10, 2
for i in range(int((tmax - tmin)/dt) + 1):
ti = tmin + dt*i
graph += plot_const_tt(ti)
graph += plot_const_tt(0, thickness=2)
show(graph, figsize=10, axes=False, frame=True)
tmin, tmax, dt = -10, 10, 2
for i in range(int((tmax - tmin)/dt) + 1):
ti = tmin + dt*i
graph += plot_const_to(ti)
graph += plot_const_to(0, thickness=2)
graph += line([(-pi,0), (0, pi)], color='black', thickness=3) \
+ line([(-pi,0), (0, -pi)], color='black', thickness=3)
show(graph, figsize=10, axes=False, frame=True)
graph.save('exk_CPdiag_M0-raw.svg', figsize=10, axes=False, frame=True)
Tc(tt, r) = KC_to_CC(tt, r, th, tph)[0].subs(m=1).simplify_full()
Xc(tt, r) = KC_to_CC(tt, r, th, tph)[1].subs(m=1).simplify_full()
show(Tc)
show(Xc)
TBL(t, r) = Tc(t - 2/(r-1) + 2*ln(abs(r-1)), r).simplify_full()
XBL(t, r) = Xc(t - 2/(r-1) + 2*ln(abs(r-1)), r).simplify_full()
show(TBL)
show(XBL)
def plot_I(n, t_values=None, r_min=1.0001, r_max=100,
color_t='dimgray', linestyle_t='-',
r_values=None, t_min=-100, t_max=100,
color_r='red', linestyle_r=':',
plot_null_geod=True, hor_thickness=3):
n2 = 2*n
res = polygon([(pi, n2*pi), (0, (n2 + 1)*pi), (-pi, n2*pi), (0, (n2 - 1)*pi)],
color='white', edgecolor='black')
if r_values is not None:
for r0 in r_values:
res += parametric_plot((Xc(tt, r0), Tc(tt, r0) + n2*pi), (tt, t_min, t_max),
color=color_r, linestyle=linestyle_r)
if t_values is not None:
for t0 in t_values:
res += parametric_plot((XBL(t0, r), TBL(t0, r) + n2*pi), (r, r_min, r_max),
color=color_t, linestyle=linestyle_t)
if plot_null_geod:
res += line([(pi/2, -pi/2 + n2*pi), (-pi/2, pi/2 + n2*pi)], color='green',
linestyle='--')
res += line([(-pi/2, -pi/2 + n2*pi), (pi/2, pi/2 + n2*pi)], color='green')
res += line([(-pi, n2*pi), (0, (n2 + 1)*pi)], color='black', thickness=hor_thickness)
res += line([(-pi, n2*pi), (0, (n2 - 1)*pi)], color='black', thickness=hor_thickness)
return res
def plot_III(n, t_values=None, r_min=-100, r_max=0.9999,
color_t='dimgray', linestyle_t='-',
r_values=None, t_min=-100, t_max=100,
color_r='red', linestyle_r=':',
plot_null_geod=True):
n2 = 2*n
res = polygon([(0, (n2 + 1)*pi), (-pi, (n2 + 2)*pi), (-2*pi, (n2 + 1)*pi), (-pi, n2*pi)],
color='cornsilk', edgecolor='black')
res += parametric_plot((Xc(tt, 0), Tc(tt, 0) + n2*pi), (tt, t_min, t_max), color='maroon',
linestyle='--', thickness=2)
if r_values is not None:
for r0 in r_values:
res += parametric_plot((Xc(tt, r0), Tc(tt, r0) + n2*pi), (tt, t_min, t_max),
color=color_r, linestyle=linestyle_r)
if t_values is not None:
for t0 in t_values:
res += parametric_plot((XBL(t0, r), TBL(t0, r) + n2*pi), (r, r_min, r_max),
color=color_t, linestyle=linestyle_t)
if plot_null_geod:
res += line([(-pi/2, pi/2 + n2*pi), (-3*pi/2, 3*pi/2 + n2*pi)], color='green',
linestyle='--')
res += line([(-3*pi/2, pi/2 + n2*pi), (-pi/2, 3*pi/2 + n2*pi)], color='green')
return res
r_val_I = [1.2, 1.4, 1.6, 1.8, 2, 4, 6, 8, 10]
r_val_III = [-10, -8, -6, -4, -2, 0.2, 0.4, 0.6, 0.8]
show(graphics_array([plot_I(0, r_values=r_val_I),
plot_III(0, r_values=r_val_III)]), axes=True)
graph = Graphics()
for n in [-2..2]:
graph += plot_I(n, r_values=r_val_I) + plot_III(n, r_values=r_val_III)
show(graph, figsize=10, axes=False, ymin=-3*pi, ymax=3*pi)
graph.save('exk_CPdiag_maximal-raw.svg', figsize=10, axes=False,
ymin=-3*pi, ymax=3*pi)
r_val_I = [1.1, 1.2, 1.3, 1.4, 1.5]
r_val_III = [0.9, 0.8, 0.7, 0.6, 0.5]
t_val = [-8, -6, -4, -2, 0, 2, 4, 6, 8]
graph = Graphics()
for n in [-1..1]:
graph += plot_I(n, t_values=t_val, r_max=1.5,
r_values=r_val_I, linestyle_r='-',
plot_null_geod=False, hor_thickness=4)
graph += plot_III(n, t_values=t_val, r_min=0.5,
r_values=r_val_III, linestyle_r='-',
plot_null_geod=False)
show(graph, figsize=8, axes=False, frame=True, ymin=-pi, ymax=2*pi)
graph += text(r"$r=\!+\infty$", (pi/2+0.2, pi/2+0.2), rotation=-45, color='black',
fontsize=16) \
+ text(r"$r=\!+\infty$", (pi/2+0.2, -pi/2-0.2), rotation=45, color='black',
fontsize=16) \
+ text(r"$r=\!+\infty$", (pi/2+0.2, 3*pi/2-0.2), rotation=45, color='black',
fontsize=16) \
+ text(r"$r=\!-\infty$", (-3*pi/2-0.2, pi/2-0.2), rotation=-45, color='black',
fontsize=16) \
+ text(r"$r=\!-\infty$", (-3*pi/2-0.2, 3*pi/2+0.2), rotation=45, color='black',
fontsize=16) \
+ text(r"$r=\!-\infty$", (-3*pi/2-0.2, -pi/2+0.2), rotation=45, color='black',
fontsize=16) \
+ text(r"$r=m$", (-pi/2-0.2, pi/2+0.2), rotation=45, color='black',
fontsize=16) \
+ text(r"$r=m$", (-pi/2+0.15, 3*pi/2+0.15), rotation=-45, color='black',
fontsize=16) \
+ text(r"$r=m$", (-pi/2+0.15, -pi/2+0.15), rotation=-45, color='black',
fontsize=16) \
+ text(r"$r=0$", (-2.4, 2), rotation=55, color='maroon',
fontsize=16) \
+ text(r"$\mathscr{H}$", (-0.25, 2.4), color='black', fontsize=20) \
+ text(r"$\mathscr{H}'$", (-0.25, -2.2), color='black', fontsize=20) \
+ text(r"$i^{\rm int}$", (-3.6, -0.1), color='orangered', fontsize=16) \
+ text(r"$\mathscr{M}_{\rm I}$", (0, 0), color='black', fontsize=24) \
+ text(r"$\mathscr{M}_{\rm III}$", (-pi, pi), color='black', fontsize=24)
graph.save("exk_NH_region.pdf", figsize=9, axes=False, ymin=-pi, ymax=2*pi)
show(graph, figsize=9, axes=False, ymin=-pi, ymax=2*pi)