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

# # Simon-Mars tensor and Kerr spacetime
# 
# This worksheet demonstrates a few capabilities of [SageManifolds](http://sagemanifolds.obspm.fr) (version 1.0, as included in SageMath 7.5) regarding the computation of the Simon-Mars tensor.
# 
# Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Worksheets/v1.0/SM_Simon-Mars_Kerr.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)


# ## Spacetime manifold
# 
# We declare the Kerr spacetime (or more precisely the part of the Kerr spacetime covered by Boyer-Lindquist coordinates) as a 4-dimensional manifold $\mathcal{M}$:

# In[4]:


M = Manifold(4, 'M', latex_name=r'\mathcal{M}')
print(M)


# The standard **Boyer-Lindquist coordinates** $(t,r,\theta,\phi)$ are introduced by declaring a chart $X$ on $\mathcal{M}$, via the method `chart()`, the argument of which is a string expressing the coordinates names, their ranges (the default is $(-\infty,+\infty)$) and their LaTeX symbols:

# In[5]:


X.<t,r,th,ph> = M.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi') 
print(X) ; X


# <h2>Metric tensor</h2>
# 
# <p>The 2 parameters $m$ and $a$ of the Kerr spacetime are declared as symbolic variables:</p>

# In[6]:


var('m, a', domain='real')


# <p>Let us introduce the spacetime metric $g$ and set its components in the coordinate frame associated with Boyer-Lindquist coordinates, which is the current manifold's default frame:</p>

# In[7]:


g = M.lorentzian_metric('g')
rho2 = r^2 + (a*cos(th))^2
Delta = r^2 -2*m*r + a^2
g[0,0] = -(1-2*m*r/rho2)
g[0,3] = -2*a*m*r*sin(th)^2/rho2
g[1,1], g[2,2] = rho2/Delta, rho2
g[3,3] = (r^2+a^2+2*m*r*(a*sin(th))^2/rho2)*sin(th)^2
g.display()


# In[8]:


g[:]


# <p>The Levi-Civita connection $\nabla$ associated with $g$:</p>

# In[9]:


nabla = g.connection() ; print(nabla)


# <p>As a check, we verify that the covariant derivative of $g$ with respect to $\nabla$ vanishes identically:</p>

# In[10]:


nabla(g).display()


# <h2>Killing vector</h2>
# <p>The default vector frame on the spacetime manifold is the coordinate basis associated with Boyer-Lindquist coordinates:</p>

# In[11]:


M.default_frame() is X.frame()


# In[12]:


X.frame()


# <p>Let us consider the first vector field of this frame:</p>

# In[13]:


xi = X.frame()[0] ; xi


# In[14]:


print(xi)


# <p>The 1-form associated to it by metric duality is</p>

# In[15]:


xi_form = xi.down(g)
xi_form.set_name('xi_form', r'\underline{\xi}')
print(xi_form) ; xi_form.display()


# <p>Its covariant derivative is</p>

# In[16]:


nab_xi = nabla(xi_form)
print(nab_xi) ; nab_xi.display()


# <p>Let us check that the Killing equation is satisfied:</p>

# In[17]:


nab_xi.symmetrize().display()


# <p>Equivalently, we check that the Lie derivative of the metric along $\xi$ vanishes:</p>

# In[18]:


g.lie_der(xi).display()


# <p>Thank to Killing equation, $\nabla_g \underline{\xi}$ is antisymmetric. We may therefore define a 2-form by $F := - \nabla_g \xi$. Here we enforce the antisymmetry by calling the function <span style="font-family: courier new,courier;">antisymmetrize()</span> on <span style="font-family: courier new,courier;">nab_xi</span>:</p>

# In[19]:


F = - nab_xi.antisymmetrize()
F.set_name('F')
print(F)
F.display()


# <p>We check that</p>

# In[20]:


F == - nab_xi


# <p>The squared norm of the Killing vector is:</p>

# In[21]:


lamb = - g(xi,xi)
lamb.set_name('lambda', r'\lambda')
print(lamb)
lamb.display()


# <p>Instead of invoking $g(\xi,\xi)$, we could have evaluated $\lambda$ by means of the 1-form $\underline{\xi}$ acting on the vector field $\xi$:</p>

# In[22]:


lamb == - xi_form(xi)


# <p>or, using index notation as $\lambda = - \xi_a \xi^a$:</p>

# In[23]:


lamb == - ( xi_form['_a']*xi['^a'] )


# <h2>Curvature</h2>
# <p>The Riemann curvature tensor associated with $g$ is</p>

# In[24]:


Riem = g.riemann()
print(Riem)


# <p>The component $R^0_{\ \, 123} = R^t_{\ \, r\theta\phi}$ is</p>

# In[25]:


Riem[0,1,2,3]


# <p>The Ricci tensor:</p>

# In[26]:


Ric = g.ricci()
print(Ric)


# <p>Let us check that the Kerr metric is a vacuum solution of Einstein equation, i.e. that the Ricci tensor vanishes identically:</p>

# In[27]:


Ric.display()


# <p><span id="cell_outer_24">The Weyl conformal curvature tensor is<br /></span></p>

# In[28]:


C = g.weyl()
print(C)


# <p>Let us exhibit two of its components $C^0_{\ \, 123}$ and $C^0_{\ \, 101}$:</p>

# In[29]:


C[0,1,2,3]


# In[30]:


C[0,1,0,1]


# <p>To form the Simon-Mars tensor, we need the fully covariant (type-(0,4) tensor) form of the Weyl tensor (i.e. $C_{\alpha\beta\mu\nu} = g_{\alpha\sigma} C^\sigma_{\ \, \beta\mu\nu}$); we get it by lowering the first index with the metric:</p>

# In[31]:


Cd = C.down(g)
print(Cd)


# <p>The (monoterm) symmetries of this tensor are those inherited from the Weyl tensor, i.e. the antisymmetry on the last two indices (position 2 and 3, the first index being at position 0):</p>

# In[32]:


Cd.symmetries()


# <p>Actually, Cd is also antisymmetric with respect to the first two indices (positions 0 and 1), as we can check:</p>

# In[33]:


Cd == Cd.antisymmetrize(0,1)


# <p>To take this symmetry into account explicitely, we set</p>

# In[34]:


Cd = Cd.antisymmetrize(0,1)


# <p>Hence we have now</p>

# In[35]:


Cd.symmetries()


# <h2>Simon-Mars tensor</h2>
# <p>The Simon-Mars tensor with respect to the Killing vector $\xi$ is a rank-3 tensor introduced by Marc Mars in 1999 (<a href="http://iopscience.iop.org/0264-9381/16/7/323/">Class. Quantum Grav. <strong>16</strong>, 2507</a>). It has the remarkable property to vanish identically if, and only if, the spacetime $(\mathcal{M},g)$ is locally isometric to a Kerr spacetime.</p>
# <p>Let us evaluate the Simon-Mars tensor by following the formulas given in Mars' article. The starting point is the self-dual complex 2-form associated with the Killing 2-form $F$, i.e. the object $\mathcal{F} := F + i \, {}^* F$, where ${}^*F$ is the Hodge dual of $F$:</p>

# In[36]:


FF = F + I * F.hodge_dual(g)
FF.set_name('FF', r'\mathcal{F}') ; print(FF)


# In[37]:


FF.display()


# <p>Let us check that $\mathcal{F}$ is self-dual, i.e. that it obeys ${}^* \mathcal{F} = -i \mathcal{F}$:</p>

# In[38]:


FF.hodge_dual(g) == - I * FF


# Let us form the right self-dual of the Weyl tensor as follows
# $$\mathcal{C}_{\alpha\beta\mu\nu} = C_{\alpha\beta\mu\nu} + \frac{i}{2} \epsilon^{\rho\sigma}_{\ \ \ \mu\nu} \, C_{\alpha\beta\rho\sigma},$$
# where $\epsilon^{\rho\sigma}_{\ \ \ \mu\nu}$ is associated to the Levi-Civita tensor $\epsilon_{\rho\sigma\mu\nu}$ and is obtained by

