Kerr_Schild.ipynb

Kernel: SageMath 9.1.beta8

Kerr-Schild coordinates on Kerr spacetime

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)

Spacetime

We declare the spacetime manifold $M$:

M = Manifold(4, 'M', structure='Lorentzian')
print(M)

4-dimensional Lorentzian manifold M

3+1 Kerr coordinates $(t,r,\theta,\phi)$

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)

Kerr metric

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

Ingoing principal null geodesics

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

Kerr-Schild form of the Kerr metric

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

Kerr-Schild coordinates $(t, x, y, z)$

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

Minkowskian expression of $f$ in terms of Kerr-Schild coordinates:

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

Expression of $k$ and $g$ in the Kerr-Schild frame:

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

Plots

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)

Plot of the $r<0$ domain

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