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

# # 3+1 Simon-Mars tensor in Kerr spacetime
# 
# This worksheet demonstrates a few capabilities of [SageManifolds](http://sagemanifolds.obspm.fr) (version 1.0, as included in SageMath 7.5) in computations regarding 3+1 slicing of Kerr spacetime. In particular, it implements the computation of the 3+1 decomposition of the Simon-Mars tensor as given in the article [arXiv:1412.6542](http://arxiv.org/abs/1412.6542).
# 
# Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Worksheets/v1.0/SM_Simon-Mars_3p1_Kerr.ipynb) to download the worksheet file (ipynb format). To run it, you must start SageMath with 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)


# <p>The two Kerr parameters:</p>

# In[5]:


var('m, a')
assume(m>0)
assume(a>0)


# <h3>Riemannian metric on $\Sigma$</h3>
# <p>The variables introduced so far satisfy the following assumptions:</p>
# 
# <p>Without any loss of generality (for $m\not =0$), we may set $m=1$:</p>

# In[6]:


m=1
assume(a<1)


# In[7]:


#a=1 # extreme Kerr


# <p>On the hypersurface $\Sigma$, we are using coordinates $(r,y,\phi)$ that are related to the standard Boyer-Lindquist coordinates $(r,\theta,\phi)$ by $y=\cos\theta$:</p>

# In[8]:


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


# <h3>Riemannian metric on $\Sigma$</h3>
# <p>The variables introduced so far obey the following assumptions:</p>

# In[9]:


assumptions()


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

# In[10]:


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


# <p>The metric $h$ induced by the spacetime metric $g$ on $\Sigma$:</p>

# In[11]:


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[12]:


gam[:]


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

# In[13]:


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


# In[14]:


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


# <h3>Extrinsic curvature of $\Sigma$</h3>
# <p>We use the formula $$ K_{ij} = \frac{1}{2N} \mathcal{L}_{\beta} \gamma_{ij}$$ which is valid for any stationary spacetime:</p>

# In[15]:


K = gam.lie_der(b) / (2*N)
K.set_name('K')
print(K) ; K.display()


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

# In[16]:


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[17]:


K[1,3] - Krp


# In[18]:


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


# In[19]:


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[20]:


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


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

# In[21]:


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[22]:


trK = Ku.trace()
print(trK)


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

# In[23]:


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


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

# In[24]:


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


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


# <p>The vacuum Hamiltonian constraint equation is $$R + K^2 -K_{ij} K^{ij} = 0 $$</p>

# In[33]:


Ham = R + trK^2 - trKK
print(Ham) ; Ham.display()


# <h3>Momentum constraint</h3>
# <p>In vaccum, the momentum constraint is $$ D_j K^j_{\ \, i} - D_i K = 0 $$</p>

# 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[1,1]


# In[39]:


KK[1,2]


# In[40]:


KK[1,3]


# In[41]:


KK[2,2]


# In[42]:


KK[2,3]


# In[43]:


KK[3,3]


# <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 $$</p>

# In[44]:


dyn = K.lie_der(b) - D(D(N)) + N*( Ric + trK*K - 2*KK )
print(dyn)
dyn.display()


# <p>Hence, we have checked that all the vacuum 3+1 Einstein equations are fulfilled.</p>
# 
# <h2>Electric and magnetic parts of the Weyl tensor</h2>
# <p>The electric part is the bilinear form $E$ given by 
# $$E_{ij} = R_{ij} + K K_{ij} - K_{ik} K^k_{\ \, j} $$</p>

# In[45]:


E = Ric + trK*K - KK
print(E)


# In[46]:


E[1,1]


# In[47]:


E[1,1].factor()


# In[48]:


E[1,2]


# In[49]:


E[1,2].factor()


# In[50]:


E[1,3]


# In[51]:


E[2,2]


# In[52]:


E[2,2].factor()


# In[53]:


E[2,3]


# In[54]:


E[3,3]


# In[55]:


E[3,3].factor()


# <p>The magnetic part is the bilinear form $B$ defined by $$ B_{ij} = \epsilon^k_{\ \, l i} D_k K^l_{\ \, j}, $$</p>
# <p>where $\epsilon^k_{\ \, l i}$ 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_{\ \, l i} = \gamma^{km} \epsilon_{m l i}$. In SageManifolds, $\epsilon$ is obtained by the command <span style="font-family: courier new,courier;">volume_form()</span> and $\epsilon^\sharp$ by the command <span style="font-family: courier new,courier;">volume_form(1)</span> (1 = 1 index raised):</p>

# In[56]:


eps = gam.volume_form() 
print(eps) ; eps.display()


# In[57]:


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


# In[58]:


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


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

# In[59]:


B1 = B.symmetrize()
B == B1


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

# In[60]:


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


# In[61]:


B[1,1]


# In[62]:


B[1,1].factor()


# In[63]:


B[1,2]


# In[64]:


B[1,2].factor()


# In[65]:


B[1,3]


# In[66]:


B[2,2]


# In[67]:


B[2,2].factor()


# In[68]:


B[2,3]


# In[69]:


B[3,3]


# In[70]:


B[3,3].factor()


# <h2>3+1 decomposition of the Simon-Mars tensor</h2>
# <p>We follow the computation presented in <a href="http://arxiv.org/abs/1412.6542">arXiv:1412.6542</a>. We start by the tensor $E^\sharp$ of components $E^i_ {\ \, j}$:</p>

# In[71]:


Eu = E.up(gam, 0) 
print(Eu)


# <p>Tensor $B^\sharp$ of components $B^i_{\ \, j}$:</p>

# In[72]:


Bu = B.up(gam, 0)
print(Bu)


# <p>1-form $\beta^\flat$ of components $\beta_i$ and its exterior derivative:</p>

# In[73]:


bd = b.down(gam)
xdb = bd.exterior_derivative()
print(xdb) ; xdb.display()


# <p>Scalar square of shift $\beta_i \beta^i$:</p>

# In[74]:


b2 = bd(b)
print(b2) ; b2.display()


# <p>Scalar $Y = E(\beta,\beta) = E_{ij} \beta^i \beta^j$:</p>

# In[75]:


Ebb = E(b,b)
Y = Ebb
print(Y) ; Y.display()


# In[76]:


Ebb.coord_function().factor()


# In[77]:


Ebb.display()


# <p>Scalar $\bar Y = B(\beta,\beta) = B_{ij}\beta^i \beta^j$:</p>

# In[78]:


Bbb = B(b,b)
Y_bar = Bbb
print(Y_bar) ; Y_bar.display()


# In[79]:


Bbb.coord_function().factor()


# <p>1-form of components $Eb_i = E_{ij} \beta^j$:</p>

# In[80]:


Eb = E.contract(b)
print(Eb) ; Eb.display()


# <p>Vector field of components $Eub^i = E^i_{\ \, j} \beta^j$:</p>

# In[81]:


Eub = Eu.contract(b)
print(Eub) ; Eub.display()


# <p>1-form of components $Bb_i = B_{ij} \beta^j$:</p>

# In[82]:


Bb = B.contract(b)
print(Bb) ; Bb.display()


# <p>Vector field of components $Bub^i = B^i_{\ \, j} \beta^j$:</p>

# In[83]:


Bub = Bu.contract(b)
print(Bub) ; Bub.display()


# <p>Vector field of components $Kub^i = K^i_{\ \, j} \beta^j$:</p>

# In[84]:


Kub = Ku.contract(b)
print(Kub) ; Kub.display()


# In[85]:


T = 2*b(N) - 2*K(b,b)
print(T) ; T.display()


# In[86]:


Db = D(b)  #  Db^i_j = D_j b^i
Dbu = Db.up(gam, 1)  # Dbu^{ij} = D^j b^i
bDb = b*Dbu  # bDb^{ijk} = b^i D^k b^j
T_bar = eps['_ijk']*bDb['^ikj']
print(T_bar) ; T_bar.display()


# In[87]:


epsb = eps.contract(b) 
print(epsb)
epsb.display()


# In[88]:


epsB = eps['_ijl']*Bu['^l_k']
print(epsB)


# In[89]:


epsB.symmetries()


# In[90]:


epsB[1,2,3]


# In[91]:


Z = 2*N*( D(N)  -K.contract(b)) + b.contract(xdb)
print(Z) ; Z.display()


# In[92]:


DNu = D(N).up(gam)
A = 2*(DNu - Ku.contract(b))*b + N*Dbu
Z_bar = eps['_ijk']*A['^kj']
print(Z_bar) ; Z_bar.display()


# In[93]:


# Test:
Dbdu = D(bd).up(gam,1).up(gam,1)  # (Db)^{ij} = D^i b^j
A = 2*b*(DNu - Ku.contract(b)) + N*Dbdu
Z_bar0 = eps['_ijk']*A['^jk']  # NB: '^jk' and not 'kj'
Z_bar0 == Z_bar


# In[94]:


W = N*Eb + epsb.contract(Bub)
print(W) ; W.display()


# In[95]:


W_bar = N*Bb - epsb.contract(Eub)
print(W_bar) ; W_bar.display()


# In[96]:


W[3].factor()


# In[97]:


W_bar[3].factor()


# In[98]:


M = - 4*Eb(Kub - DNu) - 2*(epsB['_ij.']*Dbu['^ji'])(b)
print(M) ; M.display()


# In[99]:


M_bar = 2*(eps.contract(Eub))['_ij']*Dbu['^ji'] - 4*Bb(Kub - DNu)
print(M_bar) ; M_bar.display()


# In[100]:


A = epsB['_ilk']*b['^l'] + epsB['_ikl']*b['^l'] \
    + Bu['^m_i']*epsb['_mk'] - 2*N*E
xdbE = xdb['_kl']*Eu['^k_i']
L = 2*N*epsB['_kli']*Dbu['^kl'] + 2*xdb['_ij']*Eub['^j'] \
    + 2*xdbE['_li']*b['^l'] + 2*A['_ik']*(Kub - DNu)['^k']
print(L)


# In[101]:


