This notebook demonstrates a few capabilities of SageMath in computations regarding Kerr spacetime. The corresponding tools have been developed within the SageManifolds project.
Click here to download the notebook file (ipynb format). To run it, you must start SageMath within the Jupyter notebook, via the command sage -n jupyter
NB: a version of SageMath at least equal to 8.2 is required to run this notebook:
version()
First we set up the notebook to display mathematical objects using LaTeX rendering:
%display latex
and we initialize a time counter for benchmarking:
import time
comput_time0 = time.perf_counter()
Since some computations are quite heavy, we ask for running them in parallel on 8 threads:
Parallelism().set(nproc=8)
We declare the Kerr spacetime (or more precisely the part of it covered by Boyer-Lindquist coordinates) as a 4-dimensional Lorentzian manifold $\mathcal{M}$:
M = Manifold(4, 'M', latex_name=r'\mathcal{M}', structure='Lorentzian')
print(M)
We then introduce the standard Boyer-Lindquist coordinates as a chart BL
(for Boyer-Lindquist) on $\mathcal{M}$, via the method chart()
, the argument of which is a string
(delimited by r"..."
because of the backslash symbols) expressing the coordinates names, their ranges (the default is $(-\infty,+\infty)$) and their LaTeX symbols:
BL.<t,r,th,ph> = M.chart(r"t r th:(0,pi):\theta ph:(0,2*pi):\phi")
print(BL); BL
BL[0], BL[1]
The 2 parameters $m$ and $a$ of the Kerr spacetime are declared as symbolic variables:
var('m, a', domain='real')
We get the (yet undefined) spacetime metric:
g = M.metric()
The metric is set by its components in the coordinate frame associated with Boyer-Lindquist coordinates, which is the current manifold's default frame:
rho2 = r^2 + (a*cos(th))^2
Delta = r^2 -2*m*r + a^2
g[0,0] = -(1-2*m*r/rho2)
g[0,3] = -2*a*m*r*sin(th)^2/rho2
g[1,1], g[2,2] = rho2/Delta, rho2
g[3,3] = (r^2+a^2+2*m*r*(a*sin(th))^2/rho2)*sin(th)^2
g.display()
A matrix view of the components with respect to the manifold's default vector frame:
g[:]
The list of the non-vanishing components:
g.display_comp()
The Levi-Civita connection $\nabla$ associated with $g$:
nabla = g.connection()
print(nabla)
Let us verify that the covariant derivative of $g$ with respect to $\nabla$ vanishes identically:
nabla(g) == 0
Another view of the above property:
nabla(g).display()
The nonzero Christoffel symbols (skipping those that can be deduced by symmetry of the last two indices):
g.christoffel_symbols_display()
The default vector frame on the spacetime manifold is the coordinate basis associated with Boyer-Lindquist coordinates:
M.default_frame() is BL.frame()
BL.frame()
Let us consider the first vector field of this frame:
xi = BL.frame()[0]; xi
print(xi)
The 1-form associated to it by metric duality is
xi_form = xi.down(g)
xi_form.display()
Its covariant derivative is
nab_xi = nabla(xi_form)
print(nab_xi)
nab_xi.display()
Let us check that the Killing equation is satisfied:
nab_xi.symmetrize() == 0
Similarly, let us check that $\frac{\partial}{\partial\phi}$ is a Killing vector:
chi = BL.frame()[3] ; chi
nabla(chi.down(g)).symmetrize() == 0
The Ricci tensor associated with $g$:
Ric = g.ricci()
print(Ric)
Let us check that the Kerr metric is a solution of the vacuum Einstein equation:
Ric == 0
Another view of the above property:
Ric.display()
The Riemann curvature tensor associated with $g$:
R = g.riemann()
print(R)
Contrary to the Ricci tensor, the Riemann tensor does not vanish; for instance, the component $R^0_{\ \, 123}$ is
R[0,1,2,3]
Let us check the Bianchi identity $\nabla_p R^i_{\ \, j kl} + \nabla_k R^i_{\ \, jlp} + \nabla_l R^i_{\ \, jpk} = 0$:
DR = nabla(R) # long (takes a while)
print(DR)
# from __future__ import print_function # uncomment for SageMath version < 9.0 (Python 2 based)
for i in M.irange():
for j in M.irange():
for k in M.irange():
for l in M.irange():
for p in M.irange():
print(DR[i,j,k,l,p] + DR[i,j,l,p,k] + DR[i,j,p,k,l], end=' ')
If the last sign in the Bianchi identity is changed to minus, the identity does no longer hold:
DR[0,1,2,3,1] + DR[0,1,3,1,2] + DR[0,1,1,2,3] # should be zero (Bianchi identity)
DR[0,1,2,3,1] + DR[0,1,3,1,2] - DR[0,1,1,2,3] # note the change of the second + to -
The tensor $R^\flat$, of components $R_{abcd} = g_{am} R^m_{\ \, bcd}$:
dR = R.down(g)
print(dR)
The tensor $R^\sharp$, of components $R^{abcd} = g^{bp} g^{cq} g^{dr} R^a_{\ \, pqr}$:
uR = R.up(g)
print(uR)
The Kretschmann scalar $K := R^{abcd} R_{abcd}$:
Kr_scalar = uR['^{abcd}']*dR['_{abcd}']
Kr_scalar.display()
A variant of this expression can be obtained by invoking the factor()
method on the coordinate function representing the scalar field in the manifold's default chart:
Kr = Kr_scalar.coord_function()
Kr.factor()
As a check, we can compare Kr to the formula given by R. Conn Henry, Astrophys. J. 535, 350 (2000):
Kr == 48*m^2*(r^6 - 15*r^4*(a*cos(th))^2 + 15*r^2*(a*cos(th))^4
- (a*cos(th))^6) / (r^2+(a*cos(th))^2)^6
The Schwarzschild value of the Kretschmann scalar is recovered by setting $a=0$:
Kr.expr().subs(a=0)
Let us plot the Kretschmann scalar for $m=1$ and $a=0.9$:
K1 = Kr.expr().subs(m=1, a=0.9)
plot3d(K1, (r, 1, 3), (th, 0, pi), axes_labels=['r', 'theta', 'Kr'])
print("Total elapsed time: {} s".format(time.perf_counter() - comput_time0))