DEV Community: Jyotirmaya Shivottam

GSoC 2020: Blog 5 - Adding Kerr Null Geodesics functionality to EinsteinPy

Jyotirmaya Shivottam — Sun, 30 Aug 2020 17:33:57 +0000

Null Geodesics functionality has been implemented into EinsteinPy, with PR#527, having been merged 🎉🎉. I apologize for no blogs in the past 3 weeks. There was a COVID situation here, that required multiple tests and isolation and all that it entails. This led to me foregoing an entire week. And, when that had settled, I had to take the call on abandoning the plan of numerically integrating the Geodesics equations, due to the massive error accumulation, as discussed in my last blog. A confusing fact about that, was that Mathematica could still keep the error build-up to a minimum, while Python simply could not, even with adaptive and symplectic schemes. But the symplectic schemes did bring the error down, by around 2 orders of magnitude, which gave me the idea to take a Hamiltonian approach, which would increase the number of ODEs to solve, but drop the order by 1. And, as it turns out, the Kerr Hamiltonian is separable (Carter, 1968a [1]), which makes the implementation even simpler. In this blog, I will be discussing this approach, which has finally led to proper geodesic calculations. I have also included some plots (and a cool animation) for Kerr & Schwarzschild Null-like (and Time-like) geodesics.

Some Physics...

In Chapter 33 of Gravitation [2], the authors expound on Carter's seminal paper from 1968, titled, "Global Structure of the Kerr Family of Gravitational Fields", and present some nice results from it, one of which is a derivation of the Kerr (super-)Hamiltonian, which can be written as follows (in the M-Unit system ( $G = c = M = 1$ ):

\mathcal{H} = -\frac{(a^4 (E^2 - 2 p_r^2) - 8 a E L r - 2 r (p_\theta^2 (-2 + r) + p_r^2 (-2 + r)^2 r - E^2 r^3) + a^2 (2 L^2 - 2 p_\theta^2 - 4 p_r^2 (-2 + r) r + E^2 r (2 + 3 r)) + (a^2 + (-2 + r) r) (a^2 E^2 \cos 2\theta - 2 L^2 \csc\theta^2)}{4 (a^2 + (-2 + r) r) (r^2 + a^2 \cos\theta^2))}

where,

E = -p_t

and

p_\phi

are the energy and orbital angular momentum of the test particle, respectively. Note that, this Hamiltonian is for a general test particle, i.e., it can be massive or massless. Then, the dynamical equations of motion can be derived easily, using Hamilton's principle, i.e.:

\frac{\mathrm{d}q_i}{\mathrm{d}\lambda} = \frac{\partial\mathcal{H}}{\partial p_i} \quad \frac{\mathrm{d}p_i}{\mathrm{d}\lambda} = -\frac{\partial\mathcal{H}}{\partial q_i}

I calculated these in Mathematica, and the corresponding notebooks and the Python code, making use of these, can be accessed here.

...and some plots

Unfortunately, even with these first order ODEs, the error accumulation issue in Python persisted, as can be observed in the plots below. Note that, these results were obtained with a symplectic leapfrog solver, which should, in principal keep the error build-up to a minimum.

Kerr Null-like Escape
Although, for shorter integration durations, the results were good.

Kerr Null-like Capture

After discussions with my mentors, I looked into other languages, that could help and we chose Julia, due to its excellent DifferentialEquations.jl suite and "closeness" with Python. Another key bit is that, the HamiltonianProblem type, offered by DiffEqPhysics, immensely simplifies the process of solving the system, as it uses Forward Mode Automatic Differentiation to automatically calculate the partial derivatives from the Hamiltonian. The separable nature of the Hamiltonian helps here. Considering all this, I implemented a module in Julia and voilà, the results are accurate, even for some quirky geodesics.

Kerr Null-like Capture (Plotted using `Plots.jl`)
Kerr Time-like Whirl (Plotted using `Plots.jl`)

Now came the problem of integrating the Julia code with EinsteinPy, for which I looked towards PyJulia. However, it has some issues with installation on *nix systems. So, I opted to write my own wrapper, using Python's subprocess. and, with the help of my GSoC mentor, Shreyas, packaged the Julia module and the Python wrapper into what is now einsteinpy_geodesics, an add-on module to EinsteinPy.

einsteinpy / einsteinpy-geodesics

Python wrapper for a Julia solver for geodesics in the Kerr family of spacetimes. Maintainer : @jes24

Name:	EinsteinPy Geodesics
Website:	https://docs.geodesics.einsteinpy.org/en/latest/
Version:	0.2.dev0

EinsteinPy Geodesics is an addon package for EinsteinPy, that wraps over Julia's excellent DifferentialEquations.jl suite and provides a python interface to solve for geodesics in Kerr & Schwarzschild spacetime EinsteinPy is an open source pure Python package, dedicated to problems arising in General Relativity and Gravitational Physics As with EinsteinPy, EinsteinPy Geodesics is released under the MIT license.

Documentation

Complete documentation for this module can be accessed at https://docs.geodesics.einsteinpy.org/en/latest/ (Courtesy: Read the Docs).

Requirements

EinsteinPy Geodesics requires Python >= 3.7, Julia >= 1.5 and the following Julia packages:

Julia
- DifferentialEquations.jl >= 6.15
- ODEInterfaceDiffEq.jl >= 3.7

Installation

First, ensure that, Julia is installed in your system and added to PATH. See https://julialang.org/downloads/platform/ for platform specific binaries and installation instructions. einsteinpy_geodesics also requires DifferentialEquations.jl and ODEInterfaceDiffEq.jl. You can add them, like so:

$ julia
julia> using Pkg
julia> Pkg.add("DifferentialEquations")
julia> Pkg.add("ODEInterfaceDiffEq")

Finally, einsteinpy_geodesics can…

View on GitHub

On top of this, I also overhauled the geodesic plotting module and added support for 3D animations, parametric plots and choice of spatial coordinates in 2D plots, in both Static and Interactive modes (that use matplotlib and plotly respectively). I present some of the plots, produced through the final API, below. The plots shown here, have a mix of both Static and Interactive back-ends, as well as time-like and null-like geodesics.

Kerr Null-like Geodesic
Kerr Time-like Geodesic
Kerr Frame Dragging
Schwarzschild Precession
Schwarzschild Time-like Closed Orbit
Schwarzschild Time-like Closed Orbit Parametric Plot

Until next time...

The EinsteinPy geodesics API currently provides a choice of solvers, between a python back-end and a julia back-end, through the optional einsteinpy_geodesics add-on. I will continue to work on improving the python back-end, but for now, einsteinpy_geodesics adds proper & accurate geodesic calculations to EinsteinPy, in the Kerr family of spacetimes (that includes Schwarzschild). Also, a notable aspect of the HamiltonianProblem approach is that, in principle, it should be easily extensible to Kerr-Newman geodesics, which is something, I'd like to explore, as soon as my GSoC commitment is over. I have another short blog coming up, that explains how to use the API (and has more cool plots), that will probably also be the last GSoC blog. Till then, I leave you with a nice animation, created entirely with EinsteinPy.

(Extremal) Kerr Time-like Constant Radius Orbit

References:

[1]: Carter, Brandon; Global Structure of the Kerr Family of Gravitational Fields, 1968 , Physical Review, 174(5), pp. 1559-1571

[2]: Misner, Charles W. and Thorne, K.S. and Wheeler, J.A; Gravitation, 1973, W. H. Freeman, ISBN: 978-0-7167-0344-0, 978-0-691-17779-3

GSoC 2020: Blog 4 - Update on Null Geodesics in Kerr Spacetime

Jyotirmaya Shivottam — Thu, 06 Aug 2020 10:09:20 +0000

Progress so far...

Support for calculation and graphing of Null Geodesics in Kerr (and by extension, Schwarzschild) spacetime is nearing completion (PR #527). Last week, I hit a serious obstacle, related to maximum floating point precision and accumulation of numerical errors (which is also the reason for the delayed blog). Since the Geodesic Equations are stiff ODEs, small instabilities can wreak havoc on step-size control and completely destabilize the solution. I observed this happening with my code for Null Geodesics. As the light ray approaches the black hole, the integrator can no longer choose a proper step-size and the solution becomes inaccurate. In this blog, I will be discussing this issue and how we are approaching it with the new Null Geodesics module. I also present some of the null geodesic plots, created using this module.

Stiff ODEs are evil!

Stiff ODEs and Numerical Methods have always been at loggerheads. There is no precise definition for stiff ODEs, but an important feature is that, they are prone to become unstable. The usual solution is to choose a solver, that can accommodate very small step-sizes, while keeping overall error low. SciPy provides performant wrappers for LSODA/BDF methods, that are usually suitable for stiff systems, but in our case, these methods are unhelpful, as can be seen in the image below. For comparison, I have used Mathematica to obtain geodesics for the same conditions. The only major difference, here, is the solver. The plot on the left is Mathematica-generated, while the plot on the right was generated by Python. Note that, all the plots in this post have their axes normalized to the gravitational radius, or units of $GMc2\frac{GM}{c^2}$ .

The initial conditions for the plots above are as follows. The timelike component of the initial velocity was calculated by setting $g_{ab}u^au^b = 0$ .

a = 0.9
end_lambda = 200
max_steps = 200
position = [0, 20., pi / 2, pi / 2]
velocity = [-0.2, 0., 0.002]

Other solvers (that are suited to non-stiff problems) become unstable long before the desired number of integration steps is reached. Given the lack of a proper solver, I wrote my own solver, using a step-size control scheme from the venerable Numerical Recipes (Press et al, 2007), fine-tuned to the problem. Sadly, this did not produce better results and it even failed for certain pathological higher-order orbits. Here, "high-order" implies "loopy" orbits, very close to the black hole, while "pathological" can encompass higher-order orbits to orbits, that are scattered at large angles (i.e., orbits, with sharp turning points, à la the plots above).

Then, I set out to find the reason behind the instability. Based on my tests, the stiffness comes from the singular nature of the black hole horizon (in Boyer-Lindquist coordinates), which can force the solver to choose incredibly small step-sizes, which in turn leads to more and more floating point error and over large intervals, the obtained solution becomes completely unphysical. This is what, "unstable" means here. Apart from the graphical representation of the instability through the plots, we can also see the instability numerically, through the norm of 4-Velocity of the light ray, as it evolves:

All of these values should be ~0 and on comparing with the plot on the right, it is easy to see, that the norm becomes too high as the light ray gets closer to the black hole. This tells us, that a correlation exists between the initial conditions and the instability, which is expected.

In discussions with my mentors, we explored a few solutions, such as, using another system of units or coordinate system. However, we are already using the most suitable unit and coordinate systems for numerical computation of geodesics - M-Units and Boyer-Lindquist Coordinates. I should note here, that at slightly larger initial radial distances and speeds, the code provides a good approximation to the actual solution, as can be observed in the plot and table below.

Clearly, the accumulated error over lambda is smaller in this plot.

The initial conditions for the second set of plots are as follows:

a = 0.9
end_lambda = 200
max_steps = 200
position = [0, 30., pi / 2, pi / 2] # Only difference
velocity = [-0.2, 0., 0.002]

Until next time...

Initially, I had planned to develop the Null Geodesics module, such that simulating a photon sheet would be possible through this module itself. The purpose would be applications in radiative transfer calculations, which require simulation of pathological orbits for better approximations in the strong gravity regime. But the issue of error accumulation has made it difficult to continue with this strategy. We have decided to make the current code merge-ready, while keeping the PR open, mainly because, the code performs well at larger initial distances. I have already made relevant changes to ensure the code is merge-ready. The status of the PR can be viewed at PR #527.

I am currently mulling the option of implementing the code in some low level language and then, building a wrapper around it, goal being to achieve better error-control, which I have not been able to obtain with Python/SciPy. Another approach that I am considering, is to restrict integration near the event horizon, based on step-size changes. Since multiple options are being explored and this is the last coding period, I have decided to make these blogs weekly. So my next blog should be up, next Friday. Hopefully, I will have solved this by then.

GSoC 2020: Blog 3 - Coordinates & Perturbations

Jyotirmaya Shivottam — Thu, 16 Jul 2020 20:11:49 +0000

Progress so far...

The refactor of the coordinates module of EinsteinPy was completed over the last week and PR #521 was merged. With this, EinsteinPy now supports 4-Position and 4-Velocity calculations. To better reflect these additions across EinsteinPy, a few changes were also made to einsteinpy.metric and einsteinpy.geodesic modules. All core modules of EinsteinPy are well-integrated now. As an aside, new metric and coordinate classes can now be added to the numerical side of EinsteinPy and we welcome community contributions for the same. Please refer to #525 for more details.

In my last blog, I had mentioned that this PR would also add support for Kerr-Schild (KS) coordinates and KS forms of the vacuum solutions. However, after discussions with my mentors, we have rejigged the schedule a bit and we are now prioritizing the implementation of Null Geodesics calculation and Radiative Transfer functionality in EinsteinPy, with the addition of KS coordinates becoming optional now. This has been put up as issue #526, and as before, community contributions are welcome. Nonetheless, as promised, let’s discuss what KS perturbations are and why KS coordinates are so interesting.

A case of hard differential equations

Einstein's Field Equations (EFE), at their core, describe a field theory for gravity, quite similar to how Maxwell's Equations define a field theory for Electromagnetism. In Electromagnetism, electric and magnetic fields are basic physical entities, which both affect and are affected by a matter property called charge. Similarly, in GR, there is an underlying structure of spacetime, quantified by the metric tensor, which is calculated using EFE. These equations are second order, non-linear differential equations, which makes finding solutions for general matter-energy distributions rather difficult.

However, let's say, we are only interested in cases, where the deviation from a known distribution is "small", for example, when we go from a perfectly spherical body to a spheroidal body. In these cases, we can look for an approximate "perturbed" solution (metric) to this "perturbed" distribution by writing the approximate solution as $gab(λ)=0gab(λ)+hab(λ)g_{ab}(\lambda) = {^0g_{ab}}(\lambda) + h_{ab}(\lambda)$ , where $λ\lambda$ is a parameter, ${^0g_{ab}}$ is the known exact solution for the underlying distribution and $h_{ab}$ is a linear perturbation ("Linear", as $g_{ab}$ depends on the first power of $h_{ab}$ ). Kerr-Schild perturbations are a way of linearizing EFE or to form linearized approximation to EFE. Note that, EFE are in general, non-linear. So, these approximations usually only describe the solution till the first order and higher order corrections become necessary as systems grow in complexity. However, Kerr-Schild class of perturbations is unique as all higher order corrections vanish.

Kerr-Schild Perturbation

A generalized Kerr-Schild perturbation has the form:

h_{ab}^{KS}(\lambda) = V l_a l_b

where,

V

is some scalar and

l

are one-forms or functionals, defined with respect to the underlying or background metric,

{^0g_{ab}}

. A functional can be understood as a usual function, which outputs a number, but instead of accepting numbers as inputs, it takes in functions. Now, the form of

h_{ab}^{KS}

might seem arbitrary. However, there are some reasonable constraints on this form. First among them is, for

λ→0\lambda \to 0

, we must recover

{^0g_{ab}}

. Also,

l

must be null with respect to

{^0g_{ab}}

, i.e.

gablalb=0{^0g_{ab}}l_a l_b = 0 \implies {g_{ab}}l_a l_b = 0

. And, depending on the spacetime, there are additional properties that must be satisfied.

Due to the strictly linear nature of KS perturbations, we obtain some nice features for the perturbed metric, that can even be linked back to the background metric, ${^0g_{ab}}$ . For example, the causality of an arbitrary vector, $v^a$ is guaranteed to remain the same across ${^0g_{ab}}$ and $g_{ab}$ . This implies timelike and spacelike vectors remain timelike and spacelike, respectively, even in the perturbed spacetime. This also extends to null geodesics, meaning that the two metrics have shared light rays. This property, in particular, is incredibly useful from the perspective of my project. But a more notable property is that, KS perturbations impose a structure on the Stress-Energy Tensor, $T_b^a$ .

These properties make the KS forms of the metric particularly suited to describe vacuum solutions, on a flat underlying metric, such as Schwarzschild and Kerr solutions. In these cases, a particular set of coordinates, called the (cartesian form of) Kerr-Schild coordinates, are chosen to define $V$ and $l$ . These coordinates are usually more involved than the most basic coordinate for the spacetime (such as Schwarzschild or Boyer-Lindquist Coordinates), but they have an interesting property, in that they expose the singularity structure well, with no coordinate singularities.

More interesting spacetimes, such as pp-waves and anti-de Sitter space can also be generated through KS perturbations. Not all solutions however, can be linearized to a KS form, as we only have three degrees of freedom in varying $V$ and $l$ . Moreover, as is usual with perturbation techniques, it is difficult to estimate how small $λ\lambda$ must be, so that the approximation has sufficient accuracy. In any case, Perturbation Theory is a powerful mathematical tool, that most physicists use on a regular basis. A (surprisingly) simple and general introduction can be found at Wikipedia.

Until next time...

While I am not working on KS perturbation functionality for EinsteinPy, at the moment, I have opened an issue (#526) for it. Addition of KS perturbed forms of metrics to EinsteinPy is relatively easy now, due to the structure of the Metric classes in einsteinpy.metric. Feel free to take a look and contribute. Currently, I am working on the crux of the project - Null Geodesics and Radiative Transfer in Kerr spacetime, figuring out the equations and validating them for use with SI units, with proper unit conversions. Once this is done, I will start putting them into code and hopefully, we will soon have this new functionality ready to be merged into EinsteinPy and open to people for use. I expect my next blog to be on Null Geodesics.

GSoC 2020: Blog 2 - Parameters and incompatible units

Jyotirmaya Shivottam — Tue, 30 Jun 2020 09:54:19 +0000

Progress so far...

Over the last two weeks, PR #512 was finalized and merged into EinsteinPy. With it, the big refactor of the metric module is over. As I had mentioned in my last blog, this PR adds a core structure to the metric module, with the metric.BaseMetric class, and adds some new functionalities, like support for Kerr-Schild Perturbations. Also, 7 issues, big and small, were fixed in this PR, ranging from purely semantic issues to mathematical inaccuracies. One of the important ones was Issue #514, the distinction between Rotational Length Parameter, $α\alpha$ and Dimensionless Spin Parameter, $a$ , which has had a cascade effect across the codebase. It is worth a discussion, as it has led to many changes across the numerical side of EinsteinPy.

What's all the fuss about?

The expressions for the Rotational Length Parameter, $α\alpha$ , and Spin Parameter, $a$ , are as follows:

\alpha = \frac{J}{Mc}

\frac{J}{J_{Max}} = \frac{Jc}{GM^2}

which give us:

\boxed{ \alpha = \frac{GM}{c^2}a }

Here,

J, M

denote the angular momentum and mass of the gravitating body (for example, black holes), respectively, while

G, c

denote the Gravitational Constant and speed of light in vacuum, respectively. As the expression for

a

tells us, it is a constraint on the body's angular momentum. It might seem arbitrary at first, as in classical physics, we do not bother with such constraints and

J

is not limited. However, in General Relativity, especially in the context of how black holes are formed, or more appropriately, how we hypothesize their formation, we find that, for certain solutions to the Einstein Field Equations, like the Kerr solutions (rotating black hole), a maximum limit needs to be placed on the black hole angular momentum. This is a direct consequence of the efforts to contain the causality-breaking nature of black hole singularities, especially a naked singularity. As the name suggests, naked singularities are not covered by an Event Horizon and so, they would be observable from outside the black hole. This breaks the framework of GR, as no deterministic predictions about the future evolution of spacetime can be made near a singularity. Note that, GR breaks for covered singularities as well, but causality remains intact, because the spacetime within the event horizon is not viewable or approachable from outside. To complicate things further, GR, as is, does not preclude the existence of naked singularities. All of this led to much deliberation on this issue and ultimately, in 1969, Sir Roger Penrose postulated the (Weak) Cosmic Censorship Hypothesis. This hypothesis aims to alleviate the issue with physical singularities by claiming that naked or observable singularities cannot exist in the universe, as the deterministic nature of general relativity will fall apart otherwise. It is worth mentioning that, while this is still a conjecture, we are yet to detect naked singularities through astronomical observations. Also, in order to work in the theoretical framework of GR and make physical predictions using GR, one needs to circumvent the problem of naked singularities and this is where

a

comes in for rotating black holes, like in the Kerr solution. To ensure that, the black hole singularity has an event horizon, we must constrain

J

and this is done by the spin parameter,

a

. Here, for

\leq \frac{GM^2}{c}

, the solution has an event horizon, while for

\frac{GM^2}{c}

, the solution possesses naked singularities, which we must exclude, due to reasons discussed above. This condition is obtained by solving for singularities in the Kerr solutions.

As for $α\alpha$ , it is a convenient length parameter, that helps in making equations shorter to write.

Physics aside, the issue, as mentioned in #514, was that the modules in EinsteinPy do not differentiate between $α\alpha$ and $a$ . This is due to the intermixing of Geometrized quantities with that in SI units. In geometrized units, $G = c = 1$ , making $α=a=JM\alpha = a = \frac{J}{M}$ . However, in SI, this is not the case. EinsteinPy works purely in SI and this issue is an artifact from the changes around unit conversions. This issue has been fixed now, after #512, with $a$ being one of the defining parameters in rotating spacetimes, while $α\alpha$ is calculated using $a$ .

As I mentioned earlier, this issue has caused a ripple effect across the numerical side of EinsteinPy. It has highlighted some problems, around systems of unit, in the definition of a few core methods in the metric classes, as well as the coordinates module.

Until next time...

I am currently working on fixing the aforementioned issues, as well as refactoring the coordinates module. To be specific, I am working on adding support for 4-Vector calculations and Kerr-Schild coordinates to the coordinates module. The latter would make use of the Kerr-Schild perturbation functionality of the metric module and help in modularizing the definition of new metric classes. I will discuss Kerr-Schild perturbations and the so-called Kerr-Schild form of the metric tensor in my next blog.

In the meantime, here's a piece of trivia, related to the Cosmic Censorship Hypothesis.

In 1991, John Preskill and Kip Thorne bet against Stephen Hawking that the hypothesis was false. Hawking conceded the bet in 1997, due to the discovery of some special situations, that violated the hypothesis, which he characterized as "technicalities". Hawking later reformulated the bet to exclude those technicalities. The revised bet is still open, with the prize being "clothing to cover the winner's nakedness".

The description of the original bet can be found here, while that of the new bet can be found here.

GSoC 2020: Blog 1 - Beginning of Coding Period

Jyotirmaya Shivottam — Sat, 13 Jun 2020 18:20:19 +0000

So the community bonding period of GSoC has ended and the coding period is officially underway. In my last blogpost, I had outlined the basic principles of General Relativity that go into my project. I had also mentioned, that the next blog will have details about the coding process. However, things had to be slowed down considerably, due to the announcement of closure of my academic session and some logistical issues. This has also affected my blog schedule, as I could not work on a blog, that was supposed to be up 2 weeks ago. However, I am pleased to inform that all the issues are sorted out now, leaving the rest of the summer free for me to delve into the project. Whew. I am also exceedingly grateful to my mentors, who have been understanding throughout. As for details on the code implementation for my project, I have decided to break it up across the blogs, as "Progress Reports" (bland, I know), in order to provide a better understanding of both what I am working on and my approach to it. So, read on.

Progress so far...

Over the last few discussions with my GSoC mentors, we have discovered some logical bottlenecks in some of the EinsteinPy modules, especially the way the metric & utils modules work at the moment. Currently, the utils module stores most of the necessary functions required to form a metric class, while the metric module lacks support for user-defined metrics. It also handles the calculation of particle trajectory, which logically belongs inside a geodesic module. Since my project is on Null Geodesics, these issues are major obstacles, that must be overcome, before the work on adding support for Null Geodesic calculation begins.

As such, I am refactoring these modules, so that we have logical cohesion across EinsteinPy, whilst also adding some new features, like a brand new metric class, that supports defining arbitrary metrics and also adding first order linear perturbations to the metric, written in Kerr-Schild form. I have also grouped together the utility functions present in utils in metric itself. These changes can be followed at the PR link, here. At the moment, I am reusing some of the old tests and writing some new ones for the new features. I hope to see this PR clear all tests and be merged soon.

Until next time...

After metric, comes the coordinates module, which will see some changes and code rearrangements too, albeit not as much as metric and utils. The basic work on this has already started and after the first PR is merged, I will open another PR with these changes.

We have lost around 2 weeks, due to the aforementioned issues on my side. Therefore, it is contingent on me to accelerate the pace of development now. Over the next few weeks, we should see some cool new feature additions to EinsteinPy. I will be detailing them here, as we proceed with the project.

GSoC 2020: Blog 0 - A Brief History...

Jyotirmaya Shivottam — Mon, 18 May 2020 22:14:41 +0000

>>> print("Hello, World!")
Hello, World!
>>>

Hi, I am Jyotirmaya Shivottam, a Physics Undergrad at National Institute of Science Education and Research, Bhubaneswar. I am thrilled to have been selected in Google Summer of Code 2020 to work on implementing a module for Null Geodesics in Kerr and Schwarzschild Spacetimes in EinsteinPy (EPY), under OpenAstronomy. This will (tentatively) be an 8-part blog series.

A Brief History of my brush with Numerical Relativity and EPY

I had developed an interest in Cosmology and in particular, Einstein's Theory of Relativity, after taking a course in the Special Theory of Relativity, at my university. The course content and presentation intrigued me. So, I decided to read up more on the subject and found the mathematical description of spacetime in General Relativity immensely satisfying. But soon after, it became difficult to visualize the complicated geometries. I am used to understanding things by putting them in code. So naturally, I endeavoured to do the same here. After a brief period of trying and failing to solve it myself, in true programmer fashion, I started looking up software projects that could help me. I was surprised to find that most projects were either in severely outdated FORTRAN77 (yikes!), or in CUDA C, which gets messy to set up. Then, I found EinsteinPy, which is a python package for Numerical Relativity. Their GitHub Repo had listed an issue quite similar to what I wanted to solve. I also came to know that they were participating in GSoC this year and so, I picked up and solved an issue and opened a PR. Thus began my voyage into the world of Numerical Relativity and Free & Open Source Software.

But enough about me. Let me introduce you to the project, that I will be spending this summer on.

Null Geodesics - Is it another Java Exception? (Spoiler Alert: NO!)

As mentioned before, my project is on implementing Null Geodesics in Kerr and Schwarzschild Spacetimes. Geodesics can be understood as a generalization of the concept of shortest paths, between two points, to higher dimensions. In Euclidean or flat geometries, a geodesic is a straight line and in trivial non-Euclidean geometries, like a thin spherical shell, it is an arc of a great circle. Around 115 years ago, Einstein, through his Special Theory of Relativity and then, through its successor, ten years later, the General Theory of Relativity (GR), unified our understanding of space and time into a 3+1-Dimensional spacetime, where many Euclidean axioms and equations no longer work, partly due to the special behaviour of the time dimension (hence, 3+1D). GR, in particular, sheds a lot of light on the relation between the geometry of spacetime and the distribution of matter in it and how they affect each other. Here, John Archibald Wheeler's famous quote on GR comes to mind:

"Spacetime tells matter how to move; matter tells spacetime how to curve".

This correspondence, in mathematical form, is given by Einstein's Field Equations (EFE):

\boxed{R_{\mu\nu} - {\tfrac{1}{2}}Rg_{\mu\nu} + \Lambda g_{\mu\nu} = {\frac {8\pi G}{c^{4}}}T_{\mu\nu}}

Where, you can understand the left-hand side as describing the spacetime curvature (geometry), while the right-hand side describes the matter-energy content of the spacetime. A helpful analogue might be Gauss' Divergence Theorem from Electromagnetism: $∇V=−ρ/ϵ0\nabla V = -\rho/\epsilon_0$ , where the left-hand side describes the structure of a scalar electric potential field, while the right-hand side describes the charge distribution. The potential field tells charges how to arrange themselves, while charges tell how the potential should vary.

The equation above is a Tensor Equation, wherein we solve for $gμνg_{\mu\nu}$ , the spacetime Metric Tensor or simply, the Metric. $gμνg_{\mu\nu}$ defines the geometry of the spacetime. The aforementioned Kerr and Schwarzschild Spacetimes are in fact, Kerr and Schwarzschild Metric Tensors. These fall under a special class of EFE solutions, namely "Vacuum Solutions", because they describe a spacetime, devoid of any matter or non-gravitational fields. The Kerr Metric describes the spacetime around an uncharged, rotating, axisymmetric Black Hole, while the Schwarzschild Metric is a particular simplification of Kerr Metric to non-rotating Black Holes.

Now that we have our spacetime geometry set up, we can think about how the motion of massive objects (read: having non-zero mass, e.g. neutrons) and massless objects, like photons, that make up light, would look like in the spacetime. Here's where, geodesics come in, now in 3+1D. They characterize the paths of such objects, given the metric. The corresponding Geodesic equation can be derived using the metric and the Equivalence Principle or through Variational Methods and is given as:

\boxed{ {\mathrm{d}^{2}x^{\mu} \over \mathrm{d}\lambda^{2}} + \Gamma^{\mu}_{\alpha \beta}{\mathrm{d}x^{\alpha} \over \mathrm{d}\lambda}{\mathrm{d}x^{\beta} \over \mathrm{d}\lambda} = 0 }

Here, $xαx^\alpha$ is 4-position of the object in motion (Position vector in a 3+1D spacetime), $Γαβμ\Gamma_{\alpha\beta}^{\mu}$ or the Christoffel Symbols are a convenient way of writing metric derivatives (derivatives of $gμνg_{\mu\nu}$ ) and $λ\lambda$ is called an Affine Parameter, that is helpful in describing a general geodesic. You might ask, if we are talking about a moving object, why not use time (in particular, "Proper Time", $Δτ\Delta\tau$ ), instead of some arbitrary $λ\lambda$ parameter. The reason behind this and the importance of this parameter will be apparent soon.

We have understood, what spacetime and geodesics refer to (even though, we glossed over most details :P), but what are Null Geodesics? Simply speaking, all massive objects follow time-like geodesics and all massless objects, like light, follow null-like geodesics. Before going into the distinction, we need to know, what "Proper Time" means. Proper Time ( $Δτ\Delta\tau$ ) is the time measured by a clock in an object's rest or proper frame. This is the time physically experienced by the object. Now, the distinction between paths, followed by massive and massless particles, or correspondingly, "time-like" and "null-like" geodesics, arises due to the causal structure of spacetime, in which the constancy of Speed of Light in vacuum plays a major role. Light or massless particles always move at the speed of light, at 299,792,458 m/s. This implies that we cannot construct a reference frame, where a clock is at rest with respect to the photon or massless particle, as no clock can move at that speed. Thus, no proper time ever elapses for light rays or massless particles, i.e. $Δτ=0\Delta\tau = 0$ . This is also why, we cannot use $Δτ\Delta\tau$ in the geodesic equation to obtain the paths, traced by massless particles. So, we seek another parameterization and that is given by the affine parameter, $λ\lambda$ . Do note that, for massive objects, $Δτ\Delta\tau$ works just fine, as it is well-defined there.

And so, we have assembled all the parts for getting a basic understanding of "Null Geodesics in Kerr and Schwarzschild Spacetimes". In my GSoC project, I will be implementing a numerical solver for the geodesic equation for null geodesics and I will use the computed geodesics for Radiative Transfer calculations. At the moment, the Community Bonding Period is underway. My mentors and I have set up a weekly meeting schedule, where we brainstorm on the project specifics. An outcome of those meetings is an EPE (EinsteinPy Proposal for Enhancement), where you can find the details on my proposed implementation and coolest of all, you can help shape the implementation, if you are interested. Isn't open source awesome? (Yes, it is!)

What now?

I plan to document my progress with the project here, with biweekly blogs. While this blog was physics-heavy, the plan for the next blog is to post a distilled summary of the actual code implementation. Nonetheless, I hope, you found this small introduction into the domain of Relativity interesting. If you wish to play around with some numerical relativity (especially during this pandemic, when most are home and bored), you can try out EinsteinPy. For now, I will leave here with a small tidbit:

If you've watched Chris Nolan's brilliant movie, Interstellar, and wanted to simulate and obtain a Black Hole "image" similar to Gargantua, all you would need to solve is the geodesic equation alongside a couple more complicated equations. A brilliant resource to do that is a Mathematica Notebook, hosted on Dr. Kip Thorne's personal Caltech homepage.