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

# # Walker-Penrose Killing tensor in Kerr spacetime
# 
# This notebook demonstrates a few capabilities of SageMath in computations regarding Kerr spacetime. More precisely, it focuses of the Killing tensor $K$ found by Walker & Penrose [[Commun. Math. Phys. **18**, 265 (1970)](https://doi.org/10.1007/BF01649445)].
# The employed differential geometry tools have been developed within the  [SageManifolds](https://sagemanifolds.obspm.fr) project (version 1.3, as included in SageMath 8.3).
# 
# Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Worksheets/v1.3/SM_Kerr_Killing_tensor.ipynb) 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()


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

# In[2]:


get_ipython().run_line_magic('display', 'latex')


# To speed up the computations, 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 Boyer-Lindquist domain of Kerr spacetime) as a 4-dimensional Lorentzian manifold:

# In[4]:


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


# Let us declare the **Boyer-Lindquist coordinates** 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]:


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


# In[6]:


BL[0], BL[1]


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

# In[7]:


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


# We get the (yet undefined) spacetime metric by 

# In[8]:


g = M.metric()


# <p>The metric is set by its components in the coordinate frame associated with Boyer-Lindquist coordinates, which is the current manifold's default frame:</p>

# In[9]:


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


# <p>A matrix view of the components with respect to the manifold's default vector frame:</p>

# In[10]:


g[:]


# <p>The list of the non-vanishing components:</p>

# In[11]:


g.display_comp()


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

# In[12]:


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


# <p>Let us verify that the covariant derivative of $g$ with respect to $\nabla$ vanishes identically:</p>

# In[13]:


nabla(g).display()


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

# In[14]:


M.default_frame() is BL.frame()


# In[15]:


BL.frame()


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

# In[16]:


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


# In[17]:


print(xi)


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

# In[18]:


xi_form = xi.down(g) ; xi_form.display()


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

# In[19]:


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


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

# In[20]:


nab_xi.symmetrize() == 0


# <p>Similarly, let us check that $\frac{\partial}{\partial\phi}$ is a Killing vector:</p>

# In[21]:


chi = BL.frame()[3] ; chi


# In[22]:


nabla(chi.down(g)).symmetrize() == 0


# ## Principal null vectors
# 
# We introduce the principal null vectors $k$ and $\ell$ of Kerr spacetime:

# In[23]:


k = M.vector_field(name='k')
k[:] = [(r^2+a^2)/(2*rho2), -Delta/(2*rho2), 0, a/(2*rho2)]
k.display()


# In[24]:


el = M.vector_field(name='el', latex_name=r'\ell')
el[:] = [(r^2+a^2)/Delta, 1, 0, a/Delta]
el.display()


# Let us check that $k$ and $\ell$ are null vectors:

# In[25]:


g(k,k).expr()


# In[26]:


g(el,el).expr()


# Their scalar product is $-1$:

# In[27]:


g(k,el).expr()


# Note that the scalar product (with respect to metric $g$) can also be computed by means of the method `dot`:

# In[28]:


k.dot(el).expr()


# Let us evaluate the "acceleration" of $k$, i.e. $\nabla_k k$:

# In[29]:


acc_k = nabla(k).contract(k)
acc_k.display()


# We check that $k$ is a pregeodesic vector, i.e. that $\nabla_k k = \kappa_k k$ for some scalar field $\kappa_k$:

# In[30]:


for i in [0,1,3]:
    show(acc_k[i] / k[i])


# In[31]:


kappa_k = acc_k[[0]] / k[[0]]
kappa_k.display()


# In[32]:


acc_k == kappa_k * k


# Similarly let us evaluate the "acceleration" of $\ell$:

# In[33]:


acc_l = nabla(el).contract(el)
acc_l.display()


# Hence $\ell$ is a geodesic vector.

# ## Walker-Penrose Killing tensor
# 
# We need the 1-forms associated to $k$ and $\ell$ by metric duality:

# In[34]:


uk = k.down(g)
ul = el.down(g)


# The Walker-Penrose Killing tensor $K$ is then formed as
# $$ K = \rho^2 (\underline{\ell}\otimes \underline{k} + (\underline{k}\otimes \underline{\ell}) + r^2 g $$

# In[35]:


K = rho2*(ul*uk+ uk*ul) + r^2*g
K.set_name('K')
print(K)


# In[36]:


K.display_comp()


# In[37]:


DK = nabla(K)
print(DK)


# In[38]:


DK.display_comp()


# Let us check that $K$ is a Killing tensor:

# In[39]:


DK.symmetrize().display()


# Equivalently, we may write, using index notation:

# In[40]:


DK['_(abc)'].display()