L[1]


# In[102]:


L[1].factor()


# In[103]:


L[2]


# In[104]:


L[2].factor()


# In[105]:


L[3]


# In[106]:


N2pbb = N^2 + b2
V = N2pbb*E - 2*(b.contract(E)*bd).symmetrize() + Ebb*gam \
    + 2*N*(b.contract(epsB).symmetrize())
print(V)


# In[107]:


V[1,1]


# In[108]:


V[1,1].factor()


# In[109]:


V[1,2]


# In[110]:


V[1,2].factor()


# In[111]:


V[1,3]


# In[112]:


V[2,2]


# In[113]:


V[2,2].factor()


# In[114]:


V[2,3]


# In[115]:


V[3,3]


# In[116]:


V[3,3].factor()


# In[117]:


beps = b.contract(eps)
V_bar = N2pbb*B - 2*(b.contract(B)*bd).symmetrize() + Bbb*gam \
        - 2*N*(beps['_il']*Eu['^l_j']).symmetrize()
print(V_bar)


# In[118]:


V_bar[1,1]


# In[119]:


V_bar[1,1].factor()


# In[120]:


V_bar[1,2]


# In[121]:


V_bar[1,2].factor()


# In[122]:


V_bar[1,3]


# In[123]:


V_bar[2,2]


# In[124]:


V_bar[2,2].factor()


# In[125]:


V_bar[2,3]


# In[126]:


V_bar[3,3]


# In[127]:


V_bar[3,3].factor()


# In[128]:


G = (N^2 - b2)*gam + bd*bd
print(G)


# In[129]:


G.display()


# <h3>3+1 decomposition of the real part of the Simon-Mars tensor</h3>
# <p>We follow Eqs. (77)-(80) of <a href="http://arxiv.org/abs/1412.6542">arXiv:1412.6542</a>:</p>

# In[130]:


S1 = (4*(V*Z - V_bar*Z_bar) + G*L).antisymmetrize(1,2)
print(S1)


# In[131]:


S1.display()


# In[132]:


S2 = 4*(T*V - T_bar*V_bar - W*Z + W_bar*Z_bar) + M*G \
     - N*bd*L
print(S2)


# In[133]:


S2.display()


# In[134]:


S3 = (4*(W*Z - W_bar*Z_bar) + N*bd*L).antisymmetrize()
print(S3)


# In[135]:


S3.display()


# In[136]:


S2[3,1] == -2*S3[3,1]


# In[137]:


S2[3,2] == -2*S3[3,2]


# In[138]:


S4 = 4*(T*W - T_bar*W_bar) -4*(Y*Z - Y_bar*Z_bar) + N*M*bd \
     - b2*L
print(S4)


# In[139]:


S4.display()


# <p>Hence all the tensors $S^1$, $S^2$, $S^3$ and $S^4$ involved in the 3+1 decomposition of the real part of the Simon-Mars are zero, as they should since the Simon-Mars tensor vanishes identically for the Kerr spacetime.</p>
# 
# <h3>3+1 decomposition of the imaginary part of the Simon-Mars tensor</h3>
# 
# <p>We follow Eqs. (82)-(85) of <a href="http://arxiv.org/abs/1412.6542">arXiv:1412.6542</a>.</p>

# In[140]:


epsE = eps['_ijl']*Eu['^l_k']
print(epsE)


# In[141]:


A = - epsE['_ilk']*b['^l'] - epsE['_ikl']*b['^l'] \
     - Eu['^m_i']*epsb['_mk'] - 2*N*B
xdbB = xdb['_kl']*Bu['^k_i']
L_bar = - 2*N*epsE['_kli']*Dbu['^kl'] \
        + 2*xdb['_ij']*Bub['^j'] + 2*xdbB['_li']*b['^l'] \
        + 2*A['_ik']*(Kub - DNu)['^k']
print(L_bar)


# In[142]:


L_bar.display()


# In[143]:


S1_bar = (4*(V*Z_bar + V_bar*Z) + G*L_bar).antisymmetrize(1,2)
print(S1_bar)


# In[144]:


S1_bar.display()


# In[145]:


S2_bar = 4*(T_bar*V + T*V_bar) - 4*(W*Z_bar + W_bar*Z) \
         + M_bar*G - N*bd*L_bar
print(S2_bar)


# In[146]:


S2_bar.display()


# In[147]:


S3_bar = (4*(W*Z_bar + W_bar*Z) + N*bd*L_bar).antisymmetrize()
print(S3_bar)


# In[148]:


S3_bar.display()


# In[149]:


S4_bar = 4*(T_bar*W + T*W_bar - Y*Z_bar - Y_bar*Z) \
         + M_bar*N*bd - b2*L_bar
print(S4_bar)


# In[150]:


S4_bar.display()


# <p>Hence all the tensors ${\bar S}^1$, ${\bar S}^2$, ${\bar S}^3$ and ${\bar S}^4$ involved in the 3+1 decomposition of the imaginary part of the Simon-Mars are zero, as they should since the Simon-Mars tensor vanishes identically for the Kerr spacetime.</p>