Schwarzschild spacetime

This worksheet demonstrates a few capabilities of SageManifolds (version 1.0, as included in SageMath 7.5) in computations regarding Schwarzschild spacetime.

Click here 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()
Out[1]:
'SageMath version 7.5, Release Date: 2017-01-11'

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

In [2]:
%display latex

Spacetime manifold

We declare the Schwarzschild spacetime as a 4-dimensional differentiable manifold:

In [3]:
M = Manifold(4, 'M', r'\mathcal{M}') ; M
Out[3]:
In [4]:
print(M)
4-dimensional differentiable manifold M

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.

In [5]:
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
Out[5]:

The parameter $m$ of the Schwarzschild spacetime is declared as a symbolic variable:

In [6]:
m = var('m') ; assume(m>=0)

Boyer-Lindquist coordinates

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

In [7]:
regI_II = regI.union(regII) ; regI_II
Out[7]:
In [8]:
X.<t,r,th,ph> = regI_II.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')
print(X)
Chart (R_I_union_R_II, (t, r, th, ph))
In [9]:
X
Out[9]:

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

In [10]:
X_I = X.restrict(regI, r>2*m) ; X_I
Out[10]:
In [11]:
X_II = X.restrict(regII, r<2*m) ; X_II
Out[11]:

At this stage, the manifold's atlas has 3 charts:

In [12]:
M.atlas()
Out[12]:
In [13]:
M.default_chart()
Out[13]:

Three vector frames have been defined on the manifold: the three coordinate frames:

In [14]:
M.frames()
Out[14]:
In [15]:
print(M.default_frame())
Coordinate frame (R_I_union_R_II, (d/dt,d/dr,d/dth,d/dph))
In [16]:
M.default_frame().domain()
Out[16]:

Metric tensor

The metric tensor is defined as follows:

In [17]:
g = M.lorentzian_metric('g')
print(g)
Lorentzian metric g on the 4-dimensional differentiable manifold M

The metric tensor is set by its components in the coordinate frame associated with Schwarzschild coordinates, which is the current manifold's default frame:

In [18]:
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 [19]:
g.display()
Out[19]:

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

In [20]:
v = regI.vector_field('v')
v[0] = 1
v[1] = 1 - 2*m/r
# unset components are zero   
v.display()
Out[20]:
In [21]:
v.domain()
Out[21]:
In [22]:
g.domain()
Out[22]:

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

In [23]:
s = g(v,v) ; print(s)
Scalar field g(v,v) on the Open subset R_I of the 4-dimensional differentiable manifold M
In [24]:
s.display() # v is indeed a null vector
Out[24]:

Levi-Civita Connection

The Levi-Civita connection $\nabla$ associated with $g$:

In [25]:
nab = g.connection() ; print(nab)
Levi-Civita connection nabla_g associated with the Lorentzian metric g on the 4-dimensional differentiable manifold M

Let us verify that the covariant derivative of $g$ with respect to $\nabla$ vanishes identically:

In [26]:
nab(g) == 0
Out[26]:
In [27]:
nab(g).display()
Out[27]:

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

In [28]:
g.christoffel_symbols_display()
Out[28]:

Curvature

The Riemann curvature tensor associated with $g$:

In [29]:
R = g.riemann() ; print(R)
Tensor field Riem(g) of type (1,3) on the 4-dimensional differentiable manifold M

The Weyl conformal tensor associated with $g$:

In [30]:
C = g.weyl() ; print(C)
Tensor field C(g) of type (1,3) on the 4-dimensional differentiable manifold M
In [31]:
C.display()
Out[31]:

The Ricci tensor associated with $g$:

In [32]:
Ric = g.ricci() ; print(Ric)
Field of symmetric bilinear forms Ric(g) on the 4-dimensional differentiable manifold M

Einstein equation

Let us check that the Schwarzschild metric is a solution of the vacuum Einstein equation:

In [33]:
Ric == 0
Out[33]:
In [34]:
Ric.display() # another view of the above
Out[34]:

Contrary to the Ricci tensor, the Riemann tensor does not vanish:

In [35]:
R[:]
Out[35]:
In [36]:
R.display()
Out[36]:

The nonzero components of the Riemann tensor, skipping those that can be deduced by antisymmetry:

In [37]:
R.display_comp(only_nonredundant=True)
Out[37]:
In [38]:
Ric[:]
Out[38]:

Since the Ricci tensor is zero, the Weyl tensor is of course equal to the Riemann tensor:

In [39]:
C == R
Out[39]:

Bianchi identity

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

In [40]:
DR = nab(R) ; print(DR)
Tensor field nabla_g(Riem(g)) of type (1,4) on the 4-dimensional differentiable manifold M
In [41]:
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] ,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Let us check that if we turn the first $+$ sign into a $-$ one, the Bianchi identity does no longer hold:

In [42]:
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] ,
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -12*m/(2*m*r^3 - r^4) 0 0 12*m/(2*m*r^3 - r^4) 0 0 0 0 0 0 0 0 0 0 0 0 0 -6*m/r^2 0 0 0 0 0 6*m/r^2 0 0 0 0 0 0 0 0 0 0 -6*m*sin(th)^2/r^2 0 0 0 0 0 0 0 0 6*m*sin(th)^2/r^2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -6*m/r^2 0 0 0 0 0 6*m/r^2 0 0 0 0 0 0 0 0 -6*m/r^2 0 0 6*m/r^2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -6*m*sin(th)^2/r^2 0 0 0 0 0 0 0 0 6*m*sin(th)^2/r^2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -6*m*sin(th)^2/r^2 0 0 6*m*sin(th)^2/r^2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -12*(2*m^2 - m*r)/r^5 0 0 12*(2*m^2 - m*r)/r^5 0 0 0 0 0 0 0 0 0 0 0 0 0 -6*(4*m^3 - 4*m^2*r + m*r^2)/r^4 0 0 0 0 0 6*(4*m^3 - 4*m^2*r + m*r^2)/r^4 0 0 0 0 0 0 0 0 0 0 -6*(4*m^3 - 4*m^2*r + m*r^2)*sin(th)^2/r^4 0 0 0 0 0 0 0 0 6*(4*m^3 - 4*m^2*r + m*r^2)*sin(th)^2/r^4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -6*m/r^2 0 0 6*m/r^2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6*(2*m^2 - m*r)*sin(th)^2/r 0 0 -6*(2*m^2 - m*r)*sin(th)^2/r 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -6*m*sin(th)^2/r^2 0 0 0 0 0 6*m*sin(th)^2/r^2 0 0 0 0 0 0 0 0 0 0 0 0 0 -6*(2*m^2 - m*r)*sin(th)^2/r 0 0 6*(2*m^2 - m*r)*sin(th)^2/r 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6*(2*m^2 - m*r)/r^5 0 0 0 0 0 -6*(2*m^2 - m*r)/r^5 0 0 0 0 0 0 0 0 6*(2*m^2 - m*r)/r^5 0 0 -6*(2*m^2 - m*r)/r^5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -6*m/(2*m*r^3 - r^4) 0 0 6*m/(2*m*r^3 - r^4) 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6*m*sin(th)^2/r^2 0 0 -6*m*sin(th)^2/r^2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 12*m*sin(th)^2/r^2 0 0 -12*m*sin(th)^2/r^2 0 0 0 0 0 0 0 0 6*m*sin(th)^2/r^2 0 0 0 0 0 -6*m*sin(th)^2/r^2 0 0 0 0 0 0 0 0 -6*m*sin(th)^2/r^2 0 0 6*m*sin(th)^2/r^2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6*(2*m^2 - m*r)/r^5 0 0 0 0 0 0 0 0 -6*(2*m^2 - m*r)/r^5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6*(2*m^2 - m*r)/r^5 0 0 -6*(2*m^2 - m*r)/r^5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -6*m/(2*m*r^3 - r^4) 0 0 0 0 0 6*m/(2*m*r^3 - r^4) 0 0 0 0 0 0 0 0 0 0 0 0 0 -6*m/r^2 0 0 6*m/r^2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -12*m/r^2 0 0 12*m/r^2 0 0 0 0 0 0 0 0 -6*m/r^2 0 0 0 0 0 6*m/r^2 0 0 0 0 0 0 0 0 6*m/r^2 0 0 -6*m/r^2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Kretschmann scalar

Let us first introduce tensor $R^\flat$, of components $R_{ijkl} = g_{ip} R^p_{\ \, jkl}$:

In [43]:
dR = R.down(g) ; print(dR)
Tensor field of type (0,4) on the 4-dimensional differentiable manifold M

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

In [44]:
uR = R.up(g) ; print(uR)
Tensor field of type (4,0) on the 4-dimensional differentiable manifold M

The Kretschmann scalar is $K := R^{ijkl} R_{ijkl}$:

In [45]:
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
Out[45]:

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

In [46]:
Kr = uR.contract(0, 1, 2, 3, dR, 0, 1, 2, 3)
Kr.expr()
Out[46]:

The contraction can also be performed by means of index notations:

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

Eddington-Finkelstein coordinates

Let us introduce new coordinates on the spacetime manifold: the ingoing Eddington-Finkelstein ones:

In [48]:
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
Chart (R_I_union_R_II, (v, r, th, ph))
Out[48]:

The change from Schwarzschild (Boyer-Lindquist) coordinates to the ingoing Eddington-Finkelstein ones:

In [49]:
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 [50]:
print(ch_BL_EF_I) ; ch_BL_EF_I
Change of coordinates from Chart (R_I, (t, r, th, ph)) to Chart (R_I, (v, r, th, ph))
Out[50]:
In [51]:
ch_BL_EF_I.display()
Out[51]:
In [52]:
X_EF_I = X_EF.restrict(regI) ; X_EF_I
Out[52]:
In [53]:
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 [54]:
print(ch_BL_EF_II) ; ch_BL_EF_II
Change of coordinates from Chart (R_II, (t, r, th, ph)) to Chart (R_II, (v, r, th, ph))
Out[54]:
In [55]:
ch_BL_EF_II.display()
Out[55]:
In [56]:
X_EF_II = X_EF.restrict(regII) ; X_EF_II
Out[56]:

The manifold's atlas has now 6 charts:

In [57]:
M.atlas()
Out[57]:

The default chart is 'BL':

In [58]:
M.default_chart()
Out[58]:

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

In [59]:
ch_EF_BL_I = ch_BL_EF_I.inverse() ; print(ch_EF_BL_I)
Change of coordinates from Chart (R_I, (v, r, th, ph)) to Chart (R_I, (t, r, th, ph))
In [60]:
ch_EF_BL_I.display()
Out[60]:
In [61]:
ch_EF_BL_II = ch_BL_EF_II.inverse() ; print(ch_EF_BL_II)
Change of coordinates from Chart (R_II, (v, r, th, ph)) to Chart (R_II, (t, r, th, ph))
In [62]:
ch_EF_BL_II.display()
Out[62]:

At this stage, 6 vector frames have been defined on the manifold: the 6 coordinate frames associated with the various charts:

In [63]:
M.frames()
Out[63]:

The default frame is:

In [64]:
M.default_frame()
Out[64]:

The coframes are the duals of the defined vector frames:

In [65]:
M.coframes()
Out[65]:

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

In [66]:
g.display()
Out[66]:
In [67]:
g[:]
Out[67]:

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 [68]:
g.display(X_EF_I.frame())
Out[68]:

They are also returned by the method comp(), with the frame as argument:

In [69]:
g.comp(X_EF_I.frame())[:]
Out[69]:

or, as a schortcut,

In [70]:
g[X_EF_I.frame(),:]
Out[70]:

Similarly, the metric components in the frame associated with Eddington-Finkelstein coordinates in Region II are obtained by

In [71]:
g.display(X_EF_II.frame())
Out[71]:

Note that their form is identical to that in Region I.

Plot of the Boyer-Lindquist coordinates in terms of the Eddington-Finkelstein ones

Let us perform the plot in Region I:

In [72]:
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})
Out[72]:

Tetrad of the static observer

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.

The orthonormal tetrad is defined via a tangent-space automorphism that relates it to the Boyer-Lindquist coordinate frame in Region I:

In [73]:
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[:]
Out[73]:
In [74]:
e = X_I.frame().new_frame(ch_to_stat, 'e') ; print(e)
Vector frame (R_I, (e_0,e_1,e_2,e_3))

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

In [75]:
M.frames()
Out[75]:

The first vector of the tetrad is the static observer 4-velocity:

In [76]:
print(e[0])
Vector field e_0 on the Open subset R_I of the 4-dimensional differentiable manifold M
In [77]:
e[0].display()
Out[77]:

As any 4-velocity, it is a unit timelike vector:

In [78]:
g(e[0],e[0]).expr()
Out[78]:

Let us check that the tetrad $(e_\alpha)$ is orthonormal:

In [79]:
for i in M.irange():
    for j in M.irange():
        print g(e[i],e[j]).expr() , 
    print " "
-1 0 0 0  
0 1 0 0  
0 0 1 0  
0 0 0 1  

Another view of the above result:

In [80]:
g[e,:]
Out[80]:

or, equivalently,

In [81]:
g.display(e)
Out[81]:

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

In [82]:
e[0].display(X_EF_I.frame())
Out[82]:
In [83]:
e[1].display(X_EF_I.frame())
Out[83]:

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

In [84]:
e.structure_coeff()[:]
Out[84]:

Equivalently, the Lie derivative of one vector along another one is not necessarily zero:

In [85]:
e[0].lie_der(e[1]).display(e)
Out[85]:

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

In [86]:
g.connection().curvature_form(0,1,e).display(X_I.frame())
Out[86]:

Kruskal-Szekeres coordinates

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)$:

In [87]:
M0 = regI.union(regII).union(regIII).union(regIV) ; M0
Out[87]:
In [88]:
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
Out[88]:
In [89]:
X_KS_I = X_KS.restrict(regI, [U>0, V<U, V>-U]) ; X_KS_I
Out[89]:
In [90]:
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()
Change of coordinates from Chart (R_I, (t, r, th, ph)) to Chart (R_I, (U, V, th, ph))
Out[90]:
In [91]:
X_KS_II = X_KS.restrict(regII, [V>0, V>U, V>-U]) ; X_KS_II
Out[91]:
In [92]:
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()
Change of coordinates from Chart (R_II, (t, r, th, ph)) to Chart (R_II, (U, V, th, ph))
Out[92]:

Plot of the Boyer-Lindquist coordinates in terms of the Kruskal ones

We draw the Boyer-Lindquist chart in Region I (red) and Region II (green), with lines of constant $r$ being dashed:

In [93]:
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$'])

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)$:

In [94]:
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$'])

Plot of the Eddington-Finkelstein coordinates in terms of the Kruskal ones

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)$:

In [95]:
ch_EF_KS_I = ch_BL_KS_I * ch_EF_BL_I
ch_EF_KS_I
Out[95]:
In [96]:
ch_EF_KS_I.display()
Out[96]:
In [97]:
ch_EF_KS_II = ch_BL_KS_II * ch_EF_BL_II
ch_EF_KS_II
Out[97]:
In [98]:
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$'])

There are now 9 charts defined on the spacetime manifold:

In [99]:
M.atlas()
Out[99]:
In [100]:
len(M.atlas())
Out[100]:

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

In [101]:
M.coord_changes()
Out[101]:
In [102]:
len(M.coord_changes())
Out[102]:

There are 10 vector frames (among which 9 coordinate frames):

In [103]:
M.frames()
Out[103]:
In [104]:
len(M.frames())
Out[104]:

There are 14 fields of tangent space automorphisms expressing the changes of coordinate bases and tetrad:

In [105]:
len(M.changes_of_frame())
Out[105]:

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 [106]:
g.display(X_KS_I.frame())
Out[106]:
In [107]:
g[X_KS_I.frame(),:]
Out[107]:
In [108]:
g.display(X_KS_II.frame())
Out[108]:

The first vector of the orthonormal tetrad $e$ expressed on the Kruskal-Szekeres frame:

In [109]:
e[0].display(X_KS_I.frame())
Out[109]:

The Riemann curvature tensor in terms of the Kruskal-Szekeres coordinates:

In [110]:
g.riemann().display(X_KS_I.frame())
Out[110]:
In [111]:
g.riemann().display_comp(X_KS_I.frame(), only_nonredundant=True)
Out[111]:

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

In [112]:
om = g.connection().curvature_form(0,1, X_KS_I.frame()) ; print(om)
2-form curvature (0,1) of connection nabla_g w.r.t. Coordinate frame (R_I, (d/dU,d/dV,d/dth,d/dph)) on the Open subset R_I of the 4-dimensional differentiable manifold M
In [113]:
om.display(X_KS_I.frame())
Out[113]:

Isotropic coordinates

Let us now introduce isotropic coordinates $(t,\bar{r},\theta,\varphi)$ on the spacetime manifold:

In [114]:
regI_III = regI.union(regIII) ; regI_III
Out[114]:
In [115]:
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
Chart (R_I_union_R_III, (t, ri, th, ph))
Out[115]:
In [116]:
X_iso_I = X_iso.restrict(regI, ri>m/2) ; X_iso_I
Out[116]:

The transformation from the isotropic coordinates to the Boyer-Lindquist ones:

In [117]:
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()
Change of coordinates from Chart (R_I, (t, ri, th, ph)) to Chart (R_I, (t, r, th, ph))
Out[117]:
In [118]:
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)
Check of the inverse coordinate transformation:
  t == t
  ri == ri
  th == th
  ph == ph
  t == t
  r == r
  th == th
  ph == ph
In [119]:
ch_iso_BL_I.inverse().display()
Out[119]:

At this stage, 11 charts have been defined on the manifold $\mathcal{M}$:

In [120]:
M.atlas()
Out[120]:
In [121]:
len(M.atlas())
Out[121]:

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

In [122]:
M.frames()
Out[122]:
In [123]:
len(M.frames())
Out[123]:

The components of the metric tensor in terms of the isotropic coordinates are given by

In [124]:
g.display(X_iso_I.frame(), X_iso_I)
Out[124]:

The $g_{00}$ component can be factorized:

In [125]:
g[X_iso_I.frame(), 0,0, X_iso_I]
Out[125]:
In [126]:
g[X_iso_I.frame(), 0,0, X_iso_I].factor()
Out[126]:

Let us also factorize the other components:

In [127]:
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 [128]:
g.display(X_iso_I.frame(), X_iso_I)
Out[128]:

Expression of the tetrad associated with the static observer in terms of the isotropic coordinate basis:

In [129]:
e[0].display(X_iso_I.frame(), X_iso_I)
Out[129]:
In [130]:
e[1].display(X_iso_I.frame(), X_iso_I)
Out[130]:
In [131]:
e[2].display(X_iso_I.frame(), X_iso_I)
Out[131]:
In [132]:
e[3].display(X_iso_I.frame(), X_iso_I)
Out[132]: