#!/usr/bin/env python
# coding: utf-8

# # 3+1 slicing of Kerr spacetime
# 
# This worksheet demonstrates a few capabilities of <a href="http://sagemanifolds.obspm.fr/">SageManifolds</a> (version 1.0, as included in SageMath 7.5) in computations regarding the 3+1 slicing of Kerr spacetime.
# 
# Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Worksheets/v1.0/SM_Kerr_3p1.ipynb) 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()


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

# In[2]:


get_ipython().run_line_magic('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)


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

# In[5]:


X.<r,y,ph> = Sig.chart(r'r:(1,+oo) y:(-1,1) ph:(0,2*pi):\phi')
print(X) ; X


# <h3>Riemannian metric on $\Sigma$</h3>
# 
# <p>First the two Kerr parameters:</p>

# In[6]:


var('m, a', domain='real')
assume(m>0)
assume(a>0)
assumptions()


# <p>For dealing with extreme Kerr, the following must be uncommented:</p>

# In[7]:


# m = 1 ; a = 1


# <p>Some shortcut notations:</p>

# 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()


# <p>A matrix view of the components w.r.t. coordinates $(r,y,\phi)$:</p>

# In[10]:


gam[:]


# <h3>Lapse function and shift vector</h3>

# In[11]:


N = Sig.scalar_field(sqrt(Del / BB2), name='N')
print(N)
N.display()


# In[12]:


b = Sig.vector_field('beta', latex_name=r'\beta') 
b[3] = -2*m*r*a/AA2   
# unset components are zero 
b.display()


# ### 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()


# <p>Check (comparison with known formulas):</p>

# 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


# In[15]:


K[1,3] - Krp


# In[16]:


Kyp = 2*m*r*a^3*(1-y^2)*y*sqrt(Del)/rho2^2/sqrt(BB2)
Kyp


# In[17]:


K[2,3] - Kyp


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

# In[18]:


K1 = Sig.sym_bilin_form_field('K')
K1[1,3] = Krp
K1[2,3] = Kyp
K = K1
K.display()


# In[19]:


K.display_comp()


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

# In[20]:


Ku = K.up(gam, 0)
print(Ku) ; Ku.display()


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

# In[21]:


trK = Ku.trace()
print(trK)
trK.display()


# <h3>Connection and curvature</h3>
# <p>Let us call $D$ the Levi-Civita connection associated with $\gamma$: </p>

# In[22]:


D = gam.connection(name='D')
print(D) ; D


# <p>The Ricci tensor associated with $\gamma$:</p>

# In[23]:


Ric = gam.ricci()
print(Ric) ; Ric


# In[24]:


Ric.display_comp(only_nonredundant=True)


# In[25]:


Ric[1,1]


# In[26]:


Ric[1,2]


# In[27]:


Ric[1,3]


# In[28]:


Ric[2,2]


# In[29]:


Ric[2,3]


# In[30]:


Ric[3,3]


# <p>The scalar curvature $R = \gamma^{ij} R_{ij}$:</p>

# In[31]:


R = gam.ricci_scalar(name='R')
print(R)
R.display()


# <h2>3+1 Einstein equations</h2>
# <p>Let us check that the vacuum 3+1 Einstein equations are satisfied.</p>
# <p>We start by the contraint equations:</p>
# <h3>Hamiltonian constraint</h3>
# 
# <p>Let us first evaluate the term $K_{ij} K^{ij}$:</p>

# In[32]:


Kuu = Ku.up(gam, 1)
trKK = K['_ij']*Kuu['^ij']
print(trKK) ; trKK.display()


# 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()


# ### 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()


# <h3>Dynamical Einstein equations</h3>
# <p>Let us first evaluate the symmetric bilinear form $k_{ij} := K_{ik} K^k_{\ \, j}$:</p>

# In[35]:


KK = K['_ik']*Ku['^k_j']
print(KK)


# In[36]:


KK1 = KK.symmetrize()
KK == KK1


# In[37]:


KK = KK1
print(KK)


# In[38]:


KK.set_name('(KK)')
KK.display_comp()


# <p>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()


# ## 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)


# In[41]:


E.set_name('E')
E.display_comp(only_nonzero=False)


# 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()


# In[43]:


epsu = gam.volume_form(1)
print(epsu) ; epsu.display()


# In[44]:


DKu = D(Ku)
B = epsu['^k_il']*DKu['^l_jk']
print(B)


# <p>Let us check that $B$ is symmetric:</p>

# In[45]:


B1 = B.symmetrize()
B == B1


# <p>Accordingly, we set</p>

# In[46]:


B = B1
B.set_name('B')
print(B)


# In[47]:


B.display_comp(only_nonzero=False)