This Jupyter/SageMath notebook is relative to the lectures Geometry and physics of black holes.
The involved computations are based on tools developed through the SageManifolds project.
NB: a version of SageMath at least equal to 8.8 is required to run this notebook:
version()
'SageMath version 9.1.beta8, Release Date: 2020-03-18'
First we set up the notebook to display mathematical objects using LaTeX formatting:
%display latex
To speed up computations, we ask for running them in parallel on 8 threads:
Parallelism().set(nproc=8)
We declare the spacetime manifold $M$:
M = Manifold(4, 'M', structure='Lorentzian')
print(M)
4-dimensional Lorentzian manifold M
We restrict the 3+1 Kerr patch to $r>0$, in order to introduce latter the Kerr-Schild coordinates:
X.<t,r,th,ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')
X
X.coord_range()
The Kerr parameters $m$ and $a$:
m = var('m', domain='real')
assume(m>0)
a = var('a', domain='real')
assume(a>=0)
We define the metric $g$ by its components w.r.t. the 3+1 Kerr coordinates:
g = M.metric()
rho2 = r^2 + (a*cos(th))^2
g[0,0] = -(1 - 2*m*r/rho2)
g[0,1] = 2*m*r/rho2
g[0,3] = -2*a*m*r*sin(th)^2/rho2
g[1,1] = 1 + 2*m*r/rho2
g[1,3] = -a*(1 + 2*m*r/rho2)*sin(th)^2
g[2,2] = rho2
g[3,3] = (r^2+a^2+2*m*r*(a*sin(th))^2/rho2)*sin(th)^2
g.display()
g.display_comp()
The inverse metric is pretty simple:
g.inverse()[:]
as well as the determinant w.r.t. to the 3+1 Kerr coordinates:
g.determinant().display()
g.determinant() == - (rho2*sin(th))^2
k = M.vector_field(1, -1, 0, 0, name='k')
k.display()
Let us check that $k$ is a null vector:
g(k,k).display()
Computation of $\nabla_k k$:
nabla = g.connection()
acc = nabla(k).contract(k)
acc.display()
k_form = k.down(g)
k_form.display()
Let us introduce the metric $f$ such that $$ g = f + 2 H \underline{k} \otimes \underline{k}$$ where $H = m r / \rho^2$:
H = M.scalar_field({X: m*r/rho2}, name='H')
H.display()
f = M.lorentzian_metric('f')
f.set(g - 2*H*(k_form*k_form))
f.display()
f[:]
$f$ is a flat metric:
f.riemann().display()
which proves that $g$ is a Kerr-Schild metric.
Let us check that $k$ is a null vector for $f$ as well:
f(k,k).expr()
KS.<t,x,y,z> = M.chart()
KS
X_to_KS = X.transition_map(KS, [t,
(r*cos(ph) - a*sin(ph))*sin(th),
(r*sin(ph) + a*cos(ph))*sin(th),
r*cos(th)])
X_to_KS.display()
R = sqrt((x^2 + y^2 + z^2 - a^2
+ sqrt((x^2 + y^2 + z^2 - a^2)^2 + 4*a^2*z^2))/2)
R
#X_to_KS.set_inverse(t, R, acos(z/R),
# atan2(R*y - a*x, R*x + a*y))
Check of the identity $$\frac{x^2 + y^2}{r^2 + a^2} + \frac{z^2}{r^2} = 1$$
((x^2 + y^2)/(R^2 + a^2) + z^2/R^2).simplify_full()
f.display(KS.frame())
Equivalently, we may check the following identity:
dt, dx, dy, dz = KS.coframe()[:]
f == - dt*dt + dx*dx + dy*dy + dz*dz
dx.display()
(dx*dx + dy*dy + dz*dz).display()
k.display(KS.frame())
k_form.display(KS.frame())
g.display(KS.frame())
Expression of the Killing vector $\partial/\partial\phi$ in terms of the Kerr-Schild frame:
X.frame()[3].display(KS.frame())
ap = 0.9 # value of a for the plots
rmax = 3
rcol = 'green' # color of the curves (th,ph) = const
thcol = 'red' # color of the curves (r,ph) = const
phcol = 'goldenrod' # color of the curves (r,th) = const
coordcol = {r: rcol, th: thcol, ph: phcol}
opacity = 1
surf_shift = 0.03
Numerical values of the event and Cauchy horizons:
rHp = 1 + sqrt(1 - ap^2)
rCp = 1 - sqrt(1 - ap^2)
rHp, rCp
X_plot = X.plot(KS, fixed_coords={t: 0, ph: 0}, ambient_coords=(x,y,z),
parameters={a: ap}, number_values=11,
color=coordcol, thickness=2, max_range=rmax,
label_axes=False)
X_to_KS(t, r, th, ph)
xyz_n = [s.subs({a: ap, ph: 0}) for s in X_to_KS(t, r, th, ph)[1:]]
xyz_n
xyz_H = [s.subs({r: rHp}) for s in xyz_n]
xyz_H
xyz_C = [s.subs({r: rCp}) for s in xyz_n]
xyz_C
xyz_n[1] += surf_shift # small adjustment to ensure that the surface does not cover
# the coordinate grid lines
graph = parametric_plot3d(xyz_n, (r, 0, rmax), (th, 0, pi), color='ivory',
opacity=opacity)
graph += X_plot
graph += parametric_plot3d(xyz_H, (th, 0, pi), color='black', thickness=6)
graph += parametric_plot3d(xyz_C, (th, 0, pi), color='blue', thickness=6)
graph += line([(0.03,0,0), (0.03,ap,0)], color='red', thickness=6)
graph
xyz_n = [s.subs({a: ap, ph: pi}) for s in X_to_KS(t, r, th, ph)[1:]]
xyz_H = [s.subs({r: rHp}) for s in xyz_n]
xyz_C = [s.subs({r: rCp}) for s in xyz_n]
X_plot_pi = X.plot(KS, fixed_coords={t: 0, ph: pi}, ambient_coords=(x,y,z),
parameters={a: ap}, number_values=11,
color=coordcol, thickness=2, max_range=rmax,
label_axes=False)
xyz_n[1] += surf_shift # small adjustment to ensure that the surface does not cover
# the coordinate grid lines
graph_pi = parametric_plot3d(xyz_n, (r, 0, rmax), (th, 0, pi), color='ivory',
opacity=opacity)
graph_pi += X_plot_pi
graph_pi += parametric_plot3d(xyz_H, (th, 0, pi), color='black', thickness=6)
graph_pi += parametric_plot3d(xyz_C, (th, 0, pi), color='blue', thickness=6)
graph_pi += line([(0.03,0,0), (0.03,-ap,0)], color='red', thickness=6)
graph_pi += line([(0,0,-1.1*rmax), (0,0,1.1*rmax)], color='green', thickness=4)
ph_slice = graph + graph_pi
show(ph_slice)
X.plot(KS, fixed_coords={t: 0, r: 1}, ambient_coords=(x,y,z),
parameters={a: ap}, number_values=11,
color=coordcol, label_axes=False)
The BH event horizon:
H_plot = X.plot(KS, fixed_coords={t: 0, r: rHp}, ambient_coords=(x,y,z),
parameters={a: ap}, number_values=11, color='grey',
label_axes=False)
H_plot
X.plot(KS, fixed_coords={t: 0, r: 0}, ambient_coords=(x,y,z),
parameters={a: ap}, number_values=11, color=coordcol,
label_axes=False)
rzero = X.plot(KS, fixed_coords={t: 0, r: 0}, ambient_coords=(x,y),
parameters={a: ap}, number_values={th: 17, ph: 13}, color=coordcol)
rzero += circle((0,0), ap, color='brown', thickness=3)
show(rzero, xmin=-1, xmax=1, ymin=-1, ymax=1, axes=False, frame=True,
axes_labels=[r'$x/m$', r'$y/m$'])
rzero.save('ksm_rzero_disk.pdf', xmin=-1, xmax=1, ymin=-1, ymax=1,
axes=False, frame=True, axes_labels=[r'$x/m$', r'$y/m$'])
Xth_pi2 = X.plot(KS, fixed_coords={t: 0, th: pi/2}, ambient_coords=(x,y,z),
parameters={a: ap}, number_values=11, color=coordcol,
max_range=rmax, thickness=1.5, label_axes=False)
Xth_pi2
X.plot(KS, fixed_coords={t: 0, th: pi/3}, ambient_coords=(x,y,z),
parameters={a: ap}, number_values=11, color=coordcol,
max_range=rmax, thickness=1.5, label_axes=False)
graph = X.plot(KS, fixed_coords={t: 0, th: pi/6}, ambient_coords=(x,y,z),
parameters={a: ap}, number_values=11, color=coordcol,
max_range=rmax, thickness=1.5, label_axes=False) \
+ Xth_pi2 \
+ circle((0,0,0), ap, color='pink', fill=True) \
+ circle((0,0,0), ap, color='brown', thickness=6, linestyle='dashed')
show(graph)
KS2.<t,xp,yp,zp> = M.chart(r"t xp:x' yp:y' zp:z'")
KS2
We use replace $r$ by $-r$ in the transformation formulas, because in what follows, we consider that $r$ is positive:
X_to_KS2 = X.transition_map(KS2, [t,
(-r*cos(ph) - a*sin(ph))*sin(th),
(-r*sin(ph) + a*cos(ph))*sin(th),
-r*cos(th)])
X_to_KS2.display()
X_plot2 = X.plot(KS2, fixed_coords={t: 0, ph: 0}, ambient_coords=(xp,yp,zp),
parameters={a: ap}, number_values=11,
color=coordcol, thickness=2, max_range=rmax,
label_axes=False)
surf_shift = 0
xyz_n2 = [s.subs({a: ap, ph: 0}) for s in X_to_KS2(t, r, th, ph)[1:]]
xyz_n2[1] -= surf_shift
graph2 = parametric_plot3d(xyz_n2, (r, 0, rmax), (th, 0, pi), color='pink',
opacity=opacity)
graph2 += X_plot2
X_plot2_pi = X.plot(KS2, fixed_coords={t: 0, ph: pi}, ambient_coords=(xp,yp,zp),
parameters={a: ap}, number_values=11,
color=coordcol, thickness=2, max_range=rmax,
label_axes=False)
xyz_n2 = [s.subs({a: ap, ph: pi}) for s in X_to_KS2(t, r, th, ph)[1:]]
xyz_n2[1] += surf_shift
graph2_pi = parametric_plot3d(xyz_n2, (r, 0, rmax), (th, 0, pi), color='pink',
opacity=opacity)
graph2_pi += X_plot2_pi
ph_slice2 = graph2 + graph2_pi
ph_slice2
ph_slice + ph_slice2