Symbolic and Numerical Analysis in General Relativity with Open Source Computer Algebra Systems

TL;DR

This paper compares open source computer algebra systems for tensor calculations in general relativity, demonstrating their capabilities through solving and visualizing the Klein-Gordon equation and geodesic motion.

Contribution

It provides a benchmark and detailed analysis of SageMath, Maxima, and Python-based GraviPy for tensor and algebraic computations in general relativity, highlighting their practical applications.

Findings

SageMath and SageManifolds effectively derive and visualize solutions.

Numerical Klein-Gordon solutions agree with analytical asymptotics.

Benchmark results guide future choice of tools for relativity calculations.

Abstract

We study three computer algebra systems, namely SageMath (with SageManifolds package), Maxima (with ctensor package) and Python language (with GraviPy module), which allow tensor manipulation for general relativity calculations along with general algebraic calculations. We present a benchmark of these systems using simple examples. After the general analysis, we focus on the SageMath and SageManifolds system to derive, analyze and visualize the solutions of the massless Klein-Gordon equation and geodesic motion with Hamilton-Jacobi formalism. We compare our numerical result of the Klein-Gordon equation with the asymptotic form of the analytical solution to see that they agree.

Figures

Tables

Table 1. Table 1: Benchmark results. (All numbers are in seconds).

Metric - Calculation	Python	SageMath	Maxima
Schwarzschild - Christoffel	0.000047	0.0000030	0.10
Schwarzschild - Ricci	0.000026	0.0000024	0.09
Schwarzschild - Einstein	0.000076	0.174	0.09
Kerr - Christoffel	0.000049	0.0000028	0.20
Kerr - Ricci	0.000028	0.0000024	1.67
Kerr - Einstein	0.000080	0.286	2.53

Equations

$$d{s}^2=g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }.$$

$$g_{\mu\nu}=\left[\begin{array}[]{cccc}\Big{(}1-\frac{2M}{r}\Big{)}&0&0&0\\ 0&-\Big{(}1-\frac{2M}{r}\Big{)}^{-1}&0&0\\ 0&0&-r^{2}&0\\ 0&0&0&-r^{2}sin^{2}\theta\end{array}\right].$$

$$g_{\mu\nu}=\left[\begin{array}[]{cccc}\Big{(}1-\frac{2M}{\rho^{2}}\Big{)}&0&0&\frac{2aMrsin^{2}\theta}{\rho^{2}}\\ 0&-\frac{\rho^{2}}{\Delta}&0&0\\ 0&0&-\rho^{2}&0\\ \frac{2aMrsin^{2}\theta}{\rho^{2}}&0&0&-sin^{2}\theta\big{[}(r^{2}+a^{2})+\frac{2a^{2}Mrsin^{2}\theta}{\rho^{2}}\big{]}\end{array}\right].$$

$${\rho }^2$$

$$\Delta$$

$$\frac{1}{\sqrt{-g}}{\partial }_{\mu }(\sqrt{-g}{g}^{\mu \nu }{\partial }_{\nu }\Phi )=0,$$

$$\Phi (t,r,\theta ,\varphi )=e^{-i\omega t}{e}^{ik\varphi }R(r)S(\theta ),$$

$$\text{KG}=$$

$$+$$

$$-$$

$$\text{KGradialpart}=$$

$$+\frac{\left(4M^{2}r-4M{r}^2+r^3\right)\frac{{\partial }^2}{(\partial r{)}^2}R\left(r\right)}{2M-r},$$

$$\text{KGangularpart}=-\frac{\cos \left(\theta \right)\frac{\partial }{\partial \theta }S\left(\theta \right)}{\sin \left(\theta \right)}-\frac{\left({\lambda }_{aux}\sin {\left(\theta \right)}^2-k^2\right)S\left(\theta \right)}{\sin {\left(\theta \right)}^2}-\frac{{\partial }^2}{(\partial \theta {)}^2}S\left(\theta \right).$$

$$\frac{dR}{dr}=R_{aux},$$

$$\frac{d{R}_{aux}}{dr}=-\frac{4M^2{R}_{\mathrm{aux}}+2r^2{R}_{\mathrm{aux}}+2\left({\lambda }_{aux}R-3r{R}_{\mathrm{aux}}\right)M+\left({\omega }^2{r}^3-{\lambda }_{aux}r\right)R}{4M^{2}r-4M{r}^2+r^3}.$$

$$R_{\ell}=\frac{C_{\ell}}{r}\sin\bigg{[}\omega r+2M\omega\ln(r)-\frac{\ell\pi}{2}+\arg\big{(}\Gamma(\ell+1-2iM\omega)\big{)}\bigg{]},$$

$$\frac{\partial S}{\partial \eta }-\frac{1}{2}{g}^{\mu \nu }\frac{\partial S}{\partial x^{\mu }}\frac{\partial S}{\partial x^{\nu }}=0,$$

$$\frac{d{x}^{\mu }}{d\eta }+g^{\mu \nu }\frac{\partial S}{\partial x^{\nu }}=0.$$

$$S(\eta ,t,r,\theta ,\varphi )=\frac{m^2}{2}\eta -Et+L\varphi +F(r)+G(\theta ),$$

$${\tt HJst}=\frac{m^{2}}{2}+\frac{E^{2}r}{2(2M-r)}+\frac{L^{2}}{2r^{2}}-\frac{2M-r}{2r}\bigg{(}\frac{dF(r)}{dr}\bigg{)}^{2},$$

$$\frac{dt}{d\eta }$$

$$\frac{dr}{d\eta }$$

$$\frac{d\varphi }{d\eta }$$

Full text

1 Introduction

Computer algebra systems are essential tools for theoretical physics for some decades. They are mainly important in general relativity (GR) where lengthy tensorial and differential geometry calculations are inevitable. Many of these programs have special internal or external packages for tensor manipulation and differential geometry calculations. For comprehensive reviews, see [1, 2, 3]. The review of early works can be found in [4]. These packages can also be extended for more specialized calculations [5].

Even though commercial programs dominate the area, codes written on commercial programs cannot be distributed easily since those programs may not be available for scientists with a lack of resources for purchasing such software. Besides, many of the manipulations treated in calculations do not need a sophisticated computation engine.

Open source computer algebra systems such as SageMath (also known as “Sage”) [6] and Maxima [7] provide a complete toolkit for general relativity and quantum field theory applications with their particular packages. Some freely available programming languages such as Python [8] also offer special tools for the aforementioned manipulations.

We choose to study general systems rather than specialized tensor manipulation or general relativity packages (e.g. Cadabra [9], Redberry [10], etc.). The reason for this choice is that tensor manipulation is generally an auxiliary step in GR calculations. A general program which can deal with the symbolic and numerical analysis of the results and produce graphical outputs constitute a complete calculation toolkit. In most cases, computer algebra systems (or programming languages) are supported with particular packages for tensor manipulation. We also exclude specialized tools like GYOTO [11] in our analysis.

In this work, we will first employ SageMath (with SageManifolds package [12]), Maxima (with ctensor package [13]), and Python language (with Sympy and GraviPy modules [14]) for some essential calculations in general relativity and present benchmark results for these systems.

The ever-developing open source SageMath program has gathered many utilities such as Maxima, GAP, R, and the power of Python language with well-known Python modules like NumPy, SymPy and matplotlib. SageMath can be installed on personal computers and moreover it has a powerful cloud computing server on which the user can work on projects anywhere and share them with other users easily [6]. The package SageManifolds for tensor and differential geometry calculations is included in SageMath and it does not require an additional installation process. These properties make the Sage+SageManifolds system one of the best open source choices for general relativity and quantum field theory. A comprehensive lecture note on SageManifolds is available in [15].

Analysis of the Klein–Gordon equation, mainly its radial part is the first step in most quantum gravity problems involving black holes [16]. We will study its solution using SageMath and see that the numerical result found by SageMath agrees with the asymptotic analytical result.

Computationally, it is easier to define and solve the first order differential equations rather than the second order equations. Hamilton–Jacobi formalism yields first order differential equations for the geodesic motion [17]. We will use this formalism to derive the equations and solve them as a numerical initial value problem with SageMath routines. Numerical results need to be presented in a graphical way to see their structure efficiently and SageMath system is equipped with many visualization tools. The “Examples” section in SageManifolds web page [12] also presents some applications on geodesics. However, our code provides a step-by-step procedure for the Hamilton-Jacobi approach to geodesic motion with explicit numerical analysis.

In the next section, we define the spacetimes to be used. In the third section we review three computation systems and give simple examples. We also give a benchmark of these systems as a subsection. In the fourth section we focus on the SageMath+SageManifolds system to analyze the massless Klein–Gordon equation and geodesic motion.

All codes studied in this paper can be accessed from the GitHub address [18], organized by branches.

2 Schwarzschild and Kerr solutions

We will be working for Schwarzschild and Kerr spacetimes with the metric signature $(+,-,-,-)$. In general, a line element is defined by

[TABLE]

If one takes $G_{4}=c=1$ (where $G_{4}$ is Newton’s gravitational constant in four dimensions and $c$ is the speed of light in vacuum), the Schwarzschild metric can be written as [17]

[TABLE]

We will take the order of coordinates as $\{t,r,\theta,\phi\}$ and $M$ is the mass of the black hole. The coefficient of the radial part is singular at $r=2M$ which describes the “event horizon”.

In the Boyer-Lindquist coordinates, the Kerr black hole has the metric [17]

[TABLE]

Here, $M$ is the mass and $a=J/M$ where $J$ is the angular momentum. Two functions are defined as

[TABLE]

The roots of $\Delta=0$ correspond to the locations of the Cauchy horizon (the smaller root) and the event horizon (the larger root) of the black hole.

3 Open source computer algebra systems for basic GR calculations

We will only focus on the procedures of

• •

Defining a spacetime

• •

Calculating Christoffel symbols ($\Gamma^{\mu}_{\nu\rho}$)

• •

Calculating Ricci tensor ($R_{\mu\nu}$)

• •

Calculating Einstein tensor ($G_{\mu\nu}$)

• •

Displaying components

in our examples. The examples will be studied using three computation systems, namely

• •

SageMath+SageManifolds

• •

Maxima+ctensor

• •

Python+GraviPy

The codes are enhanced with comment lines which can be helpful for analyzing the procedures step by step. Details on the procedures and further examples can be found in related references.

3.1 SageMath+SageManifolds

SageMath is an open source computer algebra system which collects many powerful open source packages and modules with a Python-like language [6]. SageManifolds was first started as an independent tensor analysis and differential geometry package to be installed on SageMath [12]. Now, SageManifolds is an internal package for SageMath and does not require an additional installation.

We start our program by reseting the SageMath environment. Then we define the four-dimensional manifold, the black hole mass and the coordinates. The coordinates are defined along with their ranges. We define the Lorentzian metric by giving its components as a matrix. Finally, we calculate and display Christoffel symbols, Ricci tensor and Einstein tensor elements.

⬇

1reset()

2# Define 4-dim. the manifold "Man":

3Man = Manifold(4, 'Man', r'\mathcal{M}')

4

5# Define the parameter "M" (mass):

6M = var('M')

7

8# Define the coordinates {t=0, r=1, theta(=th)=2, phi=3} with ranges

9# (BL = Boyer-Lidquist)

10BL.<t,r,th,ph> = Man.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')

11

12# Define the metric "g" on manifold "Man":

13g = Man.lorentzian_metric('g')

14

15# Enter the metric components:

16g[0,0] = (1-(2*M)/r)

17g[1,1] = -1/(1-(2*M)/r)

18g[2,2] = -r^2

19g[3,3] = -(r*sin(th))^2

20

21# Display the metric

22show('The metric:')

23show(g.display())

24#####################

25# Christoffel symbols

26nab = g.connection()

27# Display all components

28show(nab.display())

29# Display a single component

30show(nab[3,2,3])

31#####################

32# Ricci tensor

33Ric = g.ricci()

34# Display all components

35show(Ric.display())

36# Display a single component

37show(Ric[1,1])

38#####################

39# Einstein tensor

40ET=g.ricci()-(1/2)gg.ricci_scalar()

41# Display all components

42ET.set_name('E')

43show(ET.display())

44# Display a single component

45show(ET[1,2])

The lines between 16-19 defines the Schwarzschild metric. In order to define the Kerr metric, we need to change this part with

⬇

16a = var('a')

17rhosq=r^2+(a^2)*cos(th)^2

18Delta=r^2-2Mr+a^2

19

20g[0,0] = (1-(2Mr)/rhosq)

21g[1,1] = -rhosq/Delta

22g[2,2] = -rhosq

23g[3,3] = -(sin(th)^2)((r^2+a^2)+(2(a^2)Mr*sin(th)^2)/rhosq)

24g[0,3] = (2aMrsin(th)^2)/rhosq

It should be noted that, we defined the functions $\rho^{2}$ and $\Delta$, and the rotation parameter $a$. The Kerr black hole has non-diagonal, symmetric components unlike the Schwarzschild metric which is diagonal.

3.2 Maxima+ctensor

Maxima computer algebra system is a freely available descendant of Macsyma which is known as the first general, multipurpose computer algebra system which inspired many other systems for years [7]. V. Toth reviews tensor packages in Maxima in his article [13]. Here, we will use the ctensor package which provides component tensor manipulation.

We first kill all environmental variables, then load the ctensor package. We set ratwtlvl:false for no truncation and ratfac:true to factorize the tensor components automatically. We then define the dimension of the spacetime and the coordinates. We enter the metric as a matrix. Before the calculation of Christoffel symbols, Ricci tensor and Einstein tensor elements, we find the inverse tensor with invert.

⬇

1kill(all);

2if get('ctensor,'version)=false then load(ctensor);

3(ratwtlvl:false,ratfac:true);

4("Specify the dimension of the manifold and the coordinate labels.")$

5(dim:4,ct_coords:[t,r,theta,phi]);

6("Enter the metric.")$

7lg:matrix([(1-2M/r),0,0,0],[0,-1/(1-2M/r),0,0],[0,0,-r^2,0],[0,0,0,-(r^2)*sin(theta)^2]);

8ug:invert(lg)$

9("Compute the Christoffel symbols and display all components")$

10christof(mcs)$

11("Compute the Ricci tensor and display all components")$

12uricci(true)

13("Compute the Einstein tensor and display all components")$

14einstein(true);

15("Display single components")$

16mcs[3,4,4];

17ric[1,2];

18ein[2,2];

The lines 6-7 should be changed as

⬇

6rhosq:r^2+(a^2)*cos(theta)^2$

7Delta:r^2-2Mr+a^2$

8("Enter the general static spherically symmetric metric.")$

9lg:matrix([(1-(2Mr)/rhosq),0,0,(2aMrsin(theta)^2)/rhosq],[0,-rhosq/Delta,0,0],[0,0,-rhosq,0],[(2aMrsin(theta)^2)/rhosq,0,0,-(sin(theta)^2)((r^2+a^2)+(2(a^2)Mr*sin(theta)^2)/rhosq)]);

for the Kerr metric.

3.3 Python+GraviPy

Python is a multipurpose, object-oriented programming language [8] which can easily be expanded with modules. The module GraviPy provides tensor calculation methods and it works on a freely available symbolic analysis module SymPy [14].

We start our program by importing the GraviPy module. Then we define the four–vector of coordinates and the black hole mass. We define the metric as a diagonal matrix. We then calculate and display Christoffel symbols, Ricci tensor and Einstein tensor elements.

⬇

1from gravipy import *

2#####################

3# Coordinates (\chi is the four-vector of coordinates)

4t, r, theta, phi, M = symbols('t, r, theta, phi, M')

5x = Coordinates('\chi', [t, r, theta, phi])

6#####################

7# Metric tensor

8Metric = diag((1-2M/r), -1/(1-2M/r), -r2, -r2*sin(theta)**2)

9g = MetricTensor('g', x, Metric)

10#####################

11# Christoffel symbols

12Ga = Christoffel('Ga', g)

13# Display all components

14print(Ga(All, All, All))

15# Display a single component

16print(Ga(1,2,1))

17#####################

18# Ricci tensor

19Ri = Ricci('Ri', g)

20# Display all components

21print(Ri(All, All))

22# Display a single component

23print(Ri(1, 2))

24#####################

25# Einstein tensor

26G = Einstein('G', Ri)

27# Display all components

28print(G(All, All))

29# Display a single component

30print(G(3, 3))

31#####################

The lines 2-8 should be modified as

⬇

2#####################

3# Coordinates (\chi is the four-vector of coordinates)

4t, r, theta, phi, M, a, rhosq, Delta = symbols('t, r, theta, phi, M, a, rhosq, Delta')

5x = Coordinates('\chi', [t, r, theta, phi])

6#####################

7# Metric tensor

8rhosq = r2+(a2)*cos(theta)**2

9Delta = r2-2Mr+a2

10Metric = Matrix([[(1-(2Mr)/rhosq),0,0,(2aMrsin(theta)2)/rhosq],[0,-rhosq/Delta,0,0],[0,0,-rhosq,0],[(2aMrsin(theta)2)/rhosq,0,0,-(sin(theta)2)*((r2+a2)+(2(a**2)Mrsin(theta)**2)/rhosq)]])

in order to define the Kerr metric. The diagonal matrix definition of the Schwarzschild case is changed with a general metric.

3.4 Benchmark for the open source computer algebra systems

We performed some calculations for the Schwarzschild and Kerr spacetimes on Python+GraviPy (Python 2.7.15), SageMath+SageManifolds (SageMath 8.3) and Maxima+ctensor (Maxima 18.02.0) systems.

The Schwarzschild metric has only diagonal elements, while Kerr solution has non–diagonal elements in its metric. We aimed to measure the effect of this difference in the computations. We used the metric given in equation 2 for the Schwarzschild case and for the Kerr metric the metric is given in equation 3.

For our analysis, we used the code–block profiling method and measured the wall–clock time as it has more importance for the general user. We calculated the Christoffel symbols, Ricci tensor and Einstein tensor one by one on each system. Metric definition and displaying the results are not included in the time measurement. By averaging the timing results for 10 runs on each system, the Table 1 is generated.

In Python, we used the datetime module. This method shows the time value up to six decimals.

⬇

1from gravipy import *

2from datetime import datetime

3...

4riccitime_=datetime.now()

5Ri=Ricci('Ri',g)

6_riccitime=datetime.now()

7print "Elapsed time...", riccitime-riccitime

8...

In SageMath, we employed the inline timeit command.

⬇

1...

2timeit(g.ricci())

3...

In Maxima, we used elapsed_real_time() command which also shows time with two decimals.

⬇

1...

2t0:elapsed_real_time()$

3uricci(false);

4t1:elapsed_real_time()-t0;

5...

The test computer has an Intel(R) Core(TM) i7-4930K CPU @ 3.40GHz, 16 GB (DIMM DDR3 Synchronous 1066 MHz) RAM, NVIDIA GeForce GTX 650 Ti Boost GPU, 120GB SSD and 2 TB HDD with Ubuntu 18.04.

According to the benchmark results, Python seems to be the best choice for these calculations. However, the calculation commands lack the simplification routines that are needed for a useful result. For example, displaying even one component of the Einstein tensor of the Kerr black hole is impossible in an acceptable duration in Python.

Focusing on the Kerr solution, we see that SageMath would be the best option with its speed and yielding convenient results for further calculations and manageable outputs.

4 Scalar wave equation and geodesics with SageMath and SageManifolds

SageMath, being powered by SageManifolds, provides an easy-to-use and combined toolkit for the general user. We will focus on the SageMath+SageManifolds system to perform a set of calculations for the Schwarzschild metric.

We will first define the spacetime by declaring its variables and components. Then we will perform calculations for two examples. In the first one, we derive the Klein–Gordon equation for a massless scalar field, extract its radial part and solve this differential equation numerically. We will compare our numerical result graphically with the asymptotic form of the analytical result. In the next example, we will perform a very simple geodesic analysis for this spacetime using the Hamilton–Jacobi formalism.

Generalization of the example codes to other metrics is straightforward. The codes for the calculation of Klein–Gordon and Hamilton–Jacobi equations can be used for general metrics without modification. For the detailed calculations, we applied some metric-related information in the code to see the manipulation, equation solving and plotting capabilities of the computation system.

The reader should follow the line numbers to execute the codes without problem. The code for the definition of the spacetime is between lines (1-24). The analysis of the Klein–Gordon equation starts in line 5 and ends in line 78. Then in the next section, the study of geodesics starts again with line 1 and ends in line 38. This means that, the reader should first execute the lines (1-24) for the metric definition before each physical example.

4.1 Definition of the spacetime in SageManifolds

This part was studied before in Section 3.1, while giving examples for simple calculations in GR. Nevertheless, we will place it here in order to present a complete code structure.

⬇

1reset()

2

3# Define 4-dim. the manifold "Man":

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

5

6# Define the parameter "M" (mass):

7M = var('M')

8

9# Define the coordinates {t=0, r=1, theta(=th)=2, phi=3} with ranges

10# (BL = Boyer-Lidquist)

11BL.<t,r,th,ph> = Man.chart(r't r:(0,+oo) th:(0,pi):\theta ph:(0,2*pi):\phi')

12

13# Define the metric "g" on manifold "Man":

14g = Man.lorentzian_metric('g')

15

16# Enter the Schwarzschild metric components:

17g[0,0] = (1-(2*M)/r)

18g[1,1] = -1/(1-(2*M)/r)

19g[2,2] = -r^2

20g[3,3] = -(r*sin(th))^2

21

22# Display the metric

23show('The Schwarzschild metric:')

24show(g.display())

In the Schwarzschild case, the manifold has four dimensions and the only variable is the black hole mass $M$. We define four spacetime coordinates and enter the components of the diagonal metric.

4.2 Klein–Gordon equation in curved spacetime

Klein–Gordon equation for a massless scalar field can be written as [19]

[TABLE]

where $\Phi$ is the scalar field and it can be decomposed with the Ansatz

[TABLE]

for the Schwarzschild metric.

We will first define the variables $\omega$, $k$ and the result KG. The inverse metric and $\sqrt{-g}$ are present in the equation, thus are calculated. Then the full scalar function $\Phi$, radial and angular parts of the solution Ansatz are defined. After giving the Ansatz, we start calculating the Klein–Gordon equation in two for loops. The outer loop is over $\mu$ components and the inner one is over $\nu$ components.

⬇

5# Defining variables:

6var('omega,k,KG')

7

8# Inverse metric:

9ginv = g.inverse()

10

11# The square root of the absolute value

12# of the metric determinant:

13sqrtabsdetg=g.sqrt_abs_det().expr()

14

15# The scalar function Phi(t,r,th,phi):

16# The dependence on all coordinates

17# is provided by "(*BL)"

18Phi=function('Phi')(*BL)

19

20# The scalar field Ansatz is given here.

21# R : Radial part,

22# S : Angular part

23R=function('R')(r)

24S=function('S')(th)

25Phi=exp(-Iomegat)exp(Ik*ph)RS

26

27# Calculating the Klein-Gordon equation:

28KG=0

29for mu in range(len(BL[:])):

30 for nu in range(len(BL[:])):

31 KG=KG+diff((ginv[mu,nu].expr()sqrtabsdetgdiff(Phi,BL[nu])),BL[mu])

32

33# Displaying the Klein-Gordon equation "KG"

34show('The full Klein-Gordon equation (variable name is KG):')

35show(KG)

This part of the code is general. The reader should only change the solution Ansatz accordingly and the code can find the Klein-Gordon equation for any spacetime with any number of dimensions. The result is stored in variable KG.

In our example, the Klein–Gordon equation is found as

[TABLE]

We can find the radial and angular parts of the Klein-Gordon equation after analyzing its structure. This part is semi-automatic and the user should supply some information.

We first divide the equation by a convenient factor and extract its operands in a vector. The factor is generally the solution Ansatz multiplied by some angular or radial terms. Then we define the separation constant $\lambda_{aux}$. Using a loop over all operands, we decide which one belongs to the radial part and which one belongs to the angular part of the equation. We take the radial derivative of the components and if the result is zero, the component is angular, if it is not zero, we add it to the radial part. We use collection and simplification commands to see the results in a convenient shape.

⬇

24# We can analyze the Klein-Gordon equation

25# to see how it is separated into

26# radial and angular parts.

27

28# Common factors will be divided

29# This part should be given by the user

30divideKGby=exp(-Iomegat)exp(Ik*ph)*sin(th)RS

31finalKG=expand(KG/divideKGby)

32

33# Extract the operands in the expression:

34fkgops=finalKG.operands()

35

36# Find radial and angular parts:

37# lambda_aux is the separation constant

38var('lambda_aux')

39KGradialpart=lambda_aux

40KGangularpart=-lambda_aux

41for term in fkgops:

42 if diff(term,r)==0:

43 KGangularpart=KGangularpart+term

44 else:

45 KGradialpart=KGradialpart+term

46

47KGradialpart=expand(KGradialpart*R).simplify_full().collect(R).collect(diff(R,r)).collect(diff(R,r,r))

48KGangularpart=expand(KGangularpart*S).simplify_full().collect(S).collect(diff(S,th)).collect(diff(S,th,th))

49

50show('The radial part (variable name is KGradialpart):')

51show(KGradialpart)

52show('The angular part (variable name is KGangularpart):')

53show(KGangularpart)

The radial part is then found as

[TABLE]

and the angular part is

[TABLE]

We will concentrate on the radial part which is the first step in most quantum gravity problems involving black holes [16]. Further simplifications are obvious but SageMath’s desolve command can not give a symbolic solution for the radial equation. However, numerical solutions can be studied using related methods. We will use desolve_system_rk4 as an example and solve the radial equation numerically. This method uses a fourth–order Runge–Kutta scheme and in fact, the command desolve_system_rk4 is used as a wrapper for the Maxima command rk [6].

Numerical solvers can generally deal with first order equations. Our second order equation will yield two first order differential equations which will be solved simultaneously. We will define the first derivative of the radial function as an auxiliary function

[TABLE]

and place it in the equation 9 to have

[TABLE]

The following code finds this set of equations as radial1 and radial2. The solver can deal with the right hand sides of $\frac{df(x)}{dx}=g(f(x),x)$ type equations. Thus we isolate the derivatives by solving the differential equations as algebraic equations. radial1 and radial2 are the right hand sides of our equations to be solved.

⬇

35# The auxiliary function

36R_aux=function('R_aux')(r)

37

38# Two auxiliary equations

39radial1aux=diff(R,r)-R_aux

40radial2aux=KGradialpart.subs(diff(R,r)==R_aux).subs(diff(R,r,r)==diff(R_aux,r))

41

42# Right hand sides of the derivatives

43radial1=(solve(radial1aux==0,diff(R,r)))[0].right()

44radial2=(solve(radial2aux==0,diff(R_aux,r)))[0].right()

45

46# Getting the outputs in input format

47# (Outputs will be copied)

48print(radial1)

49print(radial2)

In their present forms, $R$ and $R_{aux}$ are defined as functions. However, in order to use desolve_system_rk4 we need to define unknowns as variables. To do this, we need to copy the outputs for radial1 and radial2 and change them accordingly. We display the outputs by print command to see them in the input form which enables us to copy them easily.

After a formal analysis of the radial part, one can see that the radial solution can be given in terms of confluent Heun functions [20]. General and confluent forms of the Heun function are encountered in many applications in physics, especially as solutions of the wave equations [21, 22, 23, 24, 25, 26].

Numerical computation of the general and confluent Heun functions are adapted by Oleg V. Motygin for GNU Octave [27, 28]. However, no freely available packages or modules can give closed symbolic solutions of these equations. Therefore we cannot compare our result with the full analytical solution. The asympotic form of the confluent Heun function is given in reference [20] as

[TABLE]

where $\lambda_{aux}=\ell(\ell+1)$.

We first import the differential equation solver from sage.calculus.desolvers. We then define the variables (including the unknown functions) and the equations and copy the equations from the outputs of the code above. We give some numerical values to the parameters arbitrarily and solve the system for some arbitrary set of initial conditions. In our example we have $R(r=0.3)=1$ and $\frac{dR}{dr}=R_{aux}(r=0.3)=0.5$. The solution is stored in radsol.

radsol has the structure $[r,R,R_{aux}]$. points=[[i,j] for i,j,k in radsol] command creates [i,j] (or $[r,R]$) points for plotting. radialsolution stores the plot of the numerical result.

In the next part, we define the asymptotic form of the analytical solution and plot it for the same parameter set (we take $C_{\ell}=250$ to match the amplitude). We display both plots together to show the agreement.

⬇

53# Import the solver

54from sage.calculus.desolvers import desolve_system_rk4

55

56# Define unknowns as variables

57var('R,R_aux,r')

58

59# Define equations

60radial1=R_aux

61radial2=-(4M^2R_aux + 2r^2R_aux + 2*(lambda_auxR - 3rR_aux)M + (omega^2r^3 - lambda_auxr)R)/(4M^2r - 4M*r^2 + r^3)

62

63# Substitute numerical values for parameters

64radial2=radial2.subs(M=0.1,omega=0.2,k=2.0,lambda_aux=2.0)

65

66# Solve the system and plot the solution

67radsol=desolve_system_rk4([radial1,radial2],[R,R_aux],ics=[0.3,1,0.5],ivar=r,end_points=100,step=0.01)

68points=[[i,j] for i,j,k in radsol]

69radialsolution=list_plot(points,axes_labels=['$r$','$R$'],legend_label='Numerical solution')

70

71# Asymptotic form of the analytic solution

72var('Cl,L')

73Rasym=Cl*(1/r)sin(omegar+2Momegalog(r)-(Lpi/2)+arg(gamma(L+1-2IM*omega)))

74Rasymnum=Rasym.subs(Cl=250,M=0.1,omega=0.2,L=1.0)

75asympplot=plot(Rasymnum,(r, 20, 100),linestyle='',marker='x',color='red',legend_label='Asymptotic solution')

76

77# Display both plots

78show(asympplot+radialsolution)

The plot 1 of the solution gives an idea on the behavior of the radial function. Numerical solution and the asymptotic form of the analytical solution are plotted together to show that they agree for large $r$.

4.3 Visualizing geodesics

The geodesic motion in a spacetime can be described by the Hamilton–Jacobi equation [17].

[TABLE]

where $S$ denotes the Hamilton’s principal function, $g^{\mu\nu}$ is the inverse metric and $\eta$ is an affine parameter. The orbital equations then found using

[TABLE]

This formalism gives first order differential equations which are needed in numerical schemes as mentioned in the previous section.

The Hamilton’s principal function is decomposed as

[TABLE]

for the Schwarzschild spacetime. For the equatorial geodesics we take $\theta=\pi/2$ and thus $G(\theta)=0$. The metric components also change according to this constraint.

Before beginning our analysis, it should be emphasized that decomposition of the Hamilton’s principle function $S$ (and the scalar field $\Phi$ in the previous section) is not evident and needs a detailed analysis in general spacetimes, which is beyond the scope of our work. Moreover, analysis of geodesics is a very detailed study and our primitive example here aims only to visualize some geodesics.

The user should execute the lines (1-24) for the definition of the spacetime before running the codes below.

We start by defining our variables and functions. We also calculate the inverse metric as needed in the equations. After giving the Hamilton’s principle function Ansatz, we calculate the Hamilton–Jacobi equation in two loops and set it to the variable HJfull. We call it HJfull as we have not put any conditions ($\theta=\pi/2$, etc.) on the equation yet.

The right hand sides (as needed by differential equation solvers) of the geodesic (orbital) equations are then calculated (using nested loops) and put in the vector geodeqnrhs. We display the equations both in the LaTeX format and as a vector.

⬇

1# Define variables, functions and calculate inverse metric

2var('eta,m,E,L,S,HJfull')

3F=function('F')(r)

4G=function('G')(th)

5ginv = g.inverse()

6

7# Define the principal function Ansatz

8S=((etam^2)/2)-Et+L*ph+F+G

9

10# Calculate the Hamilton-Jacobi equation

11HJfull=0

12for i in range(len(BL[:])):

13 for j in range(len(BL[:])):

14 HJfull=HJfull+ginv[i,j].expr()*diff(S,BL[i])*diff(S,BL[j])

15HJfull=(diff(S,eta)-(1/2)*HJfull)

16show('The Full Hamilton-Jacobi equation (variable name is HJfull):')

17show(HJfull)

18

19show('The geodesic equations in LaTeX form (variable name is geodeqnrhs):')

20geodeqnrhs=zero_vector(SR, len(BL[:]))

21

22for mu in range(len(BL[:])):

23 for nu in range(len(BL[:])):

24 geodeqnrhs[mu]=geodeqnrhs[mu]-(ginv[mu,nu].expr())*diff(S,BL[nu])

25 writeresult='D0($\eta$) = $$' %(BL[mu],latex(geodeqnrhs[mu]))

26 show(writeresult)

27

28show('Right hand sides of the geodesic equations as a vector')

29show(geodeqnrhs)

The code above is general and it can work for any spacetime if the Hamilton’s principal function is given accordingly. The codes below depend on the Schwarzschild metric.

Conventionally, we take $\theta=\pi/2$ to find the equatorial geodesics. We call the Hamilton–Jacobi equation with this constraint as HJst and the right hand sides of the orbital equations as geodeqnrhsst. We substitute $\theta=\pi/2$ and $G(\theta)=0$ in the equations found in the general code.

⬇

1# Take theta = pi/2:

2var('HJst')

3

4# Hamilton-Jacobi equation for theta=pi/2 (variable name is HJst)

5HJst=(HJfull.subs(diff(G)==0)).subs(th=pi/2)

6show(HJst)

7

8# Right hand sides of the geodesic equations for theta=pi/2

9# (variable name is geodeqnrhsst)

10geodeqnrhsst=(geodeqnrhs.subs(diff(G)==0)).subs(th=pi/2)

11show(geodeqnrhsst)

We find

[TABLE]

and the nonzero components of geodeqnrhsst are

[TABLE]

We are now ready to solve the orbital equations. Here, we will use another solver, desolve_odeint by importing from sage.calculus.desolvers. This solver uses scipy.integrate module of Python. We will plot the geodesic curve, and a black disc with a radius equal to the event horizon radius in order to indicate the black hole. Thus, we import Circle from sage.plot.circle. We define the variables and give them arbitrary numerical values. We isolate $\frac{dF(r)}{dr}$ as derofradfun from HJst and place it in $\frac{dr}{d\eta}$ equation.

