Kerr-Newman spacetime

This worksheet demonstrates a few capabilities of SageManifolds (version 1.0, as included in SageMath 7.5) in computations regarding Kerr-Newman spacetime.

Click here to download the worksheet 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 7.5 is required to run this worksheet:

In [1]:
version()
Out[1]:
'SageMath version 7.5, Release Date: 2017-01-11'

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 long, we ask for running them in parallel on 8 cores:

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

Spacetime manifold

We declare the Kerr-Newman spacetime as a 4-dimensional diffentiable manifold:

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

Let us use the standard Boyer-Lindquist coordinates on it, by first introducing the part $\mathcal{M}_0$ covered by these coordinates

In [6]:
M0 = M.open_subset('M0', r'\mathcal{M}_0')
# BL = Boyer-Lindquist
BL.<t,r,th,ph> = M0.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi') 
print(BL) ; BL
Chart (M0, (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]:

Let us introduce the spacetime metric:

In [8]:
g = M.lorentzian_metric('g')

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 = M.diff_form(2, name='F')
F.set_restriction( 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.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 differentiable 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 Open subset M0 of the 4-dimensional differentiable 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 Open subset M0 of the 4-dimensional differentiable 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 differentiable 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 differentiable 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) ; print(DR)  #long (takes a while)
Tensor field nabla_g(Riem(g)) of type (1,4) on the 4-dimensional differentiable 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) on the 4-dimensional differentiable 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 differentiable 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 differentiable 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 differentiable 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 differentiable 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, axes_labels=['r', 'theta', 'Kr'])
Out[63]: