SageMath notebooks associated to the Black Hole Lectures (https://luth.obspm.fr/~luthier/gourgoulhon/bh16)
Maximal extension of the 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()
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)
We then introduce (3+1 version of) the Kerr coordinates (t~,r,θ,φ~) 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 (−∞,+∞)) 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
KC.coord_range()
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()
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()
Regions MI and MIII
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()
Boyer-Lindquist coordinates on MI
Let us introduce on the chart of Boyer-Lindquist coordinates (t,r,θ,φ) on MI:
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
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)
Ingoing principal null geodesics
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 ∇kk=0:
nabla = g.connection()
nabla(k).contract(k).display()
Expression of k with respect to the Boyer-Lindquist frame:
k.display(BL)
Outgoing principal null geodesics
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 ℓ is a null vector:
g(el, el).expr()
Expression of ℓ with respect to the Boyer-Lindquist frame:
el.display(BL)
Computation of ∇ℓℓ:
acc = nabla(el).contract(el) acc.display()
We check that ∇ℓℓ=κℓ:
kappa = acc[0] / el[0] kappa
kappa.factor()
acc == kappa*el
Outgoing Kerr coordinates on MI
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 ℓ′=2(r−m)2r2+m2ℓ:
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()
Non-affinity coefficient of k′
acc_ok = nabla(ok).contract(ok) acc_ok.display()
kappa_ok = acc_ok[0] / ok[0] kappa_ok.factor()
We check that ∇k′k′=κk′k′:
acc_ok == kappa_ok * ok
Compactified coordinates on M
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()
Spacetime (M′,g)
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()
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(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
Plots of hypersurfaces of constant r
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)
Plots of hypersurfaces of constant t~
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)
Plots of hypersurfaces of constant t~~
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)
Maximal extension
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)
Near-horizon region
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)