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 10.0.beta8, Release Date: 2023-04-06'
%display latex
Due to spherical symmetry, we a considering a 2-dimensional cut at $(\theta,\varphi) = \mathrm{const}$:
M = Manifold(2, 'M', structure='Lorentzian')
EF.<ti, r> = M.chart(r"ti:\tilde{t} r:[0,+oo)")
EF
EF.coord_range()
I = M.open_subset('I')
We choose $\chi_{\rm s} = \pi/4$:
S.<eta,chi> = I.chart(r"eta:[0,pi]:\eta chi:[0,pi/4]:\chi")
S
S.coord_range()
S.coord_bounds()
chis = S.coord_bounds()[1][1][0]
chis
The initial areal radius of the star in units of $m$:
r0 = 2/sin(chis)^2
r0
$r_{\rm s}(\eta)$ in units of $r_0$:
rsurf(eta) = (1 + cos(eta))/2
rsurf
Matter proper time $\tau$ in units of $\sqrt{r_0^3/(8 m_0)}$:µ
tau(eta) = (eta + sin(eta))/2
tau
Proper energy density $\rho$ in units of $6 m / (\pi r_0^3)$:
rho(eta) = 1 / (1 + cos(eta))^3
rho
graph = plot(tau, (0, pi), color='purple', thickness=2, linestyle='--',
legend_label=r'$\tau \; \sqrt{2 m / r_0^3}$',
axes_labels=[r'$\eta$', r'$r_{\rm s}$, $\tau$, $\rho$'],
ticks=[pi/4, None], tick_formatter=[pi, None], fontsize=12,
axes_labels_size=1.4,
frame=True, axes=False, gridlines=True) \
+ plot(rsurf, (0, pi), color='red', thickness=2,
legend_label=r'$r_{\rm s}(\tau)/r_0 = a(\tau)\, \sin\chi_{\rm s} / r_0$') \
+ plot(rho, (0, 2), color='blue', thickness=2,
legend_label=r'$\rho(\tau) \, \pi r_0^3/(6 m)$')
graph.set_legend_options(handlelength=2, labelspacing=0.5,
font_size=11)
graph.set_aspect_ratio(0.75)
graph.show(ymax=3.2)
graph.save('lem_OS_rs_rho_eta.pdf', ymax=3.2)
graph = parametric_plot((tau(eta), eta), (eta, 0, pi), color='orange',
linestyle='--', thickness=2,
legend_label=r'$\eta$',
axes_labels=[r'$\tau \; \sqrt{2m/r_0^3}$',
r'$\eta$, $r_{\rm s}$, $\rho$'],
frame=True, axes=False, gridlines=True) \
+ parametric_plot((tau(eta), rsurf(eta)), (eta, 0, pi), color='red',
thickness=2,
legend_label=r'$r_{\rm s}(\tau)/r_0 = a(\tau)\, \sin\chi_{\rm s} / r_0$') \
+ parametric_plot((tau(eta), rho(eta)), (eta, 0, 2), color='blue',
thickness=2,
legend_label=r'$\rho(\tau) \, \pi r_0^3/(6 m)$')
graph.set_legend_options(handlelength=2, labelspacing=0.5,
font_size=11)
graph.set_aspect_ratio(0.4)
graph.show(ymax=3.2)
graph.save('lem_OS_rs_rho_tau.pdf', ymax=3.2)
$\eta$ as a function of $\tau/m$:
def etaf(tau):
def ff(et):
return float((et + sin(et))/sin(chis)^3 - tau)
return find_root(ff, float(0), float(pi))
Value of the matter proper time $\tau$ at the end of the collapse, in units of $m$:
tau_end = pi/sin(chis)^3
tau_end, n(tau_end)
Check:
etaf(tau_end) == n(pi)
plot(etaf, 0, tau_end)
graph = plot(3/(4*pi)*sin(x)^6/(1 - cos(3*x))^3, (x, 0, pi/3),
color='blue', thickness=2,
axes_labels=[r'$\chi_{\rm s}$', r'$\rho_{\rm hb} m^2$'],
ticks=[pi/8, None], tick_formatter=[pi, None],
frame=True, axes=False, gridlines=True)
graph.show(ticks=pi/9)
graph.save("lem_OS_rho_hb.pdf", ticks=pi/9)
Constant $\chi$ lines
graph = S.plot(S, ambient_coords=(chi, eta), style={eta: '-', chi: '--'},
number_values={eta: 2, chi: 9})
Constant $\tau$ lines
tau_sel = [numerical_approx(i*tau_end/8) for i in range(9)]
tau_sel
eta_sel = [etaf(tau0) for tau0 in tau_sel]
eta_sel
iso_taus = [I.curve({S: [eta0, chi]}, (chi, 0, chis)) for eta0 in eta_sel]
for iso_tau in iso_taus:
graph += iso_tau.plot(chart=S, ambient_coords=(chi, eta), color='red',
style='--')
The stellar surface:
surf = I.curve({S: [eta, chis]}, (eta, 0, pi))
graph += surf.plot(S, ambient_coords=(chi, eta), thickness=3)
The star at $\eta=0$:
init= I.curve({S: [0, chi]}, (chi, 0, chis))
graph += init.plot(S, ambient_coords=(chi, eta),
color='magenta', thickness=3)
graph.show(figsize=8)
def geod_int_out(etas):
r"""
Internal outgoing radial null geodesic ending at
(eta, chi) = (etas, chis)
"""
chi_min = 0 if etas > chis else chis - etas
return I.curve({S: [chi - chis + etas, chi]}, (chi, chi_min, chis))
def geod_int_in(etas):
r"""
Internal ingoing radial null geodesic starting at
(eta, chi) = (etas, chis)
"""
chi_min = 0 if etas < pi - chis else etas - pi + chis
return I.curve({S: [chis - chi + etas, chi]}, (chi, chi_min, chis))
hor_int = geod_int_out(pi - 2*chis)
graph += hor_int.plot(chart=S, ambient_coords=(chi, eta),
color='black', thickness=3)
matter_sing = I.curve({S: [pi, chi]}, (chi, 0, chis))
graph += line([(0 + i*chis/8, pi) for i in range(9)],
thickness=3, color='darkorange',
marker='D', markersize=10)
graph.show(frame=True, axes_pad=0, figsize=8)
EFI = EF.restrict(I)
I.atlas()
S_to_EF = S.transition_map(EFI,
[2*(sqrt(r0/2 - 1)*(eta + r0/4*(eta + sin(eta)))
+ 2*ln(cos(eta/2) + sin(eta/2)/sqrt(r0/2 - 1))),
r0*sin(chi)/sin(chis)/2*(1 + cos(eta))])
S_to_EF.display()
S_to_EF(eta, chi)
rr(eta, chi) = S_to_EF(eta, chi)[1]
rr
Adding isocontours $r = \mathrm{const}$:
graph += contour_plot(rr(eta, chi), (chi, 0, chis), (eta, 0, 3.13),
contours=[0.5, 1., 1.5, 2., 2.5, 3., 3.5],
fill=False, labels=True, label_colors='black',
label_inline=True, label_inline_spacing=10)
graph.show(frame=True, axes_pad=0, ticks=[pi/8,pi/8], tick_formatter=(pi,pi),
fontsize=16, axes_labels_size=1.2, figsize=10)
ITH = I.curve({S: [pi - 2*chi, chi]}, (chi, 0, chis))
graph += ITH.plot(chart=S, ambient_coords=(chi, eta), color='blue',
thickness=3, style=':', label_axes=False)
graph_trap = copy(graph)
A selection of null geodesics:
geod_int_out_sel = [geod_int_out(etas) for etas in eta_sel]
geod_int_in_sel = [geod_int_in(etas) for etas in eta_sel]
for geod in geod_int_out_sel:
graph += geod.plot(chart=S, ambient_coords=(chi, eta),
color='green', thickness=2)
for geod in geod_int_in_sel[:-1]:
graph += geod.plot(chart=S, ambient_coords=(chi, eta),
color='green', thickness=1, style='--')
graph.show(frame=True, axes=False, axes_pad=0, ticks=[pi/8,pi/4],
tick_formatter=(pi,pi), fontsize=16, axes_labels_size=1.2,
xmax=pi/4+0.015, ymin=-0.015, figsize=10)
graph.save('lem_OS_diag_int.pdf', frame=True, axes=False, axes_pad=0,
ticks=[pi/8,pi/4], tick_formatter=(pi,pi), fontsize=16,
axes_labels_size=1.2, xmax=pi/4+0.015, ymin=-0.015, figsize=10)
for i in range(9):
etas = pi/8*i
graph_trap += geod_int_out(etas).plot(chart=S,
ambient_coords=(chi, eta),
color='green', thickness=2)
graph_trap += geod_int_in(etas).plot(chart=S,
ambient_coords=(chi, eta),
color='green', thickness=1,
style='--')
graph_trap.show(frame=True, axes_pad=0, ticks=[pi/8,pi/4], tick_formatter=(pi,pi),
fontsize=16, axes_labels_size=1.2, figsize=10)
Adding $r$ isocontours near the trapped surface at $(\eta,\chi) = \left(\frac{11\pi}{6}, \frac{3\pi}{16} \right)$:
tsp = I((11*pi/16, 3*pi/16), chart=S)
S(tsp)
EF(tsp)
rts = n(EF(tsp)[1])
rts
graph_trap += contour_plot(rr(eta, chi), (chi, 0, chis), (eta, 0, 3.13),
contours=[rts - 0.05, rts, rts + 0.05],
fill=False, labels=True, label_colors='black',
label_inline=True, label_inline_spacing=10)
graph_trap += tsp.plot(chart=S, ambient_coords=(chi, eta), color='brown',
label=' ', size=60) \
+ text(r"$\mathscr{S}$", (0.63, 2.16), color='brown', fontsize=18)
graph_trap.show(frame=True, axes_pad=0, ticks=[pi/8,pi/8], tick_formatter=(pi,pi),
figsize=8, xmax=pi/4+0.005, ymin=1.4, ymax=2.5)
graph_trap.save('lem_OS_diag_int_zoom.pdf', frame=True, axes_pad=0, ticks=[pi/8,pi/8],
tick_formatter=(pi,pi), figsize=8, xmax=pi/4+0.005, ymin=1.4, ymax=2.5)
graph = S.plot(EFI, ambient_coords=(r, ti), style={eta: '-', chi: '--'},
number_values={eta: 2, chi: 9}, label_axes=False)
# The stellar surface:
graph += surf.plot(EFI, ambient_coords=(r, ti), fixed_coords={chi: chis},
thickness=3, label_axes=False)
# Constant tau lines:
for iso_tau in iso_taus:
graph += iso_tau.plot(chart=EFI, ambient_coords=(r, ti), color='red',
style='--', label_axes=False)
# The initial star:
graph += init.plot(chart=EFI, ambient_coords=(r, ti),
fixed_coords={chi: chis}, color='magenta',
thickness=3, label_axes=False)
graph.axes_labels([r'$r/m$', r'$\tilde{t}/m$'])
graph.show(figsize=10)
The collapse end point correspond to $\eta = \pi$; it is thus given by I((pi, 0))
and its $\tilde{t}$ coordinate is obtained via the chart EF
:
EF(I((pi, 0)))
ti_end = EF(I((pi, 0)))[0]
ti_end
Maximal value of $\tilde{t}$ for plots:
timax = 17
graph += line([(0, ti_end + i*(40 - ti_end)/80) for i in range(81)],
thickness=3, color='darkorange', marker='D', markersize=10)
graph.set_axes_range(xmin=0, xmax=10, ymin=0, ymax=timax)
graph.axes_labels([r'$r/m$', r'$\tilde{t}/m$'])
graph.show(frame=True, figsize=10, axes_pad=0)
graph += hor_int.plot(chart=EFI, ambient_coords=(r, ti), color='black',
thickness=3, label_axes=False)
hor_int(chis)
EF(hor_int(chis))
tiH = EF(hor_int(chis))[0]
tiH, n(tiH)
hor_ext = M.curve({EF: (ti, 2)}, (ti, tiH, 40))
graph += hor_ext.plot(chart=EF, ambient_coords=(r, ti), color='black',
thickness=3, label_axes=False)
Adding the inner trapping horizon:
graph += ITH.plot(EFI, ambient_coords=(r, ti), color='blue',
thickness=3, style=':', label_axes=False)
graph.set_axes_range(xmin=0, xmax=10, ymin=0, ymax=timax)
graph.show(figsize=10)
graph_delay = copy(graph)
def geod_ext_out_I(tis, rs, rm):
r"""
External outgoing radial null geodesic starting from
(ti, r) = (tis, rs) in region I
"""
return M.curve({EF: [tis + r - rs + 4*ln(abs((r - 2)/(rs - 2))), r]}, (r, rs, rm))
lamb = var('lamb', latex_name=r'\lambda')
def geod_ext_out_II(tis, rs, rm):
r"""
External outgoing radial null geodesic starting from
(ti, r) = (tis, rs) in region II
"""
return M.curve({EF: [tis - lamb - rs + 4*ln(abs((-lamb - 2)/(rs - 2))), -lamb]}, (lamb, -rs, -rm))
def geod_ext_in(tis, rs, rm):
r"""
External ingoing radial null geodesic arriving at
(ti, r) = (tis, rs)
"""
return M.curve({EF: [tis - r + rs, r]}, (r, rs, rm))
def plot_geod_out(etas, chart=EF, plot_tangent=True, plot_int=True,
color='green', thickness=1.5):
r"""
Plot of an outgoing radial null geodesic
"""
# Internal part of the geodesic
gint = geod_int_out(etas)
# External part
tis, rs = EF(gint(chis))
if rs < 2:
rm = 0
gext = geod_ext_out_II(tis, rs, rm)
pv = -0.5*rs
else:
rm = 20
gext = geod_ext_out_I(tis, rs, rm)
pv = 1.4*rs
ambc = (chart[1], chart[0])
graph = gext.plot(chart, ambient_coords=ambc, color=color,
thickness=thickness, label_axes=False)
if plot_int:
graph += gint.plot(chart.restrict(I), ambient_coords=ambc, color=color,
thickness=thickness, label_axes=False)
if plot_tangent:
try:
vtan = gext.tangent_vector_field().at(gext.domain()(n(pv)))
graph += vtan.plot(chart, ambient_coords=ambc, color=color,
scale=0.01, arrowsize=4)
except ValueError:
pass
return graph
for etas in eta_sel:
graph += plot_geod_out(etas)
graph.axes_labels([r'$r/m$', r'$\tilde{t}/m$'])
graph.set_axes_range(xmin=0, xmax=10, ymin=0, ymax=timax)
graph.show(figsize=10)
def plot_geod_in(etas, chart=EF, color='green', thickness=1):
r"""
Plot of an ingoing radial null geodesic
"""
# Internal part of the geodesic
gint = geod_int_in(etas)
# External part
tis, rs = EF(gint(chis))
rm = 20
gext = geod_ext_in(tis, rs, rm)
ambc = (chart[1], chart[0])
return gint.plot(chart.restrict(I), ambient_coords=ambc,
color=color, thickness=thickness, style='--',
label_axes=False) \
+ gext.plot(chart, ambient_coords=ambc, color=color,
thickness=thickness, style='--',
label_axes=False)
for etas in eta_sel:
graph += plot_geod_in(etas)
for ti0 in [14, 18, 22, 26]:
graph += geod_ext_in(ti0, 0, 20).plot(EF, ambient_coords=(r, ti),
color='green', thickness=1,
style='--', label_axes=False)
graph_zoom = copy(graph)
graph += text(r'$\mathscr{H}$', (2.5, 16.4), fontsize=20, color='black')
graph.axes_labels([r'$r/m$', r'$\tilde{t}/m$'])
graph.set_axes_range(xmin=0, xmax=10, ymin=-0.03, ymax=timax)
graph.show(frame=True, axes=False, axes_pad=0, figsize=10)
graph.save('lem_OS_diag_EF.pdf', frame=True, axes=False,
axes_pad=0, figsize=10)
graph_zoom += plot_geod_out(3*pi/4) + plot_geod_in(5*pi/8)
graph_zoom += tsp.plot(chart=EFI, ambient_coords=(r, ti), color='brown',
label=' ', size=50) \
+ text(r"$\mathscr{S}$", (0.8, 11.6), color='brown', fontsize=18)
graph_zoom += line([(0, ti_end + i*(15 - ti_end)/20) for i in range(21)],
thickness=3, color='darkorange', marker='D', markersize=10)
graph_zoom.axes_labels([r'$r/m$', r'$\tilde{t}/m$'])
graph_zoom.show(frame=True, figsize=8, axes_pad=0,
xmin=0, xmax=3, ymin=9, ymax=13)
graph_zoom.save('lem_OS_diag_EF_zoom.pdf',frame=True,
figsize=8, axes_pad=0, xmin=0, xmax=3,
ymin=9, ymax=13)
for etas in eta_sel:
graph_delay += plot_geod_out(etas, plot_int=False)
graph_delay += line([(10, 0), (10, 29)], thickness=3, color='maroon')
graph_delay += text(r'$\mathscr{O}$', (9.1, 3), fontsize=24, color='maroon')
graph_delay += text(r'$\mathscr{H}$', (3.1, 27), fontsize=24, color='black')
graph_delay.axes_labels([r'$r/m$', r'$\tilde{t}/m$'])
graph_delay.set_axes_range(xmin=0, xmax=10.1, ymin=-0.03, ymax=29)
graph_delay.set_aspect_ratio(1)
graph_delay.save("lem_OS_diag_delay.pdf", frame=True, axes=False, figsize=10, axes_pad=0)
graph_delay.show(frame=True, axes=False, figsize=10, axes_pad=0)
KS.<T,X> = M.chart()
KS
#t1 = 6
t1 = 5
s = (ti - t1)/4
EF_to_KS = EF.transition_map(KS, [exp(r/4)*((exp(s) + exp(-s))/2 - r/4*exp(-s)),
exp(r/4)*((exp(s) - exp(-s))/2 + r/4*exp(-s))])
EF_to_KS.display()
graph = EF.plot(KS, ambient_coords=(X,T), color='grey',
style={ti: '-', r: '--'},
ranges={ti: (0, 16), r: (0, 8)},
number_values={ti: 9, r: 9})
graph.show(xmin=-5, xmax=6, ymin=-6, ymax=6)
EF_to_KSI = EF_to_KS.restrict(I)
EF_to_KSI.display()
I.atlas()
KSI = KS.restrict(I)
KSI is I.atlas()[2]
S_to_EF.display()
I.set_simplify_function(simplify)
S_to_KS = EF_to_KSI * S_to_EF
S_to_KS.display()
graph += S.plot(KSI, ambient_coords=(X, T), style={eta: '-', chi: '--'},
number_values={eta: 2, chi: 9}, label_axes=False)
# The stellar surface:
graph += surf.plot(chart=KSI, ambient_coords=(X, T), color='red',
thickness=3, label_axes=False)
# Constant tau lines:
for iso_tau in iso_taus:
graph += iso_tau.plot(chart=KSI, ambient_coords=(X, T), color='red',
style='--', label_axes=False)
# The initial star:
graph += init.plot(chart=KSI, ambient_coords=(X, T), color='magenta',
thickness=3, label_axes=False)
graph.show(xmin=-5, xmax=6, ymin=-6, ymax=5)
Check that the matter singularity occurs at the same point in the $(T,X)$ plane:
Ts, Xs = map(numerical_approx, KS(matter_sing(0)))
Ts, Xs
Ts, Xs = map(numerical_approx, KS(matter_sing(chis/2)))
Ts, Xs
Ts, Xs = map(numerical_approx, KS(matter_sing(chis)))
Ts, Xs
Plot of the singularity in all spacetime:
Xsing = [Xs + i*(5 - Xs)/20 for i in range(21)]
sing_plot = line([(X1, sqrt(1 + X1^2)) for X1 in Xsing],
color='darkorange', thickness=3,
marker='s', markersize=7)
graph += sing_plot
Adding the horizon:
graph += hor_int.plot(chart=KSI, ambient_coords=(X, T), color='black',
thickness=3, label_axes=False)
graph += hor_ext.plot(chart=KS, ambient_coords=(X, T), color='black',
thickness=3, label_axes=False)
graph.show(xmin=-3, xmax=6, ymin=-6, ymax=5)
for etas in eta_sel:
graph += plot_geod_out(etas, chart=KS, plot_tangent=False)
graph.show(xmin=-3, xmax=5, ymin=-5, ymax=5)
for etas in eta_sel:
graph += plot_geod_in(etas, chart=KS)
graph += geod_ext_in(14, 0, 20).plot(chart=KS, ambient_coords=(X, T),
color='green', thickness=1,
style='--', label_axes=False)
graph.show(xmin=-3, xmax=5, ymin=-5, ymax=5)
graph1 = copy(graph)
graph += text(r'$\tilde{t} = 0$', (5.2, -4.2), color='grey', rotation=-45, fontsize=14)
graph += text(r'$\tilde{t} = 4m$', (5.6, -2.4), color='grey', rotation=-35, fontsize=14)
graph += text(r'$\tilde{t} = 6m$', (5.6, -0.1), color='grey', rotation=-20, fontsize=14)
graph += text(r'$\tilde{t}\!=\!8m$', (5.65, 2.72), color='grey', rotation=7, fontsize=14)
graph += text(r'$r=0$', (1, 1.6), color='grey', rotation=35, fontsize=14)
graph += text(r'$r=3m$', (2.95, 2.7), color='grey', rotation=50, fontsize=14)
graph += text(r'$r=4m$', (3.28, -1.5), color='grey', rotation=-60, fontsize=14)
graph += text(r'$r=5m$', (5.52, -3.15), color='grey', rotation=-60, fontsize=14)
graph.show(xmin=-2, xmax=6, ymin=-4.5, ymax=4.5, figsize=10,
frame=True)
graph.save('lem_OS_diag_KS.pdf', xmin=-2, xmax=6, ymin=-4.5,
ymax=4.5, figsize=10, frame=True)
graph1.show(xmin=0., xmax=3.8, ymin=0.9, ymax=3.8, aspect_ratio=1,
frame=True)
graph1 += text(r'$\tilde{t}=8m$', (3.1, 2.22), color='grey', rotation=15, fontsize=14)
graph1 += text(r'$\tilde{t}=10m$', (3.55, 3.42), color='grey', rotation=35, fontsize=14)
graph1 += text(r'$r = 4m$', (3.6, 2.47), color='grey', rotation=60, fontsize=14)
graph1 += text(r'$r = 3m$', (3.3, 3.02), color='grey', rotation=50, fontsize=14)
graph1 += text(r'$r = 0$', (2.5, 2.78), color='grey', rotation=44, fontsize=14)
graph1.show(xmin=1.7, xmax=3.7, ymin=2, ymax=3.7, aspect_ratio=1, frame=True)
graph1.save('lem_OS_diag_KS_zoom.pdf', xmin=1.7, xmax=3.7, ymin=2, ymax=3.7,
aspect_ratio=1, frame=True)