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:

In [1]:
version()
Out[1]:
'SageMath version 9.1.beta8, Release Date: 2020-03-18'

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

In [2]:
%display latex

To speed up computations, we ask for running them in parallel on 8 threads:

In [3]:
Parallelism().set(nproc=8)

Spacetime

We declare the spacetime manifold $M$:

In [4]:
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:

In [5]:
X.<t,r,th,ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')
X
Out[5]:
In [6]:
X.coord_range()
Out[6]:

The Kerr parameters $m$ and $a$:

In [7]:
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:

In [8]:
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()
Out[8]:
In [9]:
g.display_comp()
Out[9]:

The inverse metric is pretty simple:

In [10]:
g.inverse()[:]
Out[10]:

as well as the determinant w.r.t. to the 3+1 Kerr coordinates:

In [11]:
g.determinant().display()
Out[11]:
In [12]:
g.determinant() == - (rho2*sin(th))^2
Out[12]:

Ingoing principal null geodesics

In [13]:
k = M.vector_field(1, -1, 0, 0, name='k')
k.display()
Out[13]:

Let us check that $k$ is a null vector:

In [14]:
g(k,k).display()
Out[14]:

Computation of $\nabla_k k$:

In [15]:
nabla = g.connection()
acc = nabla(k).contract(k)
acc.display()
Out[15]:
In [16]:
k_form = k.down(g)
k_form.display()
Out[16]:

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$:

In [17]:
H = M.scalar_field({X: m*r/rho2}, name='H')
H.display()
Out[17]:
In [18]:
f = M.lorentzian_metric('f')
f.set(g - 2*H*(k_form*k_form))
f.display()
Out[18]:
In [19]:
f[:]
Out[19]:

$f$ is a flat metric:

In [20]:
f.riemann().display()
Out[20]:

which proves that $g$ is a Kerr-Schild metric.

Let us check that $k$ is a null vector for $f$ as well:

In [21]:
f(k,k).expr()
Out[21]:

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

In [22]:
KS.<t,x,y,z> = M.chart()
KS
Out[22]:
In [23]:
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()
Out[23]:
In [24]:
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
Out[24]:
In [25]:
#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$$

In [26]:
((x^2 + y^2)/(R^2 + a^2) + z^2/R^2).simplify_full()
Out[26]:

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

In [27]:
f.display(KS.frame())
Out[27]:

Equivalently, we may check the following identity:

In [28]:
dt, dx, dy, dz = KS.coframe()[:]
f == - dt*dt + dx*dx + dy*dy + dz*dz
Out[28]:
In [29]:
dx.display()
Out[29]:
In [30]:
(dx*dx + dy*dy + dz*dz).display()
Out[30]:

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

In [31]:
k.display(KS.frame())
Out[31]:
In [32]:
k_form.display(KS.frame())
Out[32]:
In [33]:
g.display(KS.frame())
Out[33]:

Expression of the Killing vector $\partial/\partial\phi$ in terms of the Kerr-Schild frame:

In [34]:
X.frame()[3].display(KS.frame())
Out[34]:

Plots

In [35]:
ap = 0.9  # value of a for the plots
rmax = 3
In [36]:
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}
In [37]:
opacity = 1
surf_shift = 0.03

Numerical values of the event and Cauchy horizons:

In [38]:
rHp = 1 + sqrt(1 - ap^2)
rCp = 1 - sqrt(1 - ap^2)
rHp, rCp
Out[38]:
In [39]:
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)
In [40]:
X_to_KS(t, r, th, ph)
Out[40]:
In [41]:
xyz_n = [s.subs({a: ap, ph: 0}) for s in X_to_KS(t, r, th, ph)[1:]]
xyz_n
Out[41]:
In [42]:
xyz_H = [s.subs({r: rHp}) for s in xyz_n]
xyz_H
Out[42]:
In [43]:
xyz_C = [s.subs({r: rCp}) for s in xyz_n]
xyz_C
Out[43]:
In [44]:
xyz_n[1] += surf_shift  # small adjustment to ensure that the surface does not cover
                        # the coordinate grid lines
In [45]:
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
Out[45]:
In [46]:
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)
In [47]:
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)