3+1 slicing of Kerr spacetime

This notebook demonstrates a few capabilities of SageMath in computations regarding the 3+1 slicing of Kerr 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

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

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

Spacelike hypersurface

We consider some hypersurface $\Sigma$ of a spacelike foliation $(\Sigma_t)_{t\in\mathbb{R}}$ of Kerr spacetime; we declare $\Sigma_t$ as a 3-dimensional Riemannian manifold:

In [4]:
Sig = Manifold(3, 'Sigma', latex_name=r'\Sigma', 
               structure='Riemannian', metric_name='gam', 
               metric_latex_name=r'\gamma', start_index=1)
print(Sig)
3-dimensional Riemannian manifold Sigma

On $\Sigma$, we consider the "rational-polynomial" coordinates $(r,y,\phi)$ inheritated from the standard Boyer-Lindquist coordinates $(t,r,\theta,\phi)$ of Kerr spacetime, via $y=\cos\theta$:

In [5]:
X.<r,y,ph> = Sig.chart(r'r:(1,+oo) y:(-1,1) ph:(0,2*pi):\phi')
print(X) ; X
Chart (Sigma, (r, y, ph))
Out[5]:

Riemannian metric on $\Sigma$

First the two Kerr parameters:

In [6]:
var('m, a', domain='real')
assume(m > 0)
assume(a > 0)
assume(a^2-2*m*r+r^2 > 0)  # region where Sigma is spacelike
assumptions()
Out[6]:

For dealing with extreme Kerr, the following must be uncommented:

In [7]:
# m = 1 ; a = 1

Some shortcut notations:

In [8]:
rho2 = r^2 + a^2*y^2 
Del = r^2 -2*m*r + a^2
AA2 = rho2*(r^2 + a^2) + 2*a^2*m*r*(1-y^2)
BB2 = r^2 + a^2 + 2*a^2*m*r*(1-y^2)/rho2

The metric $\gamma$ induced by the spacetime metric $g$ on $\Sigma$:

In [9]:
gam = Sig.metric()
gam[1,1] = rho2/Del
gam[2,2] = rho2/(1-y^2)
gam[3,3] = BB2*(1-y^2)
gam.display()
Out[9]:

A matrix view of the components w.r.t. coordinates $(r,y,\phi)$:

In [10]:
gam[:]
Out[10]:

Lapse function and shift vector

In [11]:
N = Sig.scalar_field(sqrt(Del / BB2), name='N')
print(N)
N.display()
Scalar field N on the 3-dimensional Riemannian manifold Sigma
Out[11]:
In [12]:
b = Sig.vector_field('beta', latex_name=r'\beta') 
b[3] = -2*m*r*a/AA2   
# unset components are zero 
b.display()
Out[12]:

Extrinsic curvature of $\Sigma$

We use the formula $$ K_{ij} = \frac{1}{2N} \mathcal{L}_{\beta} \gamma_{ij} $$ which is valid for any stationary spacetime:

In [13]:
K = gam.lie_der(b) / (2*N)
K.set_name('K')
print(K) ; K.display()
Field of symmetric bilinear forms K on the 3-dimensional Riemannian manifold Sigma
Out[13]:

Check (comparison with known formulas):

In [14]:
Krp = a*m*(1-y^2)*(3*r^4+a^2*r^2+a^2*(r^2-a^2)*y^2) / rho2^2/sqrt(Del*BB2)
Krp
Out[14]:
In [15]:
K[1,3] - Krp
Out[15]:
In [16]:
Kyp = 2*m*r*a^3*(1-y^2)*y*sqrt(Del)/rho2^2/sqrt(BB2)
Kyp
Out[16]:
In [17]:
K[2,3] - Kyp
Out[17]:

For now on, we use the expressions Krp and Kyp above for $K_{r\phi}$ and $K_{ry}$, respectively:

In [18]:
K1 = Sig.sym_bilin_form_field('K')
K1[1,3] = Krp
K1[2,3] = Kyp
K = K1
K.display()
Out[18]:
In [19]:
K.display_comp()
Out[19]:

The type-(1,1) tensor $K^\sharp$ of components $K^i_{\ \, j} = \gamma^{ik} K_{kj}$:

In [20]:
Ku = K.up(gam, 0)
print(Ku) ; Ku.display()
Tensor field of type (1,1) on the 3-dimensional Riemannian manifold Sigma
Out[20]:

We may check that the hypersurface $\Sigma$ is maximal, i.e. that $K^k_{\ \, k} = 0$:

In [21]:
trK = Ku.trace()
print(trK)
trK.display()
Scalar field zero on the 3-dimensional Riemannian manifold Sigma
Out[21]:

Connection and curvature

Let us call $D$ the Levi-Civita connection associated with $\gamma$:

In [22]:
D = gam.connection(name='D')
print(D) ; D
Levi-Civita connection D associated with the Riemannian metric gam on the 3-dimensional Riemannian manifold Sigma
Out[22]:

The Ricci tensor associated with $\gamma$:

In [23]:
Ric = gam.ricci()
print(Ric) ; Ric
Field of symmetric bilinear forms Ric(gam) on the 3-dimensional Riemannian manifold Sigma
Out[23]:
In [24]:
Ric.display_comp(only_nonredundant=True)
Out[24]:
In [25]:
Ric[1,1].factor()
Out[25]:
In [26]:
Ric[1,2].factor()
Out[26]:
In [27]:
Ric[1,3]
Out[27]:
In [28]:
Ric[2,2].factor()
Out[28]:
In [29]:
Ric[2,3]
Out[29]:
In [30]:
Ric[3,3].factor()
Out[30]:

The scalar curvature $R = \gamma^{ij} R_{ij}$:

In [31]:
R = gam.ricci_scalar(name='R')
print(R)
R.display()
Scalar field R on the 3-dimensional Riemannian manifold Sigma
Out[31]:
In [32]:
R.expr().factor()
Out[32]:

3+1 Einstein equations

Let us check that the vacuum 3+1 Einstein equations are satisfied.

We start by the contraint equations:

Hamiltonian constraint

Let us first evaluate the term $K_{ij} K^{ij}$:

In [33]:
Kuu = Ku.up(gam, 1)
trKK = K['_ij']*Kuu['^ij']
print(trKK) ; trKK.display()
Scalar field on the 3-dimensional Riemannian manifold Sigma
Out[33]:

The vacuum Hamiltonian constraint equation is $$ R + K^2 -K_{ij} K^{ij} = 0 $$

In [34]:
Ham = R + trK^2 - trKK
print(Ham) ; Ham.display()
Scalar field on the 3-dimensional Riemannian manifold Sigma
Out[34]:

Momentum constraint

In vaccum, the momentum constraint is $$ D_j K^j_{\ \, i} - D_i K = 0 $$

In [35]:
mom = D(Ku).trace(0,2) - D(trK)
print(mom)
mom.display()
1-form on the 3-dimensional Riemannian manifold Sigma
Out[35]:

Dynamical Einstein equations

Let us first evaluate the symmetric bilinear form $k_{ij} := K_{ik} K^k_{\ \, j}$:

In [36]:
KK = K['_ik']*Ku['^k_j']
print(KK)
Tensor field of type (0,2) on the 3-dimensional Riemannian manifold Sigma
In [37]:
KK1 = KK.symmetrize()
KK == KK1
Out[37]:
In [38]:
KK = KK1
print(KK)
Field of symmetric bilinear forms on the 3-dimensional Riemannian manifold Sigma
In [39]:
KK.set_name('(KK)')
KK.display_comp()
Out[39]:

In vacuum and for stationary spacetimes, the dynamical Einstein equations are $$ \mathcal{L}_\beta K_{ij} - D_i D_j N + N \left( R_{ij} + K K_{ij} - 2 K_{ik} K^k_{\ \, j}\right) = 0 $$

In [40]:
dyn = K.lie_der(b) - D(D(N)) + N*( Ric + trK*K - 2*KK )
print(dyn)
dyn.display()
Tensor field of type (0,2) on the 3-dimensional Riemannian manifold Sigma
Out[40]:

Electric and magnetic parts of the Weyl tensor

The electric part is the bilinear form $E$ given by $$ E_{ij} = R_{ij} + K K_{ij} - K_{ik} K^k_{\ \, j} $$

In [41]:
E = Ric + trK*K - KK
print(E)
Field of symmetric bilinear forms on the 3-dimensional Riemannian manifold Sigma
In [42]:
E.set_name('E')
E.display_comp(only_nonzero=False)
Out[42]:

Let us factorize the components for a shorter display:

In [43]:
for i in range(1,4):
    for j in range(i,4):
        comp = E[i,j]
        if comp:  # factorization attempted for nonzero components
            comp.factor()
E.display_comp(only_nonzero=False)
Out[43]:

The magnetic part is the bilinear form $B$ defined by $$ B_{ij} = \epsilon^k_{\ \, i l} D_k K^l_{\ \, j}, $$ where $\epsilon^k_{\ \, i l}$ are the components of the type-(1,2) tensor $\epsilon^\sharp$, related to the Levi-Civita alternating tensor $\epsilon$ associated with $\gamma$ by $\epsilon^k_{\ \, i l} = \gamma^{km} \epsilon_{m i l}$. In SageManifolds, $\epsilon$ is obtained by the command volume_form() and $\epsilon^\sharp$ by the command volume_form(1) (1 = one index raised):

In [44]:
eps = gam.volume_form() 
print(eps) ; eps.display()
3-form eps_gam on the 3-dimensional Riemannian manifold Sigma
Out[44]:
In [45]:
epsu = gam.volume_form(1)
print(epsu) ; epsu.display()
Tensor field of type (1,2) on the 3-dimensional Riemannian manifold Sigma
Out[45]:
In [46]:
DKu = D(Ku)
B = epsu['^k_il']*DKu['^l_jk']
print(B)
Tensor field of type (0,2) on the 3-dimensional Riemannian manifold Sigma

Let us check that $B$ is symmetric:

In [47]:
B1 = B.symmetrize()
B == B1
Out[47]:

Accordingly, we set

In [48]:
B = B1
B.set_name('B')
print(B)
Field of symmetric bilinear forms B on the 3-dimensional Riemannian manifold Sigma

As we did for $E$, we factorize the components of $B$ before displaying them:

In [49]:
for i in range(1,4):
    for j in range(i,4):
        comp = B[i,j]
        if comp:  # factorization attempted for nonzero components
            comp.factor()
B.display_comp(only_nonzero=False)
Out[49]: