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

# # 3+1 Einstein equations in the $\delta=2$ Tomimatsu-Sato spacetime
# 
# This worksheet demonstrates a few capabilities of SageMath in computations regarding the 3+1 slicing of the $\delta=2$ Tomimatsu-Sato spacetime. The corresponding tools have been developed within the  [SageManifolds](http://sagemanifolds.obspm.fr) project (version 1.1, as included in SageMath 8.1).
# 
# Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Worksheets/v1.1/SM_Tomimatsu-Sato_3p1.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)


# <h2>Tomimatsu-Sato spacetime</h2>
# <p>The Tomimatsu-Sato solution is an exact stationary and axisymmetric  solution of the vaccum Einstein equation, which is asymptotically flat and has a non-zero angular momentum. It has been found in 1972 by A. Tomimatsu and H. Sato [<a href="http://dx.doi.org/10.1103/PhysRevLett.29.1344">Phys. Rev. Lett. <strong>29</strong>, 1344 (1972)</a>], as a solution of the Ernst equation. It is actually the member $\delta=2$ of a larger family of solutions parametrized by a positive integer $\delta$ and exhibited by Tomimatsu and Sato in 1973 <a href="http://dx.doi.org/10.1143/PTP.50.95">[Prog. Theor. Phys. <strong>50</strong>, 95 (1973)</a>], the member $\delta=1$ being nothing but the Kerr metric. We refer to [<a href="http://dx.doi.org/10.1143/PTP.127.1057">Manko, Prog. Theor. Phys.<strong> 127</strong>, 1057 (2012)</a>] for a discussion of the properties of this solution. </p>
# <h2>Spacelike hypersurface</h2>
# <p>We consider some hypersurface $\Sigma$ of a spacelike foliation $(\Sigma_t)_{t\in\mathbb{R}}$ of $\delta=2$ Tomimatsu-Sato spacetime; we declare $\Sigma_t$ as a 3-dimensional manifold:</p>

# In[4]:


Sig = Manifold(3, 'Sigma', r'\Sigma', start_index=1)


# <p>On $\Sigma$, we consider the prolate spheroidal <span id="cell_outer_1">coordinates</span> $(x,y,\phi)$, with $x\in(1,+\infty)$, $y\in(-1,1)$ and $\phi\in(0,2\pi)$ :</p>

# In[5]:


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


# <h3>Riemannian metric on $\Sigma$</h3>
# <p>The Tomimatsu-Sato metric depens on three parameters: the integer $\delta$, the real number $p\in[0,1]$, and the total mass $m$:</p>

# In[6]:


var('d, p, m')
assume(m>0)
assumptions()


# <p>We set $\delta=2$ and choose a specific value for $p=1/5$:</p>

# In[7]:


d = 2
p = 1/5


# <p>Furthermore, without any loss of generality, we may set $m=1$ (this simply fixes some length scale):</p>

# In[8]:


m = 1


# <p>The parameter $q$ is related to $p$ by $p^2+q^2=1$:</p>

# In[9]:


q = sqrt(1-p^2)


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

# In[10]:


AA2 = (p^2*(x^2-1)^2+q^2*(1-y^2)^2)^2 \
      - 4*p^2*q^2*(x^2-1)*(1-y^2)*(x^2-y^2)^2
BB2 = (p^2*x^4+2*p*x^3-2*p*x+q^2*y^4-1)^2 \
      + 4*q^2*y^2*(p*x^3-p*x*y^2-y^2+1)^2
CC2 = p^3*x*(1-x^2)*(2*(x^4-1)+(x^2+3)*(1-y^2)) \
      + p^2*(x^2-1)*((x^2-1)*(1-y^2)-4*x^2*(x^2-y^2)) \
      + q^2*(1-y^2)^3*(p*x+1)


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

# In[11]:


gam = Sig.riemannian_metric('gam', latex_name=r'\gamma') 
gam[1,1] = m^2*BB2/(p^2*d^2*(x^2-1)*(x^2-y^2)^3)
gam[2,2] = m^2*BB2/(p^2*d^2*(y^2-1)*(-x^2+y^2)^3)
gam[3,3] = - m^2*(y^2-1)*(p^2*BB2^2*(x^2-1)
                          + 4*q^2*d^2*CC2^2*(y^2-1))/(AA2*BB2*d^2)
gam.display()


# A view of the non-vanishing components of $\gamma$ w.r.t. coordinates $(x,y,\phi)$:

# In[12]:


gam.display_comp()


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

# In[13]:


N2 = AA2/BB2 - 2*m*q*CC2*(y^2-1)/BB2*(2*m*q*CC2*(y^2-1)
                /(BB2*(m^2*(y^2-1)*(p^2*BB2^2*(x^2-1)
                +4*q^2*d^2*CC2^2*(y^2-1))/(AA2*BB2*d^2)))) 
N2.simplify_full()


# In[14]:


N = Sig.scalar_field(sqrt(N2.simplify_full()), name='N')
print(N)
N.display()


# The coordinate expression of the scalar field $N$:

# In[15]:


N.expr()


# In[16]:


b3 = 2*m*q*CC2*(y^2-1)/(BB2*(m^2*(y^2-1)*(p^2*BB2^2*(x^2-1)
        +4*q^2*d^2*CC2^2*(y^2-1))/(AA2*BB2*d^2)))
b = Sig.vector_field('beta', latex_name=r'\beta') 
b[3] = b3.simplify_full()
# unset components are zero 
b.display_comp(only_nonzero=False)


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


K = gam.lie_derivative(b) / (2*N)
K.set_name('K')
print(K)


# <p>The component $K_{13} = K_{x\phi}$:</p>

# In[18]:


K[1,3]


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

# In[19]:


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


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

# In[20]:


trK = Ku.trace()
trK


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

# In[21]:


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


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

# In[22]:


Ric = gam.ricci()
print(Ric)


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

# In[23]:


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


# The coordinate expression of the Ricci scalar is huge:

# In[24]:


R.expr()


# <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 constraint equations:</p>
# <h3>Hamiltonian constraint</h3>
# <p>Let us first evaluate the term $K_{ij} K^{ij}$:</p>

# In[25]:


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


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

# In[26]:


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


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

# In[27]:


mom = D(Ku).trace(0,2) - D(trK)
print(mom)
mom.display()


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

# In[28]:


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


# In[29]:


KK1 = KK.symmetrize()
KK == KK1


# In[30]:


KK = KK1
print(KK)


# <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[31]:


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


# <p>Hence the dynamical Einstein equations are satisfied.</p>
# 
# <p>Finally we have checked that all the 3+1 Einstein equations are satisfied by the $\delta=2$ Tomimatsu-Sato solution.</p>