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

# # Kerr-Newman spacetime
# 
# This notebook demonstrates a few capabilities of SageMath in computations regarding the Kerr-Newman spacetime. The corresponding tools have been developed within the  [SageManifolds](https://sagemanifolds.obspm.fr) project.

# *NB:* a version of SageMath at least equal to 9.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')


# and we initialize a time counter for benchmarking:

# In[3]:


import time
comput_time0 = time.perf_counter()


# Since some computations are quite heavy, we ask for running them in parallel on 8 threads:

# In[4]:


Parallelism().set(nproc=8)


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

# In[5]:


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


# We then introduce the standard <strong>Boyer-Lindquist coordinates</strong> as a chart `BL` (for *Boyer-Lindquist*) on $\mathcal{M}$, via the method `chart()`, the argument of which is a string
# (delimited by `r"..."` because of the backslash symbols) expressing the coordinates names, their ranges (the default is $(-\infty,+\infty)$) and their LaTeX symbols:

# In[6]:


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


# ## Metric tensor
# 
# The 3 parameters $m$, $a$ and $q$ of the Kerr-Newman spacetime are declared as symbolic variables:

# In[7]:


var('m a q')


# We get the (yet undefined) spacetime metric:

# In[8]:


g = M.metric()


# <p>The metric is defined 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 + q^2
g[0,0] = -1 + (2*m*r-q^2)/rho2
g[0,3] = -a*sin(th)^2*(2*m*r-q^2)/rho2
g[1,1], g[2,2] = rho2/Delta, rho2
g[3,3] = (r^2 + a^2 + (2*m*r-q^2)*(a*sin(th))^2/rho2)*sin(th)^2
g.display()


# <p><span id="cell_outer_11">The list of the non-vanishing components:</span></p>

# In[10]:


g.display_comp()


# <p>The component $g^{tt}$ of the inverse metric:</p>

# In[11]:


g.inverse()[0,0]


# <p>The lapse function:</p>

# In[12]:


N = 1/sqrt(-(g.inverse()[[0,0]])); N


# In[13]:


N.display()


# ## Electromagnetic field tensor
# 
# Let us first get the 1-forms $(\mathrm{d}t, \mathrm{d}r, \mathrm{d}\theta, \mathrm{d}\phi)$ as those forming
# the coframe associated with Boyer-Lindquist coordinates:

# In[14]:


BL.coframe()


# In[15]:


dt, dr, dth, dph = BL.coframe()[:]
dt, dph


# The electromagnetic 4-potential 1-form $A$ of the Kerr-Newman solution is

# In[16]:


A = - q*r/rho2 * (dt - a*sin(th)^2*dph)
A.set_name('A')
A.display()


# The electromagnetic field tensor $F$ is computed as the exterior derivative of the 4-potential: $F = \mathrm{d} A$. In Sage, the exterior derivative of a $p$-form is returned by the function `diff()`:

# In[17]:


F = diff(A)
F.set_name('F')
F.display()


# As a check, let us compare with Eq. (33.5) of [Misner, Thorne & Wheeler (1973)](https://press.princeton.edu/books/ebook/9781400889099/gravitation):

# In[18]:


F == q/rho2^2 * (r^2-a^2*cos(th)^2)*dr.wedge( dt - a*sin(th)^2*dph ) \
     + 2*q/rho2^2 * a*r*cos(th)*sin(th)*dth.wedge( (r^2+a^2)*dph - a*dt )


# We can get a short expression of $F$ by factoring the components:

# In[19]:


F.apply_map(factor)
F.display()


# The list of non-vanishing components of $F$:

# In[20]:


F.display_comp()


# <p>The Hodge dual of $F$:</p>

# In[21]:


star_F = F.hodge_dual(g)
star_F.display()


# <h3>Maxwell equations</h3>
# 
# <p>Let us check that $F$ obeys the two (source-free) Maxwell equations:</p>

# In[22]:


diff(F).display()


# In[23]:


diff(star_F).display()


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

# In[24]:


nabla = g.connection()
print(nabla)


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

# In[25]:


nabla(g) == 0


# <p>Another view of the above property:</p>

# In[26]:


nabla(g).display()


# <p><span id="cell_outer_13">The nonzero Christoffel symbols (skipping those that can be deduced by symmetry of the last two indices):</span></p>

# In[27]:


g.christoffel_symbols_display()


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

# In[28]:


M.default_frame() is BL.frame()


# In[29]:


BL.frame()


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

# In[30]:


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


# In[31]:


print(xi)


# <p><span id="cell_outer_35">The 1-form associated to it by metric duality is</span></p>

# In[32]:


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


# <p><span id="cell_outer_36">Its covariant derivative is</span></p>

# In[33]:


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


# <p><span id="cell_outer_37">Let us check that the vector field $\xi=\frac{\partial}{\partial t}$ obeys Killing equation:</span></p>

# In[34]:


nab_xi.symmetrize() == 0


# <p><span id="cell_outer_38">Similarly, let us check that</span> $\chi := \frac{\partial}{\partial\phi}$ is a Killing vector:</p>

# In[35]:


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


# In[36]:


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


# <p>Another way to check that $\xi$ and $\chi$ are Killing vectors is the vanishing of the Lie derivative of the metric tensor along them:</p>

# In[37]:


g.lie_derivative(xi) == 0


# In[38]:


g.lie_derivative(chi) == 0


# ## Curvature
# 
# The Ricci tensor of $g$:

# In[39]:


Ric = g.ricci()
print(Ric)


# In[40]:


Ric.display()


# We can get a shorter expression by factorizing the components:

# In[41]:


Ric.apply_map(factor)
Ric.display()


