Simulating Light in Curved Spacetime: Relativistic Raytracing

Based on the work by h1m3

Objective: To render a black hole that is physically accurate and strictly adheres to the laws of General Relativity.

To achieve this, the process is broken down into three primary phases:

1. Derivation: Establishing the equations of motion for light rays.
2. Integration: Implementing a Runge-Kutta integrator within GLSL.
3. Rendering: Creating a fragment shader that traces rays from the camera backward in time.

---

🛠️ Preparation: Metrics, Connections, and Geodesics

(Placeholder for colored_background.png)

At its most basic level, light travels in a straight line. However, "straightness" is relative to the coordinate system and the geometry of the space:

• Coordinate Distortion: In a polar system, a straight line looks linear is described by non-linear equations.
• Geometric Distortion: On a non-flat manifold (such as a torus), no coordinate system exists where a straight line appears linear.

(Placeholder for thorus_geo.png)

What is a Geodesic?

To define a "straight" path on a curved surface, we use the concept of a geodesic. A curve is considered a geodesic if a vector pointing along the curve remains parallel to itself as it is transported along that path. In essence, a geodesic is the closest approximation of a straight line on an arbitrary manifold.

Deriving the Equations via Python

Calculating geodesic equations for a 4D manifold is tedious by hand. Instead, we can utilize the sympy library, specifically the sympy.diffgeom module. While this module can be somewhat clunky and inconvenient, it is effective for deriving equations on manifolds defined by an atlas.

For this project, we use a single patch with two distinct coordinate frames:

1. Spherical: Ideal for the Schwarzschild metric.
2. Cartesian: Necessary for the camera's perspective.

Step 1: Coordinate Definitions

First, we define our symbols and the transformation logic from spherical to Cartesian coordinates:

from sympy import *
from sympy.diffgeom import *
import numpy as np

# Terminal symbols
t_, x_, y_, z_ = symbols( 't x y z' )
rho_, theta_, phi_ = symbols( '\rho \theta \phi' )

# Transformation equations
x = rho_ * sin(theta_) * sin(phi_)
y = rho_ * sin(theta_) * cos(phi_)
z = rho_ * cos(theta_)

print(latex(x_) + " = " + latex(x))
print(latex(y_) + " = " + latex(y))
print(latex(z_) + " = " + latex(z))

This results in the following LaTeX representations: $\begin{array}l x = \rho \sin{\left(\phi \right)} \sin{\left(\theta \right)} \\ y = \rho \sin{\left(\theta \right)} \cos{\left(\phi \right)} \\ z = \rho \cos{\left(\theta \right)} \\ \end{array}$

Step 2: Manifold Construction

We then establish the 4D spacetime manifold and link the coordinate systems:

spacetime = Manifold( spacetime , 4)
patch = Patch( patch , spacetime)
relation_dict = { ( spherical , cart ):[(x_, y_, z_, t_), (x, y, z, t_)] }
coord_sh = CoordSystem( spherical , patch, (t_, rho_, theta_, phi_), relation_dict)
coord_cart = CoordSystem( cart , patch, (x_, y_, z_, t_), relation_dict)

Note: sympy.diffgeom requires using the scalars returned by the CoordSystem rather than the original symbols:

c_t, c_rho, c_theta, c_phi = coord_sh.base_scalars()
# Output: \mathbf{t}, \mathbf{\rho}, \mathbf{\theta}, \mathbf{\phi}

Step 3: The Schwarzschild Metric

We define the metric tensor, which describes the geometry of spacetime around a non-rotating massive object.

Rs = symbols( r_s )
d_t, d_rho, d_theta, d_phi = coord_sh.base_oneforms()
TP = TensorProduct
metric = ( (1 - Rs/c_rho) * TP(d_t, d_t) 
           - 1 / (1 - Rs / c_rho) * TP(d_rho, d_rho) 
           - (c_rho**2) * TP(d_theta, d_theta) 
           - (c_rho**2) * (sin(c_theta)**2) * TP(d_phi, d_phi))

The resulting metric is: $\left(- \frac{r_{s}}{\mathbf{\rho}} + 1\right) \operatorname{d}t \otimes \operatorname{d}t - \sin^{2}{\left(\mathbf{\theta} \right)} \mathbf{\rho}^{2} \operatorname{d}\phi \otimes \operatorname{d}\phi - \mathbf{\rho}^{2} \operatorname{d}\theta \otimes \operatorname{d}\theta - \frac{\operatorname{d}\rho \otimes \operatorname{d}\rho}{- \frac{r_{s}}{\mathbf{\rho}} + 1}$

The Levi-Civita Connection

To find the geodesics, we need the Levi-Civita connection (the Christoffel symbols of the second kind).

• What it does: It defines how vectors in a tangent space transform as they move across the manifold.
• Key Property: It preserves the metric tensor, meaning the covariant derivative of the metric tensor is zero.

We can automate the derivation of these symbols using the metric_to_Christoffel_2nd function:

christoffel = simplify(metric_to_Christoffel_2nd(metric))
# Sorting and printing non-zero symbols...

This yields the following components: $\begin{cases} \Gamma^{1}_{00} =\frac{r_{s} \left(- r_{s} + \mathbf{\rho}\right)}{2 \mathbf{\rho}^{3}} \\ \Gamma^{0}_{01} =\frac{r_{s}}{2 \left(- r_{s} + \mathbf{\rho}\right) \mathbf{\rho}} \\ \Gamma^{0}_{10} =\frac{r_{s}}{2 \left(- r_{s} + \mathbf{\rho}\right) \mathbf{\rho}} \\ \Gamma^{1}_{11} =\frac{r_{s}}{2 \left(r_{s} - \mathbf{\rho}\right)} \end{cases}$

Summary of Implementation Progress

• Define the 4D Spacetime Manifold.
• Establish Spherical and Cartesian coordinate patches.
• Define the Schwarzschild Metric tensor.
• Derive Christoffel symbols via the Levi-Civita connection.
• Implement Runge-Kutta integrator in GLSL.
• Finalize the backward raytracing fragment shader.