We would like to solve time-dependent $\frac{dr}{dt}$ and $\frac{d\phi}{dt}$ equations, instead of the equations with affine-parameter ($\eta$) dependence. Thus, we divide $\frac{dr}{d\eta}$ and $\frac{d\phi}{d\eta}$ by $\frac{dt}{d\eta}$, and set them as our equations: geodeqn1 and geodeqn2.

We solve our equations for arbitrarily set initial conditions and put the results in the variable sol. The $0^{th}$ column of sol carries the $r$ values and the $1^{st}$ column carries the $\phi$ values. Using these, we generate $\{x,y\}$ points, where $x=rcos(\phi)$ and $y=rsin(\phi)$. The line plot of these points forms the geodesic curve. We plot this curve and the black disc (a circle with parameters fill=True,rgbcolor=‘black’) together to display the behavior.

⬇

1# Importing the solver

2from sage.calculus.desolvers import desolve_odeint

3#Importing circle for visualizing the black hole

4from sage.plot.circle import Circle

5

6var('m_aux,L_aux,E_aux,M_aux,r_initial,ph_initial,eta_end,step_size')

7

8##########################

9# Variables

10m_aux=1 # Either 1 (timelike) or 0 (null).

11M_aux=1

12L_aux=4

13E_aux=1.0

14step_size=0.1

15eta_end=500

16r_initial=2.1*M_aux

17ph_initial=0.3

18##########################

19

20# Derivative of F (the radial function)

21# (variable name is derofradfun)

22# We will use the first root

23derofradfun=solve(HJst,diff(F,r))

24

25# Define equations to solve

26# dr/dt = geodeqn1

27# dphi/dt = geodeqn2

28geodeqn1=((geodeqnrhsst[1]/geodeqnrhsst[0]).subs(diff(F,r)==derofradfun[1].rhs())).subs(E=E_aux,L=L_aux,m=m_aux,M=M_aux)

29geodeqn2=(geodeqnrhsst[3]/geodeqnrhsst[0]).subs(E=E_aux,L=L_aux,m=m_aux,M=M_aux)

30

31# Solve the equations

32sol=desolve_odeint([geodeqn1,geodeqn2],[r_initial,ph_initial],srange(0,eta_end,step_size),[r,ph])

33p=line(zip(sol[:,0]*cos(sol[:,1]),sol[:,0]*sin(sol[:,1])))

34

35# Plot the black hole as a circle

36# Show the geodesics and the circle on the same plot

37C=circle((0,0),2*M_aux,fill=True,rgbcolor='black')

38show(C+p)

The figure 2 shows an example of a null geodesic ($m=m_{aux}=0$) for an arbitrary set of parameters. Figure 3 shows the stable circular orbit for a timelike particle ($m=m_{aux}=1$) and Figure 4 shows another example of a timelike geodesic.

5 Conclusion

We studied three computer algebra systems which allow tensor manipulation for GR calculations, and showed how one can define a spacetime, calculate Christoffel symbols, Ricci tensor and Einstein tensor for this spacetime using the specialized commands of these systems. We used SageMath with its tensor manipulation and differential geometry module SageManifolds, Maxima with its tensor component manipulation package ctensor and the Python language with the GraviPy module which runs on the symbolic analysis module SymPy. A benchmark for these systems is also provided.

In the next part, our main focus was the SageMath+SageManifolds system to perform two examples: Massless Klein–Gordon equation calculations and geodesic visualization for the Schwarzschild geometry. Primarily, the general codes that can run for general spacetimes are given, and then the analysis is specialized for the Schwarzschild case.

In our first example, we derived the massless Klein–Gordon equation as a partial differential equation using our code. We then separated it to the radial and angular parts and analyzed the radial part in detail. We confirmed our numerical result graphically using the asymptotic form of the analytical result. This part involves a direct numerical integration of a confluent Heun equation resulting from the radial part of the Klein–Gordon equation and its comparison with the asymptotic analytical result which is rare in the literature.

In the second example, the orbit of a null or timelike particle around a Schwarzschild black hole was constructed from the geodesic equations. There are multiple methods of doing that, and among them, the chosen one here is the Hamilton–Jacobi formalism as it yields first order differential equations that can be treated easier than the second order equations on the computer. The orbits are visualized using the plotting tools of SageMath. Our example is far from a complete analysis of geodesics although it provides a basis for such a study.

A computation system should be chosen according to the needs of the problem. For example, if a research topic depends on symbolic manipulation of special functions of mathematical physics, commercial packages are inevitable for most of the cases in the present situation. However, many problems do not need such specialized calculations and numerical analysis is sufficient to see the result. As seen in the examples, in many cases freely available packages are capable of forming a complete system for scientific problems.

Acknowledgement

We would like to thank Profs. Neşe Özdemir, Durmuş Ali Demir and Éric Gourgoulhon for stimulating discussions. We also thank our anonymous referee for the constructive comments which helped us to improve the manuscript. This work is partially supported by TÜBİTAK, the Scientific and Technological Council of Turkey.