Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download
Project: BHLectures
Views: 628
Kernel: SageMath 9.2.beta12

Radial null geodesics in Schwarzschild spacetime

This Jupyter/SageMath worksheet is relative to the lectures Geometry and physics of black holes

These computations are based on SageManifolds (version 1.0, as included in SageMath 7.5 and higher versions)

Click here to download the worksheet file (ipynb format). To run it, you must start SageMath with the Jupyter notebook, with the command sage -n jupyter

NB: a version of SageMath at least equal to 7.5 is required to run this worksheet:

version()
'SageMath version 9.2.beta12, Release Date: 2020-09-06'

First we set up the notebook to display mathematical objects using LaTeX formatting:

%display latex

Spacetime

We declare the spacetime manifold MM:

M = Manifold(4, 'M') print(M)
4-dimensional differentiable manifold M

The Schwarzschild-Droste coordinates (t,r,θ,ϕ)(t,r,\theta,\phi):

X.<t,r,th,ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi') X
(M,(t,r,θ,ϕ))\renewcommand{\Bold}[1]{\mathbf{#1}}\left(M,(t, r, {\theta}, {\phi})\right)

The Schwarzschild metric:

g = M.lorentzian_metric('g') m = var('m') ; assume(m>=0) g[0,0], g[1,1] = -(1-2*m/r), 1/(1-2*m/r) g[2,2], g[3,3] = r^2, (r*sin(th))^2 g.display()
g=(2mr1)dtdt+(12mr1)drdr+r2dθdθ+r2sin(θ)2dϕdϕ\renewcommand{\Bold}[1]{\mathbf{#1}}g = \left( \frac{2 \, m}{r} - 1 \right) \mathrm{d} t\otimes \mathrm{d} t + \left( -\frac{1}{\frac{2 \, m}{r} - 1} \right) \mathrm{d} r\otimes \mathrm{d} r + r^{2} \mathrm{d} {\theta}\otimes \mathrm{d} {\theta} + r^{2} \sin\left({\theta}\right)^{2} \mathrm{d} {\phi}\otimes \mathrm{d} {\phi}

Radial null geodesics

The outgoing family:

var('u') outgeod = M.curve({X: [r + 2*m*ln(abs(r/(2*m)-1)) + u, r, pi/2, pi]}, (r, 0, +Infinity)) outgeod.display()
(0,+)Mr(t,r,θ,ϕ)=(2mlog(r2m1)+r+u,r,12π,π)\renewcommand{\Bold}[1]{\mathbf{#1}}\begin{array}{llcl} & \left(0, +\infty\right) & \longrightarrow & M \\ & r & \longmapsto & \left(t, r, {\theta}, {\phi}\right) = \left(2 \, m \log\left({\left| \frac{r}{2 \, m} - 1 \right|}\right) + r + u, r, \frac{1}{2} \, \pi, \pi\right) \end{array}

The ingoing family:

var('v') ingeod = M.curve({X: [-r - 2*m*ln(abs(r/(2*m)-1)) + v, r, pi/2, pi]}, (r, 0, +Infinity)) ingeod.display()
(0,+)Mr(t,r,θ,ϕ)=(2mlog(r2m1)r+v,r,12π,π)\renewcommand{\Bold}[1]{\mathbf{#1}}\begin{array}{llcl} & \left(0, +\infty\right) & \longrightarrow & M \\ & r & \longmapsto & \left(t, r, {\theta}, {\phi}\right) = \left(-2 \, m \log\left({\left| \frac{r}{2 \, m} - 1 \right|}\right) - r + v, r, \frac{1}{2} \, \pi, \pi\right) \end{array}
graph = Graphics() for u0 in range(-20, 20, 2): graph += outgeod.plot(ambient_coords=(r,t), prange=(0.01, 1.99), parameters={m: 1, u: u0}, color='green', style='-', thickness=1, label_axes=False) graph += outgeod.plot(ambient_coords=(r,t), prange=(2.01, 8), parameters={m: 1, u: u0}, color='green', style='-', thickness=1, label_axes=False) graph += ingeod.plot(ambient_coords=(r,t), prange=(0.01, 1.99), parameters={m: 1, v: u0}, color='green', style='--', thickness=1, label_axes=False) graph += ingeod.plot(ambient_coords=(r,t), prange=(2.01, 8), parameters={m: 1, v: u0}, color='green', style='--', thickness=1, label_axes=False) show(graph, axes_labels=[r"$r/m$", r"$t/m$"], aspect_ratio=1, ymin=-4, ymax=4)
Image in a Jupyter notebook

Plot of some light cones:

fout(r,u) = r + 2*ln(abs(r/2-1) + 1e-10) + u fin(r,v) = -r - 2*ln(abs(r/2-1) + 1e-10) + v def cone1(r,t): return t<fout(r,2) and t>fin(r,0) and r>1 def cone2(r,t): return t>fout(r,2) and t>fin(r,0) and t<1.8 def cone3(r,t): return t>fout(r,-2) and t>fin(r,4) and t<1.8 def cone4(r,t): return t>fout(r,-6) and t>fin(r,8) and t<1.8 graph += region_plot(cone1, (1,2), (0,2), incol='yellow') graph += region_plot(cone2, (2,8), (0,2), incol='yellow') graph += region_plot(cone3, (2,8), (0,2), incol='yellow') graph += region_plot(cone4, (2,8), (0,2), incol='yellow') show(graph, axes_labels=[r"$r/m$", r"$t/m$"], aspect_ratio=1, ymin=-4, ymax=4)
Image in a Jupyter notebook
graph.save("sch_rad_null_geod.pdf", aspect_ratio=1, ymin=-4, ymax=4)

Eddington-Finkelstein coordinates

The ingoing Eddington-Finkelstein chart:

X_EF.<te,r,th,ph> = M.chart(r'te:\tilde{t} r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi') X_EF
(M,(t~,r,θ,ϕ))\renewcommand{\Bold}[1]{\mathbf{#1}}\left(M,({\tilde{t}}, r, {\theta}, {\phi})\right)
X_to_X_EF = X.transition_map(X_EF, [t+2*m*ln(abs(r/(2*m)-1)), r, th, ph]) X_to_X_EF.display()
{t~=2mlog(r2m1)+tr=rθ=θϕ=ϕ\renewcommand{\Bold}[1]{\mathbf{#1}}\left\{\begin{array}{lcl} {\tilde{t}} & = & 2 \, m \log\left({\left| \frac{r}{2 \, m} - 1 \right|}\right) + t \\ r & = & r \\ {\theta} & = & {\theta} \\ {\phi} & = & {\phi} \end{array}\right.
X_to_X_EF.inverse().display()
{t=2mlog(2)+2mlog(m)2mlog(2m+r)+t~r=rθ=θϕ=ϕ\renewcommand{\Bold}[1]{\mathbf{#1}}\left\{\begin{array}{lcl} t & = & 2 \, m \log\left(2\right) + 2 \, m \log\left(m\right) - 2 \, m \log\left({\left| -2 \, m + r \right|}\right) + {\tilde{t}} \\ r & = & r \\ {\theta} & = & {\theta} \\ {\phi} & = & {\phi} \end{array}\right.
ingeod.display()
(0,+)Mr(t,r,θ,ϕ)=(2mlog(r2m1)r+v,r,12π,π)r(t~,r,θ,ϕ)=(r+v,r,12π,π)\renewcommand{\Bold}[1]{\mathbf{#1}}\begin{array}{llcl} & \left(0, +\infty\right) & \longrightarrow & M \\ & r & \longmapsto & \left(t, r, {\theta}, {\phi}\right) = \left(-2 \, m \log\left({\left| \frac{r}{2 \, m} - 1 \right|}\right) - r + v, r, \frac{1}{2} \, \pi, \pi\right) \\ & r & \longmapsto & \left({\tilde{t}}, r, {\theta}, {\phi}\right) = \left(-r + v, r, \frac{1}{2} \, \pi, \pi\right) \end{array}
outgeod.display()
(0,+)Mr(t,r,θ,ϕ)=(2mlog(r2m1)+r+u,r,12π,π)r(t~,r,θ,ϕ)=(4mlog(2)4mlog(m)+4mlog(2m+r)+r+u,r,12π,π)\renewcommand{\Bold}[1]{\mathbf{#1}}\begin{array}{llcl} & \left(0, +\infty\right) & \longrightarrow & M \\ & r & \longmapsto & \left(t, r, {\theta}, {\phi}\right) = \left(2 \, m \log\left({\left| \frac{r}{2 \, m} - 1 \right|}\right) + r + u, r, \frac{1}{2} \, \pi, \pi\right) \\ & r & \longmapsto & \left({\tilde{t}}, r, {\theta}, {\phi}\right) = \left(-4 \, m \log\left(2\right) - 4 \, m \log\left(m\right) + 4 \, m \log\left({\left| -2 \, m + r \right|}\right) + r + u, r, \frac{1}{2} \, \pi, \pi\right) \end{array}

Plot of the radial null geodesics in terms of the ingoing Eddington-Finkelstein coordinates:

graph = Graphics() for u0 in range(-20, 20, 2): graph += outgeod.plot(chart=X_EF, ambient_coords=(r,te), prange=(0.01, 1.99), parameters={m: 1, u: u0}, color='green', style='-', thickness=1, label_axes=False) graph += outgeod.plot(chart=X_EF, ambient_coords=(r,te), prange=(2.01, 8), parameters={m: 1, u: u0}, color='green', style='-', thickness=1, label_axes=False) graph += ingeod.plot(chart=X_EF, ambient_coords=(r,te), prange=(0.01, 8), parameters={m: 1, v: u0}, color='green', style='--', thickness=1, label_axes=False) show(graph, axes_labels=[r"$r/m$", r"$\tilde{t}/m$"], aspect_ratio=1, ymin=-4, ymax=4)
Image in a Jupyter notebook
fout(r,u) = r + 4*ln(abs(r/2-1) + 1e-10) + u fin(r,v) = -r + v def cone1(r,t): return t<fout(r,4) and t>fin(r,2) and t<2 def cone2(r,t): return t>fout(r,2) and t>fin(r,0) and t<1.8 def cone3(r,t): return t>fout(r,-2) and t>fin(r,4) and t<1.5 def cone4(r,t): return t>fout(r,-10) and t>fin(r,8) and t<1.8 graph += region_plot(cone1, (1,2), (0,2), incol='yellow') #graph += region_plot(cone2, (2,8), (0,2), incol='yellow') graph += region_plot(cone3, (2,8), (0,2), incol='yellow') graph += region_plot(cone4, (2,8), (0,2), incol='yellow') show(graph, axes_labels=[r"$r/m$", r"$\tilde{t}/m$"], aspect_ratio=1, ymin=-4, ymax=4)
Image in a Jupyter notebook
graph.save("sch_rad_null_geod_EF.pdf", aspect_ratio=1, ymin=-4, ymax=4)

Plot of the hypersurfaces t=constt=\mathrm{const} in terms of the ingoing Eddington-Finkelstein coordinates:

graph = Graphics() for t0 in range(-8, 10): graph += X.plot(X_EF, ranges={r:(2.01,8)}, fixed_coords={t:t0,th:pi/2,ph:0}, ambient_coords=(r,te), style={t:'--', r:'-'}, parameters={m:1}, color='blue') graph += X.plot(X_EF, ranges={r:(0.01, 1.99)}, fixed_coords={t:t0,th:pi/2,ph:0}, ambient_coords=(r,te), style={t:'--', r:'-'}, parameters={m:1}, color='blue') show(graph, axes_labels=[r"$r/m$", r"$\tilde{t}/m$"], aspect_ratio=1, ymin=-4, ymax=4)
Image in a Jupyter notebook
graph.save("sch_SD_slices.pdf", axes_labels=[r"$r/m$", r"$\tilde{t}/m$"], aspect_ratio=1, ymin=-4, ymax=4)