Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download

SageMath notebooks associated to the Black Hole Lectures (https://luth.obspm.fr/~luthier/gourgoulhon/bh16)

Project: BHLectures
Views: 20094
Kernel: SageMath 10.0.beta8

Oppenheimer-Snyder collapse

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 (θ,φ)=const(\theta,\varphi) = \mathrm{const}:

M = Manifold(2, 'M', structure='Lorentzian') EF.<ti, r> = M.chart(r"ti:\tilde{t} r:[0,+oo)") EF

(M,(t~,r))\displaystyle \left(M,({\tilde{t}}, r)\right)

EF.coord_range()

t~: (,+);r: [0,+)\displaystyle {\tilde{t}} :\ \left( -\infty, +\infty \right) ;\quad r :\ \left[ 0 , +\infty \right)

Interior

I = M.open_subset('I')

We choose χs=π/4\chi_{\rm s} = \pi/4:

S.<eta,chi> = I.chart(r"eta:[0,pi]:\eta chi:[0,pi/4]:\chi") S

(I,(η,χ))\displaystyle \left(I,({\eta}, {\chi})\right)

S.coord_range()

η: [0,π];χ: [0,14π]\displaystyle {\eta} :\ \left[ 0 , \pi \right] ;\quad {\chi} :\ \left[ 0 , \frac{1}{4} \, \pi \right]

S.coord_bounds()

(((0,True),(π,True)),((0,True),(14π,True)))\displaystyle \left(\left(\left(0, \mathrm{True}\right), \left(\pi, \mathrm{True}\right)\right), \left(\left(0, \mathrm{True}\right), \left(\frac{1}{4} \, \pi, \mathrm{True}\right)\right)\right)

chis = S.coord_bounds()[1][1][0] chis

14π\displaystyle \frac{1}{4} \, \pi

The initial areal radius of the star in units of mm:

r0 = 2/sin(chis)^2 r0

4\displaystyle 4

Global evolution curves

rs(η)r_{\rm s}(\eta) in units of r0r_0:

rsurf(eta) = (1 + cos(eta))/2 rsurf

η  12cos(η)+12\displaystyle {\eta} \ {\mapsto}\ \frac{1}{2} \, \cos\left({\eta}\right) + \frac{1}{2}

Matter proper time τ\tau in units of r03/(8m0)\sqrt{r_0^3/(8 m_0)}

tau(eta) = (eta + sin(eta))/2 tau

η  12η+12sin(η)\displaystyle {\eta} \ {\mapsto}\ \frac{1}{2} \, {\eta} + \frac{1}{2} \, \sin\left({\eta}\right)

Proper energy density ρ\rho in units of 6m/(πr03)6 m / (\pi r_0^3):

rho(eta) = 1 / (1 + cos(eta))^3 rho

η  1(cos(η)+1)3\displaystyle {\eta} \ {\mapsto}\ \frac{1}{{\left(\cos\left({\eta}\right) + 1\right)}^{3}}

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)
Image in a Jupyter notebook
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)
Image in a Jupyter notebook

η\eta as a function of τ/m\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 mm:

tau_end = pi/sin(chis)^3 tau_end, n(tau_end)

(22π,8.88576587631673)\displaystyle \left(2 \, \sqrt{2} \pi, 8.88576587631673\right)

Check:

etaf(tau_end) == n(pi)

True\displaystyle \mathrm{True}

plot(etaf, 0, tau_end)
Image in a Jupyter notebook

Plot of ρ\rho at the black hole's birth as a function of χs\chi_s

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)
Image in a Jupyter notebook

Plot of the interior region in the (η,χ)(\eta,\chi) plane

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

[0.000000000000000,1.11072073453959,2.22144146907918,3.33216220361878,4.44288293815837,5.55360367269796,6.66432440723755,7.77504514177714,8.88576587631673]\displaystyle \left[0.000000000000000, 1.11072073453959, 2.22144146907918, 3.33216220361878, 4.44288293815837, 5.55360367269796, 6.66432440723755, 7.77504514177714, 8.88576587631673\right]

eta_sel = [etaf(tau0) for tau0 in tau_sel] eta_sel

[0.0,0.1969852777521357,0.39790775618592655,0.6073795752895496,0.8317111935799412,1.0810453707259005,1.3752523669869274,1.7683424316081915,3.141592653589793]\displaystyle \left[0.0, 0.1969852777521357, 0.39790775618592655, 0.6073795752895496, 0.8317111935799412, 1.0810453707259005, 1.3752523669869274, 1.7683424316081915, 3.141592653589793\right]

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 η=0\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)
Image in a Jupyter notebook

Null radial geodesics in the interior

Functions initializing the internal null geodesics

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

The event horizon

hor_int = geod_int_out(pi - 2*chis) graph += hor_int.plot(chart=S, ambient_coords=(chi, eta), color='black', thickness=3)

The final singularity

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)
Image in a Jupyter notebook

Eddington-Finkelstein-type coordinates for the interior

EFI = EF.restrict(I) I.atlas()

[(I,(η,χ)),(I,(t~,r))]\displaystyle \left[\left(I,({\eta}, {\chi})\right), \left(I,({\tilde{t}}, r)\right)\right]

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

{t~=4η+4log(cos(12η)+sin(12η))+2sin(η)r=22(cos(η)+1)sin(χ)\displaystyle \left\{\begin{array}{lcl} {\tilde{t}} & = & 4 \, {\eta} + 4 \, \log\left(\cos\left(\frac{1}{2} \, {\eta}\right) + \sin\left(\frac{1}{2} \, {\eta}\right)\right) + 2 \, \sin\left({\eta}\right) \\ r & = & 2 \, \sqrt{2} {\left(\cos\left({\eta}\right) + 1\right)} \sin\left({\chi}\right) \end{array}\right.

S_to_EF(eta, chi)

(4η+4log(cos(12η)+sin(12η))+2sin(η),22cos(η)sin(χ)+22sin(χ))\displaystyle \left(4 \, {\eta} + 4 \, \log\left(\cos\left(\frac{1}{2} \, {\eta}\right) + \sin\left(\frac{1}{2} \, {\eta}\right)\right) + 2 \, \sin\left({\eta}\right), 2 \, \sqrt{2} \cos\left({\eta}\right) \sin\left({\chi}\right) + 2 \, \sqrt{2} \sin\left({\chi}\right)\right)

rr(eta, chi) = S_to_EF(eta, chi)[1] rr

(η,χ)  22cos(η)sin(χ)+22sin(χ)\displaystyle \left( {\eta}, {\chi} \right) \ {\mapsto} \ 2 \, \sqrt{2} \cos\left({\eta}\right) \sin\left({\chi}\right) + 2 \, \sqrt{2} \sin\left({\chi}\right)

Adding isocontours r=constr = \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)
Image in a Jupyter notebook
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)
Image in a Jupyter notebook

Trapped surfaces

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)
Image in a Jupyter notebook

Adding rr isocontours near the trapped surface at (η,χ)=(11π6,3π16)(\eta,\chi) = \left(\frac{11\pi}{6}, \frac{3\pi}{16} \right):

tsp = I((11*pi/16, 3*pi/16), chart=S)
S(tsp)

(1116π,316π)\displaystyle \left(\frac{11}{16} \, \pi, \frac{3}{16} \, \pi\right)

EF(tsp)

(114π+4log(cos(1132π)+sin(1132π))+2sin(516π),22cos(516π)sin(316π)+22sin(316π))\displaystyle \left(\frac{11}{4} \, \pi + 4 \, \log\left(\cos\left(\frac{11}{32} \, \pi\right) + \sin\left(\frac{11}{32} \, \pi\right)\right) + 2 \, \sin\left(\frac{5}{16} \, \pi\right), -2 \, \sqrt{2} \cos\left(\frac{5}{16} \, \pi\right) \sin\left(\frac{3}{16} \, \pi\right) + 2 \, \sqrt{2} \sin\left(\frac{3}{16} \, \pi\right)\right)

rts = n(EF(tsp)[1]) rts

0.698372454547306\displaystyle 0.698372454547306

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)
Image in a Jupyter notebook

Plots in terms of Eddington-Finkelstein coordinates

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)
Image in a Jupyter notebook

Coordinate t~\tilde{t} of the collapse end point

The collapse end point correspond to η=π\eta = \pi; it is thus given by I((pi, 0)) and its t~\tilde{t} coordinate is obtained via the chart EF:

EF(I((pi, 0)))

(4π,0)\displaystyle \left(4 \, \pi, 0\right)

ti_end = EF(I((pi, 0)))[0] ti_end

4π\displaystyle 4 \, \pi

Drawing the singularity

Maximal value of t~\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)
Image in a Jupyter notebook

The event horizon

graph += hor_int.plot(chart=EFI, ambient_coords=(r, ti), color='black', thickness=3, label_axes=False)
hor_int(chis)

Point on the 2-dimensional Lorentzian manifold M\displaystyle \text{Point on the 2-dimensional Lorentzian manifold M}

EF(hor_int(chis))

(2π+2log(2)+2,2)\displaystyle \left(2 \, \pi + 2 \, \log\left(2\right) + 2, 2\right)

tiH = EF(hor_int(chis))[0] tiH, n(tiH)

(2π+2log(2)+2,9.66947966829948)\displaystyle \left(2 \, \pi + 2 \, \log\left(2\right) + 2, 9.66947966829948\right)

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)
Image in a Jupyter notebook
graph_delay = copy(graph)

Functions initializing the external null geodesics

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)
Image in a Jupyter notebook
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)
Image in a Jupyter notebook
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)
Image in a Jupyter notebook

Plot for the time delay effect

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)
Image in a Jupyter notebook

Plots in Kruskal-Szekeres coordinates

KS.<T,X> = M.chart() KS

(M,(T,X))\displaystyle \left(M,(T, X)\right)

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

{T=14(re(14t~+54)2e(14t~54)2e(14t~+54))e(14r)X=14(re(14t~+54)+2e(14t~54)2e(14t~+54))e(14r)\displaystyle \left\{\begin{array}{lcl} T & = & -\frac{1}{4} \, {\left(r e^{\left(-\frac{1}{4} \, {\tilde{t}} + \frac{5}{4}\right)} - 2 \, e^{\left(\frac{1}{4} \, {\tilde{t}} - \frac{5}{4}\right)} - 2 \, e^{\left(-\frac{1}{4} \, {\tilde{t}} + \frac{5}{4}\right)}\right)} e^{\left(\frac{1}{4} \, r\right)} \\ X & = & \frac{1}{4} \, {\left(r e^{\left(-\frac{1}{4} \, {\tilde{t}} + \frac{5}{4}\right)} + 2 \, e^{\left(\frac{1}{4} \, {\tilde{t}} - \frac{5}{4}\right)} - 2 \, e^{\left(-\frac{1}{4} \, {\tilde{t}} + \frac{5}{4}\right)}\right)} e^{\left(\frac{1}{4} \, r\right)} \end{array}\right.

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)
Image in a Jupyter notebook
EF_to_KSI = EF_to_KS.restrict(I) EF_to_KSI.display()

{T=14(re(14t~+54)2e(14t~54)2e(14t~+54))e(14r)X=14(re(14t~+54)+2e(14t~54)2e(14t~+54))e(14r)\displaystyle \left\{\begin{array}{lcl} T & = & -\frac{1}{4} \, {\left(r e^{\left(-\frac{1}{4} \, {\tilde{t}} + \frac{5}{4}\right)} - 2 \, e^{\left(\frac{1}{4} \, {\tilde{t}} - \frac{5}{4}\right)} - 2 \, e^{\left(-\frac{1}{4} \, {\tilde{t}} + \frac{5}{4}\right)}\right)} e^{\left(\frac{1}{4} \, r\right)} \\ X & = & \frac{1}{4} \, {\left(r e^{\left(-\frac{1}{4} \, {\tilde{t}} + \frac{5}{4}\right)} + 2 \, e^{\left(\frac{1}{4} \, {\tilde{t}} - \frac{5}{4}\right)} - 2 \, e^{\left(-\frac{1}{4} \, {\tilde{t}} + \frac{5}{4}\right)}\right)} e^{\left(\frac{1}{4} \, r\right)} \end{array}\right.

I.atlas()

[(I,(η,χ)),(I,(t~,r)),(I,(T,X))]\displaystyle \left[\left(I,({\eta}, {\chi})\right), \left(I,({\tilde{t}}, r)\right), \left(I,(T, X)\right)\right]

KSI = KS.restrict(I) KSI is I.atlas()[2]

True\displaystyle \mathrm{True}

S_to_EF.display()

{t~=4η+4log(cos(12η)+sin(12η))+2sin(η)r=22(cos(η)+1)sin(χ)\displaystyle \left\{\begin{array}{lcl} {\tilde{t}} & = & 4 \, {\eta} + 4 \, \log\left(\cos\left(\frac{1}{2} \, {\eta}\right) + \sin\left(\frac{1}{2} \, {\eta}\right)\right) + 2 \, \sin\left({\eta}\right) \\ r & = & 2 \, \sqrt{2} {\left(\cos\left({\eta}\right) + 1\right)} \sin\left({\chi}\right) \end{array}\right.

I.set_simplify_function(simplify) S_to_KS = EF_to_KSI * S_to_EF S_to_KS.display()

{T=12(2(cos(η)+1)e(η12sin(η)+54)sin(χ)cos(12η)+sin(12η)(cos(12η)+sin(12η))e(η+12sin(η)54)e(η12sin(η)+54)cos(12η)+sin(12η))e(122(cos(η)+1)sin(χ))X=12(2(cos(η)+1)e(η12sin(η)+54)sin(χ)cos(12η)+sin(12η)+(cos(12η)+sin(12η))e(η+12sin(η)54)e(η12sin(η)+54)cos(12η)+sin(12η))e(122(cos(η)+1)sin(χ))\displaystyle \left\{\begin{array}{lcl} T & = & -\frac{1}{2} \, {\left(\frac{\sqrt{2} {\left(\cos\left({\eta}\right) + 1\right)} e^{\left(-{\eta} - \frac{1}{2} \, \sin\left({\eta}\right) + \frac{5}{4}\right)} \sin\left({\chi}\right)}{\cos\left(\frac{1}{2} \, {\eta}\right) + \sin\left(\frac{1}{2} \, {\eta}\right)} - {\left(\cos\left(\frac{1}{2} \, {\eta}\right) + \sin\left(\frac{1}{2} \, {\eta}\right)\right)} e^{\left({\eta} + \frac{1}{2} \, \sin\left({\eta}\right) - \frac{5}{4}\right)} - \frac{e^{\left(-{\eta} - \frac{1}{2} \, \sin\left({\eta}\right) + \frac{5}{4}\right)}}{\cos\left(\frac{1}{2} \, {\eta}\right) + \sin\left(\frac{1}{2} \, {\eta}\right)}\right)} e^{\left(\frac{1}{2} \, \sqrt{2} {\left(\cos\left({\eta}\right) + 1\right)} \sin\left({\chi}\right)\right)} \\ X & = & \frac{1}{2} \, {\left(\frac{\sqrt{2} {\left(\cos\left({\eta}\right) + 1\right)} e^{\left(-{\eta} - \frac{1}{2} \, \sin\left({\eta}\right) + \frac{5}{4}\right)} \sin\left({\chi}\right)}{\cos\left(\frac{1}{2} \, {\eta}\right) + \sin\left(\frac{1}{2} \, {\eta}\right)} + {\left(\cos\left(\frac{1}{2} \, {\eta}\right) + \sin\left(\frac{1}{2} \, {\eta}\right)\right)} e^{\left({\eta} + \frac{1}{2} \, \sin\left({\eta}\right) - \frac{5}{4}\right)} - \frac{e^{\left(-{\eta} - \frac{1}{2} \, \sin\left({\eta}\right) + \frac{5}{4}\right)}}{\cos\left(\frac{1}{2} \, {\eta}\right) + \sin\left(\frac{1}{2} \, {\eta}\right)}\right)} e^{\left(\frac{1}{2} \, \sqrt{2} {\left(\cos\left({\eta}\right) + 1\right)} \sin\left({\chi}\right)\right)} \end{array}\right.

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)
Image in a Jupyter notebook

Check that the matter singularity occurs at the same point in the (T,X)(T,X) plane:

Ts, Xs = map(numerical_approx, KS(matter_sing(0))) Ts, Xs

(3.39037541861741,3.23954402334112)\displaystyle \left(3.39037541861741, 3.23954402334112\right)

Ts, Xs = map(numerical_approx, KS(matter_sing(chis/2))) Ts, Xs

(3.39037541861741,3.23954402334112)\displaystyle \left(3.39037541861741, 3.23954402334112\right)

Ts, Xs = map(numerical_approx, KS(matter_sing(chis))) Ts, Xs

(3.39037541861741,3.23954402334112)\displaystyle \left(3.39037541861741, 3.23954402334112\right)

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)
Image in a Jupyter notebook
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)
Image in a Jupyter notebook
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)
Image in a Jupyter notebook
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)
Image in a Jupyter notebook
graph1.show(xmin=0., xmax=3.8, ymin=0.9, ymax=3.8, aspect_ratio=1, frame=True)
Image in a Jupyter notebook
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)
Image in a Jupyter notebook