# In[39]:


eps = g.volume_form(2)  # 2 = the first 2 indices are contravariant
print(eps)
eps.symmetries()


# <p>The right self-dual Weyl tensor is then</p>

# In[40]:


CC = Cd + I/2*( eps['^rs_..']*Cd['_..rs'] )
CC.set_name('CC', r'\mathcal{C}') ; print(CC)


# In[41]:


CC.symmetries()


# In[42]:


CC[0,1,2,3]


# <p>The Ernst 1-form $\sigma_\alpha = 2 \mathcal{F}_{\mu\alpha} \, \xi^\mu$ (0 = contraction on the first index of $\mathcal{F}$):</p>

# In[43]:


sigma = 2*FF.contract(0, xi)


# <p>Instead of invoking the function <span style="font-family: courier new,courier;">contract()</span>, we could have used the index notation to denote the contraction:</p>

# In[44]:


sigma == 2*( FF['_ma']*xi['^m'] )


# In[45]:


sigma.set_name('sigma', r'\sigma') ; print(sigma)
sigma.display()


# <p>The symmetric bilinear form $\gamma = \lambda \, g + \underline{\xi}\otimes\underline{\xi}$:</p>

# In[46]:


gamma = lamb*g + xi_form * xi_form
gamma.set_name('gamma', r'\gamma') ; print(gamma)
gamma.display()


# ### Final computation leading to the Simon-Mars tensor:
# 
# First we evaluate
# $$ S^{(1)}_{\alpha\beta\gamma} = 4 \mathcal{C}_{\mu\alpha\nu\beta} \, \xi^\mu \, \xi^\nu \, \sigma_\gamma $$

# In[47]:


S1 = 4*( CC.contract(0,xi).contract(1,xi) ) * sigma
print(S1)


# Then we form the tensor
# $$ S^{(2)}_{\alpha\beta\gamma} = - \gamma_{\alpha\beta} \, \mathcal{C}_{\rho\gamma\mu\nu} \, \xi^\rho \, \mathcal{F}^{\mu\nu} $$
# by first computing $\mathcal{C}_{\rho\gamma\mu\nu} \, \xi^\rho$:

# In[48]:


xiCC = CC['_.r..']*xi['^r']
print(xiCC)


# <p>We use the index notation to perform the double contraction $\mathcal{C}_{\gamma\rho\mu\nu} \mathcal{F}^{\mu\nu}$:</p>

# In[49]:


FFuu = FF.up(g)


# In[50]:


S2 = gamma * ( xiCC['_.mn']*FFuu['^mn'] )
print(S2)
S2.symmetries()


# The Simon-Mars tensor with respect to $\xi$ is obtained by antisymmetrizing $S^{(1)}$ and $S^{(2)}$ on their last two indices and adding them:
# $$ S_{\alpha\beta\gamma} = S^{(1)}_{\alpha[\beta\gamma]} + S^{(2)}_{\alpha[\beta\gamma]} $$
# We use the index notation for the antisymmetrization:

# In[51]:


S1A = S1['_a[bc]']
S2A = S2['_a[bc]']


# <p>An equivalent writing would have been (the last two indices being in position 1 and 2):</p>

# In[52]:


# S1A = S1.antisymmetrize(1,2)
# S2A = S2.antisymmetrize(1,2)


# <p>The Simon-Mars tensor is</p>

# In[53]:


S = S1A + S2A
S.set_name('S') ; print(S)
S.symmetries()


# In[54]:


S.display()


# **We thus recover the fact that the Simon-Mars tensor vanishes identically in Kerr spacetime.**
# 
# To check that the above computation was not trival, here is the component 112=$rr\theta$ for each of the two parts of the Simon-Mars tensor:

# In[55]:


S1A[1,1,2]


# In[56]:


S2A[1,1,2]


# In[57]:


S1A[1,1,2] + S2A[1,1,2]