This notebook demonstrates a few capabilities of SageMath in computations regarding Kerr spacetime. The corresponding tools have been developed within the SageManifolds project.
NB: a version of SageMath at least equal to 8.2 is required to run this notebook:
version()
'SageMath version 9.5.beta2, Release Date: 2021-09-26'
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)
4-dimensional Lorentzian manifold 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
Chart (M, (t, r, th, ph))
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)
Levi-Civita connection nabla_g associated with the Lorentzian metric g on the 4-dimensional Lorentzian manifold M
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)
Vector field ∂/∂t on the 4-dimensional Lorentzian manifold M
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()
Tensor field of type (0,2) on the 4-dimensional Lorentzian manifold M
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)
Field of symmetric bilinear forms Ric(g) on the 4-dimensional Lorentzian manifold M
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)
Tensor field Riem(g) of type (1,3) on the 4-dimensional Lorentzian manifold M
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)
Tensor field nabla_g(Riem(g)) of type (1,4) on the 4-dimensional Lorentzian manifold M
# 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=' ')
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
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)
Tensor field of type (0,4) on the 4-dimensional Lorentzian manifold M
The tensor $R^\sharp$, of components $R^{abcd} = g^{bp} g^{cq} g^{dr} R^a_{\ \, pqr}$:
uR = R.up(g)
print(uR)
Tensor field of type (4,0) on the 4-dimensional Lorentzian manifold M
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))
Total elapsed time: 449.14473272299983 s