# A matrix view of the components:

# In[42]:


Ric[:]


# <p>Let us check that in the Kerr case, i.e. when $q=0$, the Ricci tensor is zero:</p>

# In[43]:


Ric_Kerr = Ric.copy()
Ric_Kerr.apply_map(lambda f: f.subs({q: 0}))
Ric_Kerr[:]


# ### Riemann tensor
# 
# The Riemann curvature tensor of $g$:

# In[44]:


R = g.riemann()
R.apply_map(factor)
print(R)


# <p>The component $R^0_{\ \, 101}$ of the Riemann tensor is</p>

# In[45]:


R[0,1,0,1]


# <p>The expression in the uncharged limit (Kerr spacetime) is</p>

# In[46]:


R[0,1,0,1].expr().subs(q=0).factor()


# <p>while in the non-rotating limit (Reissner-Nordström spacetime), it is</p>

# In[47]:


R[0,1,0,1].expr().subs(a=0).factor()


# <p>In the Schwarzschild limit, it reduces to</p>

# In[48]:


R[0,1,0,1].expr().subs(a=0, q=0)


# <p>Obviously, it vanishes in the flat space limit:</p>

# In[49]:


R[0,1,0,1].expr().subs(m=0, a=0, q=0)


# ### Ricci scalar
# 
# The Ricci scalar $R = g^{ab} R_{ab}$ of the Kerr-Newman spacetime vanishes identically:

# In[50]:


g.ricci_scalar().display()


# <h2>Einstein equation</h2>
# <p>The Einstein tensor is</p>

# In[51]:


G = Ric - 1/2*g.ricci_scalar()*g
print(G)


# <p>Since the Ricci scalar is zero, the Einstein tensor reduces to the Ricci tensor:</p>

# In[52]:


G == Ric


# The invariant $F_{ab} F^{ab}$ of the electromagnetic field:

# In[53]:


Fuu = F.up(g)
F2 = F['_ab']*Fuu['^ab']
print(F2)


# In[54]:


F2.display()


# <p>The energy-momentum tensor of the electromagnetic field:</p>

# In[55]:


Fud = F.up(g,0)
T = 1/(4*pi)*( F['_k.']*Fud['^k_.'] - 1/4*F2 * g )
T.apply_map(factor)
print(T)


# In[56]:


T[:]


# Check that the Einstein equation is satisfied:

# In[57]:


G == 8*pi*T


# <h3>Bianchi identity</h3>
# 
# <p>Let us check the Bianchi identity $\nabla_p R^i_{\ \, j kl} + \nabla_k R^i_{\ \, jlp} + \nabla_l R^i_{\ \, jpk} = 0$:</p>

# In[58]:


DR = nabla(R)   # long (takes a while)
print(DR)  


# In[59]:


for i in M.irange():
    for j in M.irange():
        for k in M.irange():
            for l in M.irange():
                for p in M.irange():
                    print(DR[i,j,k,l,p] + DR[i,j,l,p,k] + DR[i,j,p,k,l], end=' ')


# <p>If the last sign in the Bianchi identity is changed to minus, the identity does no longer hold:</p>

# In[60]:


DR[0,1,2,3,1] + DR[0,1,3,1,2] + DR[0,1,1,2,3]  # should be zero (Bianchi identity)


# In[61]:


DR[0,1,2,3,1] + DR[0,1,3,1,2] - DR[0,1,1,2,3]  # note the change of the second + to -


# ### Kretschmann scalar
# 
# The tensor $R^\flat$, of components $R_{abcd} = g_{am} R^m_{\ \, bcd}$:

# In[62]:


dR = R.down(g)
print(dR)


# The tensor $R^\sharp$, of components $R^{abcd} = g^{bp} g^{cq} g^{dr} R^a_{\ \, pqr}$:

# In[63]:


uR = R.up(g)
print(uR)


# The Kretschmann scalar $K := R^{abcd} R_{abcd}$:

# In[64]:


Kr_scalar = uR['^abcd']*dR['_abcd']
Kr_scalar.display()


# A variant of this expression can be obtained by invoking the method `factor()` on the coordinate function representing the scalar field in the manifold's default chart:

# In[65]:


Kr = Kr_scalar.coord_function()
Kr.factor()


# <p>As a check, we can compare Kr to the formula given by R. Conn Henry, <a href="http://iopscience.iop.org/0004-637X/535/1/350/fulltext/">Astrophys. J. <strong>535</strong>, 350 (2000)</a>:</p>

# In[66]:


Kr == 8/(r^2+(a*cos(th))^2)^6 *( 
          6*m^2*(r^6 - 15*r^4*(a*cos(th))^2 + 15*r^2*(a*cos(th))^4 - (a*cos(th))^6) 
        - 12*m*q^2*r*(r^4 - 10*(a*r*cos(th))^2 + 5*(a*cos(th))^4) 
        + q^4*(7*r^4 - 34*(a*r*cos(th))^2 + 7*(a*cos(th))^4) )


# <p>The Schwarzschild value of the Kretschmann scalar is recovered by setting $a=0$ and $q=0$:</p>

# In[67]:


Kr.expr().subs(a=0, q=0)


# <p>Let us plot the Kretschmann scalar for $m=1$, $a=0.9$ and $q=0.5$:</p>

# In[68]:


K1 = Kr.expr().subs(m=1, a=0.9, q=0.5)
plot3d(K1, (r,1,3), (th, 0, pi), axes_labels=['r', 'theta', 'Kr'])


# In[69]:


print("Total elapsed time: {} s".format(time.perf_counter() - comput_time0))


# *NB:* most of the computational time is spent in checking the Bianchi identity.