Kerr-Newman spacetime

This notebook demonstrates a few capabilities of SageMath in computations regarding Kerr-Newman spacetime. The corresponding tools have been developed within the SageManifolds project (version 1.3, as included in SageMath 8.3).

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:

In [1]:
version()
Out[1]:
'SageMath version 8.3, Release Date: 2018-08-03'

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

In [2]:
%display latex

We also define a viewer for 3D plots (use 'threejs' or 'jmol' for interactive 3D graphics):

In [3]:
viewer3D = 'threejs' # must be 'threejs', jmol', 'tachyon' or None (default)

Since some computations are quite heavy, we ask for running them in parallel on 8 threads:

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

Spacetime manifold

We declare the Kerr-Newman spacetime (or more precisely the part of it covered by Boyer-Lindquist coordinates) as a 4-dimensional Lorentzian manifold $\mathcal{M}$:

In [5]:
M = Manifold(4, 'M', latex_name=r'\mathcal{M}', structure='Lorentzian')

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:

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

Metric tensor

The 3 parameters $m$, $a$ and $q$ of the Kerr-Newman spacetime are declared as symbolic variables:

In [7]:
var('m a q')
Out[7]:

We get the (yet undefined) spacetime metric:

In [8]:
g = M.metric()

The metric is defined by its components in the coordinate frame associated with Boyer-Lindquist coordinates, which is the current manifold's default frame:

In [9]:
rho2 = r^2 + (a*cos(th))^2
Delta = r^2 -2*m*r + a^2 + q^2
g[0,0] = -1 + (2*m*r-q^2)/rho2
g[0,3] = -a*sin(th)^2*(2*m*r-q^2)/rho2
g[1,1], g[2,2] = rho2/Delta, rho2
g[3,3] = (r^2 + a^2 + (2*m*r-q^2)*(a*sin(th))^2/rho2)*sin(th)^2
g.display()
Out[9]:

The list of the non-vanishing components:

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

The component $g^{tt}$ of the inverse metric:

In [11]:
g.inverse()[0,0]
Out[11]:

The lapse function:

In [12]:
N = 1/sqrt(-(g.inverse()[[0,0]])); N
Out[12]:
In [13]:
N.display()
Out[13]:

Electromagnetic field tensor

Let us first introduce the 1-form basis associated with Boyer-Lindquist coordinates:

In [14]:
dBL = BL.coframe(); dBL
Out[14]:

The electromagnetic field tensor $F$ is formed as [cf. e.g. Eq. (33.5) of Misner, Thorne & Wheeler (1973)]

In [15]:
F = q/rho2^2 * (r^2-a^2*cos(th)^2)* dBL[1].wedge( dBL[0] - a*sin(th)^2* dBL[3] ) + \
    2*q/rho2^2 * a*r*cos(th)*sin(th)* dBL[2].wedge( (r^2+a^2)* dBL[3] - a* dBL[0] ) 
F.set_name('F')
F.display()
Out[15]:

The list of non-vanishing components:

In [16]:
F.display_comp()
Out[16]:

The Hodge dual of $F$:

In [17]:
star_F = F.hodge_dual(g)
star_F.display()
Out[17]:

Maxwell equations

Let us check that $F$ obeys the two (source-free) Maxwell equations:

In [18]:
F.exterior_derivative().display()
Out[18]:
In [19]:
star_F.exterior_derivative().display()
Out[19]:

Levi-Civita Connection

The Levi-Civita connection $\nabla$ associated with $g$:

In [20]:
nab = g.connection()
print(nab)
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:

In [21]:
nab(g) == 0
Out[21]:

Another view of the above property:

In [22]:
nab(g).display()
Out[22]:

The nonzero Christoffel symbols (skipping those that can be deduced by symmetry of the last two indices):

In [23]:
g.christoffel_symbols_display()
Out[23]:

Killing vectors

The default vector frame on the spacetime manifold is the coordinate basis associated with Boyer-Lindquist coordinates:

In [24]:
M.default_frame() is BL.frame()
Out[24]:
In [25]:
BL.frame()
Out[25]:

Let us consider the first vector field of this frame:

In [26]:
xi = BL.frame()[0] ; xi
Out[26]:
In [27]:
print(xi)
Vector field d/dt on the 4-dimensional Lorentzian manifold M

The 1-form associated to it by metric duality is

In [28]:
xi_form = xi.down(g)
xi_form.display()
Out[28]:

Its covariant derivative is

In [29]:
nab_xi = nab(xi_form)
print(nab_xi)
nab_xi.display()
Tensor field of type (0,2) on the 4-dimensional Lorentzian manifold M
Out[29]:

Let us check that the vector field $\xi=\frac{\partial}{\partial t}$ obeys Killing equation:

In [30]:
nab_xi.symmetrize() == 0
Out[30]:

Similarly, let us check that $\chi := \frac{\partial}{\partial\phi}$ is a Killing vector:

In [31]:
chi = BL.frame()[3] ; chi
Out[31]:
In [32]:
nab(chi.down(g)).symmetrize() == 0
Out[32]:

Another way to check that $\xi$ and $\chi$ are Killing vectors is the vanishing of the Lie derivative of the metric tensor along them:

In [33]:
g.lie_derivative(xi) == 0
Out[33]:
In [34]:
g.lie_derivative(chi) == 0
Out[34]:

Curvature

The Ricci tensor associated with $g$:

In [35]:
Ric = g.ricci()
print(Ric)
Field of symmetric bilinear forms Ric(g) on the 4-dimensional Lorentzian manifold M
In [36]:
Ric.display()
Out[36]:
In [37]:
Ric[:]
Out[37]:

Let us check that in the Kerr case, i.e. when $q=0$, the Ricci tensor is zero:

In [38]:
Ric[:].subs(q=0)
Out[38]:

The Riemann curvature tensor associated with $g$:

In [39]:
R = g.riemann()
print(R)
Tensor field Riem(g) of type (1,3) on the 4-dimensional Lorentzian manifold M

The component $R^0_{\ \, 101}$ of the Riemann tensor is

In [40]:
R[0,1,0,1]
Out[40]:

The expression in the uncharged limit (Kerr spacetime) is

In [41]:
R[0,1,0,1].expr().subs(q=0).simplify_rational()
Out[41]:

while in the non-rotating limit (Reissner-Nordström spacetime), it is

In [42]:
R[0,1,0,1].expr().subs(a=0).simplify_rational()
Out[42]:

In the Schwarzschild limit, it reduces to

In [43]:
R[0,1,0,1].expr().subs(a=0, q=0).simplify_rational()
Out[43]:

Obviously, it vanishes in the flat space limit:

In [44]:
R[0,1,0,1].expr().subs(m=0, a=0, q=0)
Out[44]:

Bianchi identity

Let us check the Bianchi identity $\nabla_p R^i_{\ \, j kl} + \nabla_k R^i_{\ \, jlp} + \nabla_l R^i_{\ \, jpk} = 0$:

In [45]:
DR = nab(R)   #long (takes a while)
print(DR)  
Tensor field nabla_g(Riem(g)) of type (1,4) on the 4-dimensional Lorentzian manifold M
In [46]:
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] ,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 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:

In [47]:
DR[0,1,2,3,1] + DR[0,1,3,1,2] + DR[0,1,1,2,3] # should be zero (Bianchi identity)
Out[47]:
In [48]:
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 -
Out[48]:

Ricci scalar

The Ricci scalar $R = g^{ab} R_{ab}$ of the Kerr-Newman spacetime vanishes identically:

In [49]:
g.ricci_scalar().display()
Out[49]:

Einstein equation

The Einstein tensor is

In [50]:
G = Ric - 1/2*g.ricci_scalar()*g
print(G)
Field of symmetric bilinear forms Ric(g)-unnamed metric on the 4-dimensional Lorentzian manifold M

Since the Ricci scalar is zero, the Einstein tensor reduces to the Ricci tensor:

In [51]:
G == Ric
Out[51]:

The invariant $F_{ab} F^{ab}$ of the electromagnetic field:

In [52]:
Fuu = F.up(g)
F2 = F['_ab']*Fuu['^ab'] ; print(F2)
Scalar field on the 4-dimensional Lorentzian manifold M
In [53]:
F2.display()
Out[53]:

The energy-momentum tensor of the electromagnetic field:

In [54]:
Fud = F.up(g,0)
T = 1/(4*pi)*( F['_k.']*Fud['^k_.'] - 1/4*F2 * g ); print(T)
Tensor field of type (0,2) on the 4-dimensional Lorentzian manifold M
In [55]:
T[:]
Out[55]:

Check of the Einstein equation:

In [56]:
G == 8*pi*T
Out[56]:

Kretschmann scalar

The tensor $R^\flat$, of components $R_{abcd} = g_{am} R^m_{\ \, bcd}$:

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

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

In [59]:
Kr_scalar = uR['^ijkl']*dR['_ijkl']
Kr_scalar.display()
Out[59]:

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:

In [60]:
Kr = Kr_scalar.coord_function()
Kr.factor()
Out[60]:

As a check, we can compare Kr to the formula given by R. Conn Henry, Astrophys. J. 535, 350 (2000):

In [61]:
Kr == 8/(r^2+(a*cos(th))^2)^6 *( 
          6*m^2*(r^6 - 15*r^4*(a*cos(th))^2 + 15*r^2*(a*cos(th))^4 - (a*cos(th))^6) 
        - 12*m*q^2*r*(r^4 - 10*(a*r*cos(th))^2 + 5*(a*cos(th))^4) 
        + q^4*(7*r^4 - 34*(a*r*cos(th))^2 + 7*(a*cos(th))^4) )
Out[61]:

The Schwarzschild value of the Kretschmann scalar is recovered by setting $a=0$ and $q=0$:

In [62]:
Kr.expr().subs(a=0, q=0)
Out[62]:

Let us plot the Kretschmann scalar for $m=1$, $a=0.9$ and $q=0.5$:

In [63]:
K1 = Kr.expr().subs(m=1, a=0.9, q=0.5)
plot3d(K1, (r,1,3), (th, 0, pi), viewer=viewer3D, online=True,
       axes_labels=['r', 'theta', 'Kr'])
Out[63]: