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.3.beta1, Release Date: 2023-12-10'
First we set up the notebook to display mathematical objects using LaTeX formatting:
%display latex
We declare the spacetime manifold $M$:
M = Manifold(4, 'M', structure='Lorentzian')
print(M)
4-dimensional Lorentzian manifold M
We introduce coordinates $(v,r,\theta,\varphi)$ analogous to the ingoing null Eddington-Finkelstein coordinates in Schwarzschild spacetime, i.e. such that $v$ is constant along ingoing radial null geodesics:
XN.<v,r,th,ph> = M.chart(r'v r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\varphi:periodic')
XN
XN.coord_range()
The metric tensor corresponding to the Vaidya solution is:
m = function('m')
g = M.metric()
g[0,0] = -(1 - 2*m(v)/r)
g[0,1] = 1
g[2,2] = r^2
g[3,3] = (r*sin(th))^2
g.display()
g.inverse().display()
g.inverse()[:]
The Ricci tensor is
Ric = g.ricci()
Ric.display()
It has zero trace, i.e. the Ricci scalar vanishes:
g.ricci_scalar().expr()
The Riemann tensor:
Riem = g.riemann()
Riem.display_comp(only_nonredundant=True)
The Kretschmann scalar $K = R_{abcd} R^{abcd}$:
K = Riem.down(g)['_{abcd}'] * Riem.up(g)['^{abcd}']
K.expr()
XN.coframe()
dv = XN.coframe()[0]
dv.display()
k = - dv.up(g)
k.set_name('k')
k.display()
Check that 𝑘 is a null vector:
g(k, k).expr()
Check that $k$ is a geodesic vector field, i.e. fulfils $\nabla_k k = 0$:
nabla = g.connection()
acc = nabla(k).contract(k)
acc.display()
xi = M.vector_field(v, r, 0, 0, name='xi',
latex_name=r'\xi')
xi.display()
Lg = g.lie_derivative(xi) - 2*g
Lg.display()
If $m(v)$ is a linear function, the above result is identically zero, showing that $\mathcal{L}_{\xi} g = 2 g$ in that case, i.e. that $\xi$ is a homethetic Killing vector.
Let us introduce a new chart $(t,r,\theta,\varphi)$ such that the advanced time $t+r$ is $v$: $v = t + r$; this is the analog of ingoing Eddington-Finkelstein (IEF) coordinates in Schwarzschild spacetime.
X.<t,r,th,ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\varphi:periodic')
X
X.coord_range()
We declare the transition map between the $(t,r,\theta,\varphi)$ and $(v,r,\theta,\varphi)$ coordinates:
X_to_XN = X.transition_map(XN, (t + r, r, th, ph))
X_to_XN.display()
X_to_XN.inverse().display()
Expression of the metric tensor in the IEF coordinates:
g.display(X)
From now on, we set the IEF chart X
to be the default one on $M$:
M.set_default_chart(X)
M.set_default_frame(X.frame())
Then g.display(X)
can be substituted by g.display()
:
g.display()
g.inverse()[:]
xi.display()
The Ricci tensor in terms of the IEF coordinates:
Ric.display()
The notation $\frac{\partial m}{\partial(r+t)}$ to denote $\frac{\mathrm{d}m}{\mathrm{d}v}$ is quite unfortunate (this shall be improved in a future version). The display of the corresponding symbolic expression is slightly better, $\mathrm{D}_0(m)$ standing for the derivative of function $m$ with respect to its first (index $0$) and unique argument, i.e. $\mathrm{D}_0(m) = \frac{\mathrm{d}m}{\mathrm{d}v}$:
Ric[0,0].expr()
The Ricci scalar is vanishing:
g.ricci_scalar().display()
The energy-momentum vector ensuring that the Einstein equation is fulfilled is then:
T = 1/(8*pi)*Ric
T.set_name('T')
T.display()
Since $v=t+r$, we have $\mathrm{d}v = \mathrm{d}t + \mathrm{d}r$:
dv.display()
The derivative of the function $m(v)$:
mp(v) = diff(m(v), v)
mp(v)
T == 1/(4*pi)*mp(t+r)/r^2 * dv*dv
The future-directed null vector along the ingoing null geodesics:
k.display()
Let us consider the vector field:
l = M.vector_field(1, (r - 2*m(t+r))/(r+2*m(t+r)), 0, 0,
name='l', latex_name=r'\ell')
l.display()
It is a null vector:
g(l,l).display()
Moreover $\ell$ is a pregeodesic vector field, i.e. it obeys $\nabla_\ell \ell = \kappa \ell$:
acc = nabla(l).contract(l)
acc.display()
kappa = acc[0]/l[0]
kappa
acc == kappa*l
The outgoing radial null geodesics are the field lines of $\ell$; they thus obey to $$ \frac{\mathrm{d}r}{\mathrm{d}t} = \frac{\ell^r}{\ell^t}.$$ Hence the value of $\frac{\mathrm{d}r}{\mathrm{d}t}$:
drdt = (l[1] / l[0]).expr()
drdt
Let us choose a simple function $m(v)$, based on one of the following smoothstep functions:
S0(x) = x
S1(x) = -2*x^3 + 3*x^2
S2(x) = 6*x^5 - 15*x^4 + 10*x^3
S(x) = S0(x)
#S(x) = S2(x)
m_0 = var('m_0')
v_0 = var('v_0')
h(v) = (1+sgn(v))/2 # the Heaviside function
mS(v) = m_0*(h(v)*h(v_0 - v)*S(v/v_0) + h(v - v_0))
mS(v)
NB: we don't use Sage's predefined heaviside
function, since it is incompatible with SciPy numerical integrators.
The singularity is entirely hidden under the event horizon for large energy densities of the radiation shell, i.e. for small values of the shell width $v_0$. The precise criterion is $v_0 < 16 m$ for $S(x) = S_0(x)$. We select here $v_0 = 3m$:
v0 = 3
m_num(v) = mS(v).subs({m_0: 1, v_0: v0})
m_num(v)
graph = plot(m_num(v), (v, -v0, 2*v0), color='blue',
legend_label=r'$S(x) = x$', thickness=3,
axes_labels=[r'$v/m$', r'$M(v)/m$'],
frame=True, axes=False, gridlines=True)
graph += plot(S2(v/v0) , (v, 0, v0), color='red',
legend_label=r'$S(x) = 6x^5 - 15x^4 + 10x^3$',
thickness=3)
graph += line([(v0, 0), (v0, 1)], linestyle='--', color='grey')
graph += text(r'$v_0$', (v0, -0.08), color='black', fontsize=16)
graph.set_axes_range(ymin=0)
graph.save("vai_mass_function.pdf")
graph
We plug the function $m(v)$ into the expression of $\frac{\mathrm{d}r}{\mathrm{d}t}$ along the outgoing radial null geodesics found above:
drdt1 = drdt.substitute_function(m, m_num)
drdt1
and we perform a numerical integration:
fdrdt = fast_callable(drdt1, vars=[t, r])
from scipy.integrate import ode
def solve_ode(t0, t1, r0=0.001, dt=0.02):
t0 = RDF(t0)
t1 = RDF(t1)
r0 = RDF(r0)
dt = RDF(dt)
forward = 1 if t1 > t0 else -1
rd = ode(solve_ode.rhs)
rd.set_initial_value(r0, t0)
rd.set_integrator('dopri5')
sol = []
ts = t0
rs = r0
while rd.successful() and forward*ts < forward*t1 and rs > 0:
sol.append((ts, rs))
ts = rd.t + forward*dt
rs = RDF(rd.integrate(ts)[0])
return sol
solve_ode.rhs = fdrdt
Plot parameters:
tmax = v0 + 1.5
ymin = -4
ymax = tmax - 0.5
rmax = 6
Outgoing radial null geodesics from the Minkowski region:
def geod_line(sol):
return line([(s[1], s[0]) for s in sol], color='green')
def outgeods_from_Mink(t0, tmax, dt0, t0_max=0):
geods = []
print('t0 : ', end='')
while t0 < t0_max:
sol = solve_ode(t0, tmax, r0=1e-8, dt=0.05)
print(t0, end=', ')
geods.append(geod_line(sol))
t0 += dt0
return geods
outgeods = outgeods_from_Mink(-10., tmax, 1., t0_max=-3.)
t0 : -10.0000000000000, -9.00000000000000, -8.00000000000000, -7.00000000000000, -6.00000000000000, -5.00000000000000, -4.00000000000000,
dt0 = 0.5 if S == S0 else 0.49
outgeods += outgeods_from_Mink(-3., tmax, dt0)
t0 : -3.00000000000000, -2.50000000000000, -2.00000000000000, -1.50000000000000, -1.00000000000000, -0.500000000000000,
Drawing of the radiation region (yellow rays):
def draw_radiation_region(v0, rmax, nl):
graph = Graphics()
for i in range(nl+1):
t0 = float(i)/float(nl)*v0
graph += line([(0, t0), (rmax, t0 - rmax)], color='yellow')
return graph
graph = draw_radiation_region(v0, rmax, 20)
Adding the outgoing radial null geodesics:
for geod in outgeods:
graph += geod
show(graph, aspect_ratio=1, xmax=rmax, ymin=ymin, ymax=ymax,
axes_labels=[r'$r/m$', r'$t/m$'], figsize=10)
The ingoing null geodesics:
def ingoing_geods(tmin, tmax, rmax, v0, dt_out, dt_in):
geods = []
t0 = tmin
while t0 < 0:
geods.append(line([(0, t0), (rmax, t0 - rmax)], color='green',
linestyle='--'))
t0 += dt_out
t0 = 0
while t0 <= v0:
geods.append(line([(0, t0), (rmax, t0 - rmax)], color='green',
linestyle='--'))
t0 += dt_in
t0 = v0 + dt_out
while t0 < tmax + rmax:
geods.append(line([(0, t0), (rmax, t0 - rmax)], color='green',
linestyle='--'))
t0 += dt_out
return geodsd𝑟
ingeods = ingoing_geods(-6., tmax, rmax, v0, 1., 0.5)
for geod in ingeods:
graph += geod
def curvature_sing(tmax, nb):
h = tmax/(nb - 1)
return line([(0, h*i) for i in range(nb)],
thickness=3, color='darkorange',
marker='D', markersize=8)
graph += curvature_sing(tmax, 21)
def event_hor(tmax, tmin=-4.):
sol = solve_ode(tmax, tmin, r0=2.)
hor = line([(s[1], s[0]) for s in sol], color='black', thickness=4)
# Refinement near the horizon birth point:
t0, r0 = sol[-2][0], sol[-2][1]
sol = solve_ode(t0, tmin, r0=r0, dt=0.001)
hor += line([(s[1], s[0]) for s in sol], color='black', thickness=4)
thb = sol[-1][0]
dthb = sol[-2][0] - thb
print("Coordinate t at the event horizon birth [unit: m]: ")
print(thb, " +/- ", dthb)
return hor
graph += event_hor(tmax)
pA = (2., v0 - 2.) # point A
Coordinate t at the event horizon birth [unit: m]: -2.600999999999992 +/- 0.0009999999999998899
Check of the determination of $t_{\rm hb}$ by comparison with the analytic formula for $S(x) = S_0(x) := x$:
if S == S0:
aa = 1/sqrt(16/v0 - 1)
thb = numerical_approx(-4*exp(-2*aa*atan
(aa)))
print("Analytic value of coordinate t at the event horizon birth [unit: m]:")
print(thb)
pB = (-thb/2, thb/2) # point B
pC = (0, thb) # point C
Analytic value of coordinate t at the event horizon birth [unit: m]: -2.60135096372454
def trapping_hor(mv, tmax):
return parametric_plot((2*mv, v - 2*mv), (v, 0, tmax + 2),
color='red', thickness=2)
graph += trapping_hor(m_num(v), tmax)
graph_wo_vectors = copy(graph)
The vectors $k$ and $\ell$ at some point $p$:
p = M((1, 4, pi/2, 0), chart=X)
l.at(p).display()
graph += k.at(p).plot(ambient_coords=(r, t), color='green', fontsize=16)
graph += l.at(p).plot(ambient_coords=(r, t), color='green', fontsize=16,
parameters={m(5): 1}, label_offset=0.15)
if S == S0:
graph += circle(pA, 0.08, fill=True, color='black') \
+ text(r"$A$", (pA[0]+0.2, pA[1]+0.1), fontsize=16, color='black')
graph += circle(pB, 0.08, fill=True, color='black') \
+ text(r"$B$", (pB[0]+0.25, pB[1]), fontsize=16, color='black')
graph += circle(pC, 0.08, fill=True, color='black') \
+ text(r"$C$", (pC[0]+0.25, pC[1]), fontsize=16, color='black')
show(graph, aspect_ratio=1, xmax=rmax, ymin=ymin, ymax=ymax,
axes_labels=[r'$r/m$', r'$t/m$'], figsize=10)
filename = "vai_diag_S0.pdf" if S == S0 else "vai_diag_S2.pdf"
graph.save(filename, aspect_ratio=1, xmax=rmax, ymin=ymin,
ymax=ymax, axes_labels=[r'$r/m$', r'$t/m$'], figsize=10)
A zoom on the trapping horizon in its dynamical part: notice that the "outgoing" null geodesics cross it with a vertical tangent, in agreement with the cross-sections of the trapping horizon being marginally trapped surfaces.
show(graph_wo_vectors, aspect_ratio=1, ymin=-0.8, ymax=1.2, xmax=2.5,
axes_labels=[r'$r/m$', r'$t/m$'])
The naked singularity is obtained for small energy densities of the radiation shell, i.e. for large values of the shell width $v_0$. The precise criterion is $v_0 > 16 m$ for $S(x) = S_0(x)$. We select here $v_0 = 18m$:
v0 = 18
m_num(v) = mS(v).subs({m_0: 1, v_0: v0})
m_num(v)
drdt1 = drdt.substitute_function(m, m_num)
fdrdt = fast_callable(drdt1, vars=[t, r])
solve_ode.rhs = fdrdt
#tmax = v0 + 4.5
tmax = v0 + 3.5
ymax = tmax - 0.5
rmax = 16
Outgoing radial null geodesics from the Minkowski region:
outgeods = outgeods_from_Mink(-rmax - 3, tmax, 2.)
t0 : -19, -17.0000000000000, -15.0000000000000, -13.0000000000000, -11.0000000000000, -9.00000000000000, -7.00000000000000, -5.00000000000000, -3.00000000000000, -1.00000000000000,
Outgoing radial null geodesics in the black hole region:
def outgeods_in_BH(tmax):
geods = []
sol = solve_ode(tmax, -1, r0=1)
geods.append(geod_line(sol))
sol = solve_ode(2/3*tmax, -1, r0=1)
geods.append(geod_line(sol))
sol = solve_ode(2/3*tmax, tmax, r0=1)
geods.append(geod_line(sol))
sol = solve_ode(1/3*tmax, -1, r0=0.5)
geods.append(geod_line(sol))
sol = solve_ode(1/3*tmax, tmax, r0=0.5)
geods.append(geod_line(sol))
return geods
outgeods += outgeods_in_BH(tmax)
if S == S0:
a0 = 2 / v0
rC = n(4/(1 - sqrt(1 - 8*a0)))
tC = v0 - rC
pC = (rC, tC) # point C
sol = solve_ode(tC, -1, r0=rC)
cauchy_hor = line([(s[1], s[0]) for s in sol],
color='blue', thickness=2)
sol = solve_ode(tC, tmax, r0=rC)
cauchy_hor += line([(s[1], s[0]) for s in sol],
color='blue', thickness=2)
else:
print("approximate CH")
sol = solve_ode(-1.e-3, tmax, r0=1e-8)
cauchy_hor = line([(s[1], s[0]) for s in sol],
color='blue', thickness=2)
Outgoing radial null geodesics emerging from the initial singularity:
if S == S0:
rB = n(4/(1 + sqrt(1 - 8*a0)))
pB = (rB, v0 - rB) # point B
else:
rB = 2.5
rC = 6.5
nl = 5
dr = (rC - rB)/(nl - 1)
print("r0 : ", end='')
for i in range(nl):
rr = rB + i*dr
print(rr, end=', ')
tt = v0 - rr
sol = solve_ode(tt, -1, r0=rr)
outgeods.append(geod_line(sol))
sol = solve_ode(tt, tmax, r0=rr)
outgeods.append(geod_line(sol))
r0 : 3.00000000000000, 3.75000000000000, 4.50000000000000, 5.25000000000000, 6.00000000000000,
Drawing of the radiation region (yellow rays):
graph = draw_radiation_region(v0, rmax, 40)
Adding the outgoing radial null geodesics:
for geod in outgeods:
graph += geod
graph += cauchy_hor
show(graph, aspect_ratio=1, xmax=rmax, ymin=-4, ymax=ymax,
axes_labels=[r'$r/m$', r'$t/m$'], figsize=10)
show(graph, aspect_ratio=1, ymin=0, ymax=1.4, xmax=1,
axes_labels=[r'$r/m$', r'$t/m$'])
The ingoing null geodesics:
ingeods = ingoing_geods(-6., tmax, rmax, v0, 2., 2.)
for geod in ingeods:
graph += geod
graph += curvature_sing(tmax, 21)
graph += event_hor(tmax)
pA = (2, v0 - 2) # point A
Coordinate t at the event horizon birth [unit: m]: 0.001000000000369725 +/- 0.001
graph += trapping_hor(m_num(v), tmax)
graph_wo_points = copy(graph)
if S == S0:
graph += circle((0,0), 0.2, fill=True, color='black') \
+ text(r"$O$", (0.5, -0.6), fontsize=16, color='black')
graph += circle(pA, 0.2, fill=True, color='black') \
+ text(r"$A$", (pA[0]+0.7, pA[1]+0.6), fontsize=16, color='black')
graph += circle(pB, 0.2, fill=True, color='green') \
+ text(r"$B$", (pB[0]+0.8, pB[1]+0.2), fontsize=16, color='green')
graph += circle(pC, 0.2, fill=True, color='blue') \
+ text(r"$C$", (pC[0]+0.8, pC[1]+0.2), fontsize=16)
show(graph, aspect_ratio=1, xmax=rmax, ymin=ymin, ymax=ymax,
axes_labels=[r'$r/m$', r'$t/m$'], figsize=10)
filename = "vai_diag_naked_S0.pdf" if S == S0 else "vai_diag_naked_S2.pdf"
graph.save(filename, aspect_ratio=1, xmax=rmax, ymin=ymin,
ymax=ymax, axes_labels=[r'$r/m$', r'$t/m$'], figsize=10)
A zoom on the initial singularity:
graph = graph_wo_points
graph += circle((0,0), 0.05, fill=True, color='black') \
+ text(r"$O$", (0.12, -0.15), fontsize=20, color='black')
show(graph, aspect_ratio=1, ymin=-0.8, ymax=2, xmax=2,
#show(graph, aspect_ratio=0.05, ymin=0.5, ymax=3, xmax=0.2,
axes_labels=[r'$r/m$', r'$t/m$'])