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

# # Schwarzschild spacetime
# 
# This notebook illustrates a few capabilities of SageMath in computations regarding Schwarzschild spacetime. The corresponding tools have been developed within the  [SageManifolds](http://sagemanifolds.obspm.fr) project.
# 
# Click [here](https://raw.githubusercontent.com/sagemanifolds/SageManifolds/master/Notebooks/SM_Schwarzschild.ipynb) to download the notebook file (ipynb format). To run it, you must start SageMath with the Jupyter interface, 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')


# and we initialize a time counter for benchmarking:

# In[3]:


import time
comput_time0 = time.perf_counter()


# ## Spacetime manifold
# 
# We declare the Schwarzschild spacetime as a 4-dimensional Lorentzian manifold:

# In[4]:


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


# In[5]:


print(M)


# <p>The spacetime manifold $\mathcal{M}$ can be split into 4 regions, corresponding to the 4 quadrants in the Kruskal diagram.Let us denote by $\mathcal{R}_{\mathrm{I}}$ to $\mathcal{R}_{\mathrm{IV}}$ the interiors of these 4 regions. $\mathcal{R}_{\mathrm{I}}$ and  $\mathcal{R}_{\mathrm{III}}$  are asymtotically flat regions outside the event horizon;  $\mathcal{R}_{\mathrm{II}}$ is inside the future event horizon and $\mathcal{R}_{\mathrm{IV}}$ is inside the past event horizon.</p>

# In[6]:


regI = M.open_subset('R_I', r'\mathcal{R}_{\mathrm{I}}')
regII = M.open_subset('R_II', r'\mathcal{R}_{\mathrm{II}}')
regIII = M.open_subset('R_III', r'\mathcal{R}_{\mathrm{III}}')
regIV = M.open_subset('R_IV', r'\mathcal{R}_{\mathrm{IV}}')
regI, regII, regIII, regIV


# <p><span id="cell_outer_8">The parameter $m$ of </span><span id="cell_outer_8">the Schwarzschild spacetime is declared as a symbolic variable:</span></p>

# In[7]:


m = var('m') ; assume(m>=0)


# <h2>Boyer-Lindquist coordinates</h2>
# <p>The standard Boyer-Lindquist coordinates (also called Schwarzschild coordinates) are defined on $\mathcal{R}_{\mathrm{I}}\cup \mathcal{R}_{\mathrm{II}}$</p>

# In[8]:


regI_II = regI.union(regII) ; regI_II


# In[9]:


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


# In[10]:


X


# <p>We naturally introduce two subcharts as the restrictions of the chart $X$ to regions $\mathcal{R}_{\mathrm{I}}$ and $\mathcal{R}_{\mathrm{II}}$ respectively. Since, in terms of the Boyer-Lindquist coordinates, $\mathcal{R}_{\mathrm{I}}$ (resp. $\mathcal{R}_{\mathrm{II}}$) is defined by $r>2m$ (resp. $r<2m$), we set</p>

# In[11]:


X_I = X.restrict(regI, r>2*m) ; X_I


# In[12]:


X_II = X.restrict(regII, r<2*m) ; X_II


# <p>At this stage, the manifold's atlas has 3 charts:</p>

# In[13]:


M.atlas()


# In[14]:


M.default_chart()


# <p>Three vector frames have been defined on the manifold: the three coordinate frames:</p>

# In[15]:


M.frames()


# In[16]:


print(M.default_frame())


# In[17]:


M.default_frame().domain()


# ## Metric tensor
# 
# The metric tensor is obtained as follows:

# In[18]:


g = M.metric()
print(g)


# One has to set its components in the coordinate frame associated with Schwarzschild coordinates, which is the current manifold's default frame:

# In[19]:


g[0,0], g[1,1] = -(1-2*m/r), 1/(1-2*m/r)
g[2,2], g[3,3] = r^2, (r*sin(th))^2


# In[20]:


g.display()


# <p>As an example, let us consider a vector field defined only on $\mathcal{R}_{\mathrm{I}}$:</p>

# In[21]:


v = regI.vector_field('v')
v[0] = 1
v[1] = 1 - 2*m/r
# unset components are zero   
v.display()


# In[22]:


v.domain()


# In[23]:


g.domain()


# <p>Since $\mathcal{R}_{\mathrm{I}}\subset \mathcal{M}$, it is possible to apply $g$ to $v$:</p>

# In[24]:


s = g(v,v)
print(s)


# In[25]:


s.display() # v is indeed a null vector


# <div id="cell_text_59" class="text_cell">
# <h2>Levi-Civita Connection</h2>
# </div>
# <div id="insert_new_cell_before27" class="insert_new_cell"> </div>
# <div id="cell_text_27" class="text_cell">
# <p>The Levi-Civita connection $\nabla$ associated with $g$:</p>
# </div>

# In[26]:


nab = g.connection()
print(nab)


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

# In[27]:


nab(g) == 0


# In[28]:


nab(g).display()


# <p>The nonzero Christoffel symbols of $g$ with respect to Schwarzschild coordinates, skipping those that can be deduced by symmetry:</p>

# In[29]:


g.christoffel_symbols_display()


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

# In[30]:


R = g.riemann()
print(R)


# <p>The Weyl conformal tensor associated with $g$:</p>

# In[31]:


C = g.weyl()
print(C)


# In[32]:


C.display()


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

# In[33]:


Ric = g.ricci()
print(Ric)


# <h3>Einstein equation</h3>
# <p>Let us check that the Schwarzschild metric is a solution of the vacuum Einstein equation:</p>

# In[34]:


Ric == 0


# In[35]:


Ric.display() # another view of the above


# <p>Contrary to the Ricci tensor, the Riemann tensor does not vanish:</p>

# In[36]:


R[:]


# In[37]:


R.display()


# <p>The nonzero components of the Riemann tensor, skipping those that can be deduced by antisymmetry:</p>

# In[38]:


R.display_comp(only_nonredundant=True)


# In[39]:


Ric[:]


# <p>Since the Ricci tensor is zero, the Weyl tensor is of course equal to the Riemann tensor:</p>

# In[40]:


C == R


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


DR = nab(R)
print(DR)


# In[42]:


from __future__ import print_function  # because the print below is valid only in Python3
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>Let us check that if we turn the first $+$ sign into a $-$ one, the Bianchi identity does no longer hold:</p>

# In[43]:


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=' ')


# <h3>Kretschmann scalar</h3>
# <p>Let us first introduce tensor $R^\flat$, of components $R_{ijkl} = g_{ip} R^p_{\ \, jkl}$:</p>

# In[44]:


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


# <p>and tensor $R^\sharp$, of components $R^{ijkl} = g^{jp} g^{kq} g^{lr} R^i_{\ \, pqr}$:</p>

# In[45]:


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


# <p>The Kretschmann scalar is $K := R^{ijkl} R_{ijkl}$:</p>

# In[46]:


Kr = 0
for i in M.irange():
    for j in M.irange():
        for k in M.irange():
            for l in M.irange():
                Kr += uR[i,j,k,l]*dR[i,j,k,l]
Kr


# Instead of the above loops, the Kretschmann scalar can also be computed by means of the method `contract()`, asking that the contraction takes place on all indices (positions 0, 1, 2, 3):

# In[47]:


Kr = uR.contract(0, 1, 2, 3, dR, 0, 1, 2, 3)
Kr.expr()


# <p>The contraction can also be performed by means of index notations:</p>

# In[48]:


Kr = uR['^{ijkl}']*dR['_{ijkl}']
Kr.expr()


# <h2>Eddington-Finkelstein coordinates</h2>
# <p>Let us introduce new coordinates on the spacetime manifold: the ingoing Eddington-Finkelstein ones:</p>

# In[49]:


X_EF.<v,r,th,ph> = regI_II.chart(r'v r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\varphi')
print(X_EF); X_EF


# <p>The change from Schwarzschild (Boyer-Lindquist) coordinates to the ingoing Eddington-Finkelstein ones:</p>

# In[50]:


ch_BL_EF_I = X_I.transition_map(X_EF, [t+r+2*m*ln(r/(2*m)-1), r, th, ph], 
                                restrictions2=r>2*m)


# In[51]:


print(ch_BL_EF_I) ; ch_BL_EF_I


# In[52]:


ch_BL_EF_I.display()


# In[53]:


X_EF_I = X_EF.restrict(regI) ; X_EF_I


# In[54]:


ch_BL_EF_II = X_II.transition_map(X_EF, [t+r+2*m*ln(1-r/(2*m)), r, th, ph], 
                                  restrictions2=r<2*m)


# In[55]:


print(ch_BL_EF_II) ; ch_BL_EF_II


# In[56]:


ch_BL_EF_II.display()


# In[57]:


X_EF_II = X_EF.restrict(regII) ; X_EF_II


# <p>The manifold's atlas has now 6 charts:</p>

# In[58]:


M.atlas()


# <p>The default chart is 'BL':</p>

# In[59]:


M.default_chart()


# The change from Eddington-Finkelstein coordinates to the Schwarzschild (Boyer-Lindquist) ones, computed as the inverse of `ch_BL_EF`:

# In[60]:


ch_EF_BL_I = ch_BL_EF_I.inverse()
print(ch_EF_BL_I)


# In[61]:


ch_EF_BL_I.display()


# In[62]:


ch_EF_BL_II = ch_BL_EF_II.inverse()
print(ch_EF_BL_II)


# In[63]:


ch_EF_BL_II.display()


# <p>At this stage, 6 vector frames have been defined on the manifold: the 6 coordinate frames associated with the various charts:</p>

# In[64]:


M.frames()


# <p>The default frame is:</p>

# In[65]:


M.default_frame()


# <p>The coframes are the duals of the defined vector frames:</p>

# In[66]:


M.coframes()


# <p>If not specified, tensor components are assumed to refer to the manifold's default frame. For instance, for the metric tensor:</p>

# In[67]:


g.display()


# In[68]:


g[:]


# The tensor components in the frame associated with Eddington-Finkelstein coordinates in Region I are obtained by providing the frame to the function `display()`:

# In[69]:


g.display(X_EF_I.frame())


# <p>They are also returned by the method comp(), with the frame as argument:</p>

# In[70]:


g.comp(X_EF_I.frame())[:]


# <p>or, as a schortcut,</p>

# In[71]:


g[X_EF_I.frame(),:]


# <p>Similarly, the metric<span id="cell_outer_53"> components in the frame associated with Eddington-Finkelstein coordinates in Region II are obtained by<br /></span></p>

# In[72]:


g.display(X_EF_II.frame())


# <p><span id="cell_outer_56">Note that their form is identical to that in Region I.</span></p>
# 
# <h3>Plot of the Boyer-Lindquist coordinates in terms of the Eddington-Finkelstein ones</h3>
# <p>Let us perform the plot in Region I:</p>

# In[73]:


X_I.plot(X_EF_I, ranges={t:(0,8), r:(2.1,10)}, fixed_coords={th:pi/2,ph:0}, 
         ambient_coords=(r,v), style={t:'--', r:'-'}, parameters={m:1})


# <h2>Tetrad of the static observer</h2>
# <p>Let us introduce the orthonormal tetrad $(e_\alpha)$ associated with the static observers in Schwarzschild spacetime, i.e. the observers whose worldlines are parallel to the timelike Killing vector in the Region I.</p>
# <p>The orthonormal tetrad is defined via a tangent-space automorphism that relates it to the Boyer-Lindquist coordinate frame in Region I:</p>

# In[74]:


ch_to_stat = regI.automorphism_field()
ch_to_stat[0,0], ch_to_stat[1,1] = 1/sqrt(1-2*m/r), sqrt(1-2*m/r)
ch_to_stat[2,2], ch_to_stat[3,3] = 1/r, 1/(r*sin(th))
ch_to_stat[:]


# In[75]:


e = X_I.frame().new_frame(ch_to_stat, 'e')
print(e)


# <p>At this stage, 7 vector frames have been defined on the manifold $\mathcal{M}$:</p>

# In[76]:


M.frames()


# <p>The first vector of the tetrad is the static observer 4-velocity:</p>

# In[77]:


print(e[0])


# In[78]:


e[0].display()


# <p>As any 4-velocity, it is a unit timelike vector:</p>

# In[79]:


g(e[0],e[0]).expr()


# Equivalently, we may use the method `dot` to get the scalar product (since $g$ is the default metric on $\mathcal{M}$):

# In[80]:


e[0].dot(e[0]).expr()


# <p>Let us check that the tetrad $(e_\alpha)$ is orthonormal:</p>

# In[81]:


for i in M.irange():
    for j in M.irange():
        print(e[i].dot(e[j]).expr(), end=' ') 
    print("")


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

# In[82]:


g[e,:]


# <p>or, equivalently,</p>

# In[83]:


g.display(e)


# <p>The expression of the 4-velocity $e_0$ and the vector $e_1$ in terms of the frame associated with Eddington-Finkelstein coordinates:</p>

# In[84]:


e[0].display(X_EF_I.frame())


# In[85]:


e[1].display(X_EF_I.frame())


# <p>Contrary to vectors of a coordinate frame, the vectors of the tetrad $e$ do not commute: their structure coefficients are not identically zero:</p>

# In[86]:


e.structure_coeff()[:]


# <p>Equivalently, the Lie derivative of one vector along another one is not necessarily zero:</p>

# In[87]:


e[0].lie_der(e[1]).display(e)


# <p>The curvature 2-form $\Omega^0_{\ \, 1}$ associated with the tetrad $(e_\alpha)$:</p>

# In[88]:


g.connection().curvature_form(0,1,e).display(X_I.frame())


# <h2>Kruskal-Szekeres coordinates</h2>
# <p>Let us now introduce the Kruskal-Szekeres coordinates $(U,V,\theta,\varphi)$ on the spacetime manifold, via the standard transformation expressing them in terms of the Boyer-Lindquist coordinates $(t,r,\theta,\varphi)$:</p>

# In[89]:


M0 = regI.union(regII).union(regIII).union(regIV) ; M0


# In[90]:


X_KS.<U,V,th,ph> = M0.chart(r'U V th:(0,pi):\theta ph:(0,2*pi):\varphi')
X_KS.add_restrictions(V^2 < 1 + U^2)
X_KS


# In[91]:


X_KS_I = X_KS.restrict(regI, [U>0, V<U, V>-U]) ; X_KS_I


# In[92]:


ch_BL_KS_I = X_I.transition_map(X_KS_I, [sqrt(r/(2*m)-1)*exp(r/(4*m))*cosh(t/(4*m)), 
                                         sqrt(r/(2*m)-1)*exp(r/(4*m))*sinh(t/(4*m)), 
                                         th, ph])
print(ch_BL_KS_I)
ch_BL_KS_I.display()


# In[93]:


X_KS_II = X_KS.restrict(regII, [V>0, V>U, V>-U]) ; X_KS_II


# In[94]:


ch_BL_KS_II = X_II.transition_map(X_KS_II, [sqrt(1-r/(2*m))*exp(r/(4*m))*sinh(t/(4*m)), 
                                            sqrt(1-r/(2*m))*exp(r/(4*m))*cosh(t/(4*m)), 
                                            th, ph])
print(ch_BL_KS_II)
ch_BL_KS_II.display()


# <h3>Plot of the Boyer-Lindquist coordinates in terms of the Kruskal ones</h3>
# <p>We draw the Boyer-Lindquist chart in Region I (red) and Region II (green), with lines of constant $r$ being dashed:</p>

# In[95]:


graphI = X_I.plot(X_KS, ranges={t:(-12,12), r:(2.001,5)}, number_values={t:17, r:9}, 
                  fixed_coords={th:pi/2,ph:0}, ambient_coords=(U,V), 
                  style={t:'--', r:'-'}, parameters={m:1})
graphII = X_II.plot(X_KS, ranges={t:(-12,12), r:(0.001,1.999)}, number_values={t:17, r:9}, 
                    fixed_coords={th:pi/2,ph:0}, ambient_coords=(U,V), 
                    style={t:'--', r:'-'}, color='green', parameters={m:1})
show(graphI+graphII, xmin=-3, xmax=3, ymin=-3, ymax=3, axes_labels=['$U$', '$V$'])


# <p>We may add to the graph the singularity $r=0$ as a Boyer-Lindquist chart plot with $(r,\theta,\phi)$ fixed at $(0, \pi/2, 0)$.  Similarly, we add the event horizon $r=2m$ as a Boyer-Lindquist chart plot with $(r,\theta,\phi)$ fixed at $(2.00001m, \pi/2, 0)$:</p>

# In[96]:


graph_r0 = X_II.plot(X_KS, fixed_coords={r:0, th:pi/2, ph:0}, ambient_coords=(U,V),
                     color='yellow', thickness=3, parameters={m:1})
graph_r2 = X_I.plot(X_KS, ranges={t:(-40,40)}, fixed_coords={r:2.00001, th:pi/2, ph:0},
                    ambient_coords=(U,V), color='black', thickness=2, parameters={m:1})
show(graphI+graphII+graph_r0+graph_r2, xmin=-3, xmax=3, ymin=-3, ymax=3, 
     axes_labels=['$U$', '$V$'])


# <h3>Plot of the Eddington-Finkelstein coordinates in terms of the Kruskal ones</h3>
# <p>We first get the change of coordinates $(v,r,\theta,\phi) \mapsto (U,V,\theta,\phi)$ by composing the change $(v,r,\theta,\phi) \mapsto (t,r,\theta,\phi)$ with $(t,r,\theta,\phi) \mapsto (U,V,\theta,\phi)$: </p>

# In[97]:


ch_EF_KS_I = ch_BL_KS_I * ch_EF_BL_I
ch_EF_KS_I


# In[98]:


ch_EF_KS_I.display()


# In[99]:


ch_EF_KS_II = ch_BL_KS_II * ch_EF_BL_II
ch_EF_KS_II


# In[100]:


graphI_EF = X_EF_I.plot(X_KS, ranges={v:(-12,12), r:(2.001,5)}, number_values={v:17, r:9},
                        fixed_coords={th:pi/2,ph:0}, ambient_coords=(U,V),
                        style={v:'--', r:'-'}, parameters={m:1})
graphII_EF = X_EF_II.plot(X_KS, ranges={v:(-12,12), r:(0.001,1.999)}, 
                          number_values={v:17, r:9},
                          fixed_coords={th:pi/2,ph:0}, ambient_coords=(U,V),
                          style={v:'--', r:'-'}, color='green', parameters={m:1})
show(graphI_EF+graphII_EF+graph_r0+graph_r2, xmin=-3, xmax=3, ymin=-3, ymax=3,
     axes_labels=['$U$', '$V$'])


# <p>There are now 9 charts defined on the spacetime manifold:</p>

# In[101]:


M.atlas()


# In[102]:


len(M.atlas())


# <p>There are 8 explicit coordinate changes (the coordinate change KS $\rightarrow$ BL is not known in explicit form):</p>

# In[103]:


M.coord_changes()


# In[104]:


len(M.coord_changes())


# <p>There are 10 vector frames (among which 9 coordinate frames):</p>

# In[105]:


M.frames()


# In[106]:


len(M.frames())


# <p>There are 14 fields of tangent space automorphisms expressing the changes of coordinate bases and tetrad:</p>

# In[107]:


len(M.changes_of_frame())


# Thanks to these changes of frames, the components of the metric tensor with respect to the Kruskal-Szekeres can be computed by the method `display()` and are found to be:

# In[108]:


g.display(X_KS_I.frame())


# In[109]:


g[X_KS_I.frame(),:]


# In[110]:


g.display(X_KS_II.frame())


# <p>The first vector of the orthonormal tetrad $e$ expressed on the <span id="cell_outer_86"><span id="cell_outer_121">Kruskal-Szekeres</span></span> frame:</p>

# In[111]:


e[0].display(X_KS_I.frame())


# <p>The Riemann curvature tensor in terms of the <span id="cell_outer_121">Kruskal-Szekeres</span> coordinates:</p>

# In[112]:


g.riemann().display(X_KS_I.frame())


# In[113]:


g.riemann().display_comp(X_KS_I.frame(), only_nonredundant=True)


# <p>The curvature 2-form $\Omega^0_{\ \, 1}$ associated to the Kruskal-Szekeres coordinate frame:</p>

# In[114]:


om = g.connection().curvature_form(0,1, X_KS_I.frame())
print(om)


# In[115]:


om.display(X_KS_I.frame())


# <h2>Isotropic coordinates</h2>
# <p>Let us now introduce isotropic coordinates $(t,\bar{r},\theta,\varphi)$ on the spacetime manifold:</p>

# In[116]:


regI_III = regI.union(regIII) ; regI_III


# In[117]:


X_iso.<t,ri,th,ph> = regI_III.chart(r't ri:(0,+oo):\bar{r} th:(0,pi):\theta ph:(0,2*pi):\varphi')
print(X_iso) ; X_iso


# In[118]:


X_iso_I = X_iso.restrict(regI, ri>m/2) ; X_iso_I


# <p>The transformation from the isotropic coordinates to the Boyer-Lindquist ones:</p>

# In[119]:


assume(2*ri>m) # we consider only region I
ch_iso_BL_I = X_iso_I.transition_map(X_I, [t, ri*(1+m/(2*ri))^2, th, ph])
print(ch_iso_BL_I)
ch_iso_BL_I.display()


# In[120]:


assume(r>2*m) # we consider only region I
ch_iso_BL_I.set_inverse(t, (r-m+sqrt(r*(r-2*m)))/2, th, ph, verbose=True)


# In[121]:


ch_iso_BL_I.inverse().display()


# <p><span id="cell_outer_90">At this stage, 11 charts have been defined on the manifold $\mathcal{M}$:</span></p>

# In[122]:


M.atlas()


# In[123]:


len(M.atlas())


# <p>12 vector frames have been defined on $\mathcal{M}$: 11 coordinate bases and the tetrad $(e_\alpha)$:</p>

# In[124]:


M.frames()


# In[125]:


len(M.frames())


# <p>The components of the metric tensor in terms of the isotropic coordinates are given by</p>

# In[126]:


g.display(X_iso_I.frame(), X_iso_I)


# <p>The $g_{00}$ component can be factorized:</p>

# In[127]:


g[X_iso_I.frame(), 0,0, X_iso_I]


# In[128]:


g[X_iso_I.frame(), 0,0, X_iso_I].factor()


# Let us also factor the other components:

# In[129]:


for i in range(1,4):
    g[X_iso_I.frame(), i,i, X_iso_I].factor()


# The output of the `display()` command looks nicer:

# In[130]:


g.display(X_iso_I.frame(), X_iso_I)


# <p>Expression of the tetrad associated with the static observer in terms of the isotropic coordinate basis:</p>

# In[131]:


e[0].display(X_iso_I.frame(), X_iso_I)


# In[132]:


e[1].display(X_iso_I.frame(), X_iso_I)


# In[133]:


e[2].display(X_iso_I.frame(), X_iso_I)


# In[134]:


e[3].display(X_iso_I.frame(), X_iso_I)


# In[135]:


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