3+1 slicing of Kerr spacetime

This worksheet demonstrates a few capabilities of SageManifolds (version 1.0, as included in SageMath 7.5) in computations regarding the 3+1 slicing of Kerr 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.1, Release Date: 2017-01-15'

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 manifold:

In [4]:
Sig = Manifold(3, 'Sigma', r'\Sigma', start_index=1)
print(Sig)
3-dimensional differentiable 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)
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.riemannian_metric('gam', latex_name=r'\gamma') 
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 differentiable 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 differentiable 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 differentiable 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 differentiable 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 differentiable 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 differentiable manifold Sigma
Out[23]:
In [24]:
Ric.display_comp(only_nonredundant=True)
Out[24]:
In [25]:
Ric[1,1]
Out[25]:
In [26]:
Ric[1,2]
Out[26]:
In [27]:
Ric[1,3]
Out[27]:
In [28]:
Ric[2,2]
Out[28]:
In [29]:
Ric[2,3]
Out[29]:
In [30]:
Ric[3,3]
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 differentiable manifold Sigma
Out[31]:

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 [32]:
Kuu = Ku.up(gam, 1)
trKK = K['_ij']*Kuu['^ij']
print(trKK) ; trKK.display()
Scalar field on the 3-dimensional differentiable manifold Sigma
Out[32]:

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

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

Momentum constraint

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

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

Dynamical Einstein equations

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

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

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 [39]:
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 differentiable manifold Sigma
Out[39]:

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 [40]:
E = Ric + trK*K - KK
print(E)
Field of symmetric bilinear forms +Ric(gam)-(KK) on the 3-dimensional differentiable manifold Sigma
In [41]:
E.set_name('E')
E.display_comp(only_nonzero=False)
Out[41]:

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 [42]:
eps = gam.volume_form() 
print(eps) ; eps.display()
3-form eps_gam on the 3-dimensional differentiable manifold Sigma
Out[42]:
In [43]:
epsu = gam.volume_form(1)
print(epsu) ; epsu.display()
Tensor field of type (1,2) on the 3-dimensional differentiable manifold Sigma
Out[43]:
In [44]:
DKu = D(Ku)
B = epsu['^k_il']*DKu['^l_jk']
print(B)
Tensor field of type (0,2) on the 3-dimensional differentiable manifold Sigma

Let us check that $B$ is symmetric:

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

Accordingly, we set

In [46]:
B = B1
B.set_name('B')
print(B)
Field of symmetric bilinear forms B on the 3-dimensional differentiable manifold Sigma
In [47]:
B.display_comp(only_nonzero=False)
Out[47]: