Python (SymPy) Programming for finding Schwarzschild Solution to Einstein Equation

This article demonstrates a way to check if the Schwarzschild metrc satisfies the Einstein equation using a Python library "SymPy". In general, in order to examine the Einstein equation, one needs to compute various objects, for instance, partial differentials of the spacetime metric, the Christoffel symbols, Riemann curvature tensors, Ricci tensors and so on. It will be of great benefit if you could use your computer for such symbolic calculations and a Python's library "SymPy" is such a one.

The following article describes some basic usage for SymPy:

pianofisica.hatenablog.com

This article mainly takes a look into 4-dimensional case, but also 5-dim, and without showing the code, 6-dim. cases as well. It is an advantage of using computer that you can extend the situation easily.　

(Pseudo-)Riemannian Geometry for 4-dimensional Spacetime
Schwarzschild Sol. in 4D
- Einstein Equation
Schwarzschild Sol. in 5D
Schwarzschild Sol. in 6D

(Pseudo-)Riemannian Geometry for 4-dimensional Spacetime

Hereafter, a pair of the same indices appearing in the same term on lower/upper is to be summed over 0 to 3 unless otherwise mentioned (Einstein convention).

Coordinates

A coordinate system

$\qquad\displaystyle{\xi^i=\left(\begin{array}{c} \xi_0 \\ \xi_1 \\ \xi_2 \\ \xi_3 \end{array} \right)=\left(\begin{array}{c} t \\ r \\ \theta \\ \phi \end{array} \right)}$

is implemented in SymPy as

import sympy as sp

sp.var('t, r, theta, phi')

xi = sp.Matrix([ [t],[r],[theta],[phi] ])

for the later convenience to extend, let D be the dimension of spacetime:

D = len(xi)

Spacetime Metric

The spacetime metric $g_{ij}$ associated with the coordinate system $\{x^i\}$ defines an invariant distance $ds^2$ via

$\qquad\displaystyle{ds^2=g_{ij}\,d\xi^i d\xi^j }$

Since the metric is a 2nd rank tensor, it can be represented as a matirx

$\qquad\displaystyle{g_{ij}=\left(\begin{array}{cccc} g_{00} & g_{01} & g_{02} & g_{03} \\ g_{10} & g_{11} & g_{12} & g_{13} \\ g_{20} & g_{21} & g_{22} & g_{23} \\ g_{30} & g_{31} & g_{32} & g_{33} \end{array} \right)}$
　
Here assume the metric is given in the following form

$\qquad\displaystyle{ds^2=-f(r)dt^2 + \frac{1}{f(r)}dr^2 + r^2d\theta^2 + r^2\sin^2\theta d\phi^2 }$

That is, their non-trivial components read

$\qquad\displaystyle{g_{tt}=-f(r)\qquad g_{rr}=\frac{1}{f(r)} \qquad g_{\theta\theta} = r^2 \qquad g_{\phi\phi} = r^2\sin^2\theta }$

where $f=f(r)$ is a function of $r$ and its expression shall be determined by the Einstein equation later. The metric is input by

f = sp.Function('f')(r)

g = sp.Matrix([
[-f, 0, 0, 0],
[0, 1/f, 0, 0],
[0, 0, r**2, 0],
[0, 0, 0, r**2*sp.sin(theta)**2] ])

The metric inverse $g^{ij}$ that satisfies $g^{ik}g_{kj}=\delta^i_j$ is obtained by

ginv = g.inv()

Christoffel Symbols

The Christoffel symbols which specifies the Levi-Civita connection is given by

$\qquad\displaystyle{\Gamma^i_{jk}=\frac12g^{in}\left(\frac{\partial g_{nj}}{\partial x^k}+\frac{\partial g_{nk}}{\partial x^j}-\frac{\partial g_{jk}}{\partial x^n}\right)}$

which can be computed by

def Gamma(k,i,j):
    gm = 0
    for n in range(D):
        gm = gm + sp.Rational(1,2)*ginv[k,n]*(g[n,i].diff(xi[j])\
                +g[n,j].diff(xi[i])-g[i,j].diff(xi[n]))
    return gm

for example

Gamma(1,3,3)

gives

$\qquad\displaystyle{\Gamma^1_{33}=-rf(r)\sin^2\theta}$
　

Riemann Curvature Tensor

Riemann curvature is defined by

$\qquad\displaystyle{R^i_{jkl}=\frac{\partial\Gamma^i_{lj}}{\partial x^k}-\frac{\partial\Gamma^i_{kj}}{\partial x^l}+\Gamma^n_{lj}\Gamma^i_{kn}-\Gamma^n_{kj}\Gamma^i_{ln}}$

def R(i,j,k,l):
    rm = Gamma(i,l,j).diff(xi[k]) - Gamma(i,k,j).diff(xi[l]) 
    for n in range(D):
        rm = rm + Gamma(n,l,j)*Gamma(i,k,n) - Gamma(n,k,j)*Gamma(i,l,n)
    return rm

for instance

R(1,3,1,3)

shows

$\qquad\displaystyle{R^1_{313}=-\frac{r}{2}f'(r)\sin^2\theta}$

Ricci Tensor

Ricci tensor is defined by

$\qquad\displaystyle{R_{jk}=R^n_{jnk}}$

which shall be computed by

def Ricc(j,k):
    ric = 0
    for n in range(D):
        ric = ric + R(n,j,n,k)
    return ric

Since this is a 2nd-rank tensor, let it be represented as a matrix

RiccMat = sp.Matrix([
[Ricc(0,0), Ricc(0,1), Ricc(0,2), Ricc(0,3)],
[Ricc(1,0), Ricc(1,1), Ricc(1,2), Ricc(1,3)],
[Ricc(2,0), Ricc(2,1), Ricc(2,2), Ricc(2,3)],
[Ricc(3,0), Ricc(3,1), Ricc(3,2), Ricc(3,3)] ])
sp.simplify(RiccMat)

then one shall find

${R_{jk}=\left(\begin{array}{cccc}\frac{\left(r \frac{d^{2}}{d r^{2}} f{\left(r \right)} + 2 \frac{d}{d r} f{\left(r \right)}\right) f{\left(r \right)}}{2 r} & 0 & 0 & 0\\0 & \frac{- \frac{r \frac{d^{2}}{d r^{2}} f{\left(r \right)}}{2} - \frac{d}{d r} f{\left(r \right)}}{r f{\left(r \right)}} & 0 & 0\\0 & 0 & - r \frac{d}{d r} f{\left(r \right)} - f{\left(r \right)} + 1 & 0\\0 & 0 & 0 & \left(- r \frac{d}{d r} f{\left(r \right)} - f{\left(r \right)} + 1\right) \sin^{2}{\left(\theta \right)}\end{array}\right)}$

Ricci Scalar

Ricci scalar is an invariant given by

$\qquad\displaystyle{R=g^{jk}R_{jk}}$

R = sp.simplify((ginv*RiccMat).trace())
R

results in

$\qquad\displaystyle{R=\frac{- r^{2} \frac{d^{2}}{d r^{2}} f{\left(r \right)} - 4 r \frac{d}{d r} f{\left(r \right)} - 2 f{\left(r \right)} + 2}{r^{2}}}$

Einstein Tensor

Einstein tensor is a combination of the quantities mentioned above:

$\qquad\displaystyle{G_{ij}=R_{ij}-\frac12Rg_{ij}}$

GMat = sp.simplify(RiccMat - sp.Rational(1,2)*g*R)
GMat

shows

$\displaystyle{G_{jk}=\left(\begin{array}{cccc}\frac{\left(- r \frac{d}{d r} f{\left(r \right)} - f{\left(r \right)} + 1\right) f{\left(r \right)}}{r^{2}} & 0 & 0 & 0\\0 & \frac{r \frac{d}{d r} f{\left(r \right)} + f{\left(r \right)} - 1}{r^{2} f{\left(r \right)}} & 0 & 0\\0 & 0 & \frac{r \left(r \frac{d^{2}}{d r^{2}} f{\left(r \right)} + 2 \frac{d}{d r} f{\left(r \right)}\right)}{2} & 0\\0 & 0 & 0 & \frac{r \left(r \frac{d^{2}}{d r^{2}} f{\left(r \right)} + 2 \frac{d}{d r} f{\left(r \right)}\right) \sin^{2}{\left(\theta \right)}}{2}\end{array}\right)}$
　

Schwarzschild Sol. in 4D

Einstein Equation

The fundamental equation for determining the spacetime metric $g_{ij}$ is the Einstein equation

$\qquad\displaystyle{G_{ij}=0}$

where no matter-fields neither cosmological constant is assumed. The equation requires for the unknown function $f=f(r)$ to satisfy

$\qquad\displaystyle{rf''+2f'=0\\ -rf'-f+1=0}$

Since the first equation automatically follows from the second, only the second equation is to be considered:

$\qquad\displaystyle{\frac{f'}{1-f}=\frac{1}{r}}$

whose solution is easily found as

$\qquad\displaystyle{f=1-\frac{a}{r}}$

where $a$ is a constant. The vacuum solution to the Einstein equation given by this function $f$ is known as the Schwarzschild solution. One of the feature of this solution is that it shows a singularity at $r=a$ whose radius is called Schwarzschild radius.

Schwarzschild Sol. in 5D

Let it apply the same calculation to 5-dim. spacetime with coordinate system $\xi=(t,r,\theta,\phi,\psi)$ whose metric shall take the form

$\displaystyle{g_{ij}=\left(\begin{array}{ccccc} g_{tt} & g_{tr} & g_{t\theta} & g_{t\phi} & g_{t\psi} \\ g_{rt} & g_{rr} & g_{r\theta} & g_{r\phi} & g_{r\psi}\\ g_{\theta t} & g_{\theta r} & g_{\theta\theta} & g_{\theta\phi} & g_{\theta\psi} \\ g_{\phi t} & g_{\phi r} & g_{\phi\theta} & g_{\phi\phi} & g_{\phi\psi} \\ g_{\psi t} & g_{\psi r} & g_{\psi\theta} & g_{\psi\phi} & g_{\psi\psi} \end{array} \right)=\left(\begin{array}{ccccc} -f(r) & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{f(r)} & 0 & 0 & 0 \\ 0 & 0 & r^2 & 0 & 0 \\ 0 & 0 & 0 & r^2\sin^2\theta &0 \\ 0 & 0 & 0 &0 & r^2\sin^2\theta\sin^2\phi \end{array} \right)}$

%reset
import sympy as sp

sp.var('t, r, theta, phi, psi')
xi = sp.Matrix([ [t],[r],[theta],[phi],[psi] ])
D = len(xi)

f = sp.Function('f')(r)
g = sp.Matrix([
[-f, 0, 0, 0, 0],
[0, 1/f, 0, 0, 0],
[0, 0, r**2, 0, 0],
[0, 0, 0, r**2*sp.sin(theta)**2, 0],
[0, 0, 0, 0, r**2*sp.sin(theta)**2*sp.sin(phi)**2] ])

ginv = g.inv()

def Gamma(k,i,j):
    gm = 0
    for n in range(D):
        gm = gm + sp.Rational(1,2)*ginv[k,n]*(g[n,i].diff(xi[j])\
                +g[n,j].diff(xi[i])-g[i,j].diff(xi[n]))
    return gm

def R(i,j,k,l):
    rm = Gamma(i,l,j).diff(xi[k]) - Gamma(i,k,j).diff(xi[l]) 
    for n in range(D):
        rm = rm + Gamma(n,l,j)*Gamma(i,k,n) - Gamma(n,k,j)*Gamma(i,l,n)
    return rm

def Ricc(j,k):
    ric = 0
    for n in range(D):
        ric = ric + R(n,j,n,k)
    return ric

RiccMat = sp.Matrix([
[Ricc(0,0), Ricc(0,1), Ricc(0,2), Ricc(0,3), Ricc(0,4)],
[Ricc(1,0), Ricc(1,1), Ricc(1,2), Ricc(1,3), Ricc(1,4)],
[Ricc(2,0), Ricc(2,1), Ricc(2,2), Ricc(2,3), Ricc(2,4)],
[Ricc(3,0), Ricc(3,1), Ricc(3,2), Ricc(3,3), Ricc(3,4)],
[Ricc(4,0), Ricc(4,1), Ricc(4,2), Ricc(4,3), Ricc(4,4)] ])

R = sp.simplify((ginv*RiccMat).trace())

GMat = sp.simplify(RiccMat - sp.Rational(1,2)*g*R)
GMat

shows the Einstein tensor is

$\displaystyle{G_{ij}=\left(\begin{array}{ccccc} -\frac{3f(rf'+2f-2)}{2r^2} & 0 & 0 & 0 & 0 \\ 0 & \frac{3(rf'+2f-2)}{2r^2f} & 0 & 0 &0 \\ 0 & 0 & \frac{r^2f''+4rf'+2f-2}{2} & 0 & 0 \\ 0 & 0 & 0 & \frac{r^2f''+4rf'+2f-2}{2}\sin^2\theta & 0 \\ 0 & 0 & 0 & 0 & \frac{r^2f''+4rf'+2f-2}{2}\sin^2\theta\sin^2\varphi \end{array} \right)}$

The Einstein equation $G_{ij}=0$ imposes

$\qquad\displaystyle{rf'+2f-2=0\\ r^2f''+4rf'+2f-2=0}$

on the function $f$ . This case as well it reduces to a single equation

$\qquad\displaystyle{\frac{f'}{1-f}=\frac{2}{r}}$

which is solved by

$\qquad\displaystyle{f(r)=1-\frac{a}{r^2}}$

where $a$ is a constant of integration. The metric specified by this $f(r)$

$\displaystyle{g_{ij}=\left(\begin{array}{ccccc} -\left(1-\frac{a}{r^2}\right) & 0 & 0 & 0 & 0 \\ 0 & \left(1-\frac{a}{r^2}\right)^{-1} & 0 & 0 & 0 \\ 0 & 0 & r^2 & 0 & 0 \\ 0 & 0 & 0 & r^2\sin^2\theta &0 \\ 0 & 0 & 0 &0 & r^2\sin^2\theta\sin^2\varphi \end{array} \right)}$

is called 5-dimensional Schwarzschild metric.

Schwarzschild Sol. in 6D

It shall be straightforward to extend the previous analysis to the case of 6-dimensional spacetime $\xi^i=(t,r,\theta,\phi,\psi,\chi)$ with

$\displaystyle{g_{ij}=\left(\begin{array}{cccccc} -f(r) & 0 & 0 & 0 & 0 & 0 \\ 0 & \frac{1}{f(r)} & 0 & 0 & 0 & 0 \\ 0 & 0 & r^2 & 0 & 0 & 0 \\ 0 & 0 & 0 & r^2\sin^2\theta &0 & 0 \\ 0 & 0 & 0 &0 & r^2\sin^2\theta\sin^2\varphi & 0 \\ 0 & 0 & 0 &0 & 0 & r^2\sin^2\theta\sin^2\varphi\sin^2\psi \end{array} \right)}$