Black hole light Simulation

Before

After

To simulate the curvature of space around a black hole we need a function that describes this curvature. In order to keep the problem as simple as possible we remove all constants and complex terms from the equation and reduce it to a 2D problem. As a non-charged non-rotating black hole is a perfect sphere this problem can easily be represented in 2D using Polar coordinates.

\frac{d^{2} u}{d φ^{2}} - \frac{3}{2} \cdot R_{S} \cdot u^{2} + u = 0

where

R_{S}

is the schwarzschild radius of the black hole and

u = \frac{1}{r}

the inverse of the radius from the center of the black hole.

Simulation implementation

Numerical integration

The previously shown ODE can easely be rewritten to look more like a second order ODE:

\frac{d^{2} u}{d φ^{2}} = \frac{3}{2} \cdot R_{S} \cdot u^{2} - u

In Python this can be implemented as follows:

python

def acc(u):
    return 3/2.0 * R_S * u**2 - u

We can integrate this function by using the second order version of the RK4 solver. This solver can be implemented as follows:

python

# second order RK4 integration
def RK4_2(X, V, h, f):
    # RK4 coefficients
    k1_v = h * f(X)
    k1_r = h *   V

    k2_v = h * f(X + 0.5 * k1_r)
    k2_r = h *  (V + 0.5 * k1_v)

    k3_v = h * f(X + 0.5 * k2_r)
    k3_r = h *  (V + 0.5 * k2_v)

    k4_v = h * f(X + k3_r)
    k4_r = h *  (V + k3_v)

    # Update values
    X = X + (k1_r + 2 * k2_r + 2 * k3_r + k4_r) / 6
    V = V + (k1_v + 2 * k2_v + 2 * k3_v + k4_v) / 6

    return X, V

By calling this solver multiple time we can integrate the underlying ODE. The full integration then looks like this:

python

def integrate(u0, phi0, w0, h, numsteps, R_min, R_max):
    """
    u0, phi0, w0 are the initial conditions of the photon.
    h, numsteps are the step size and maximum number of steps to take.
    R_min and R_max are the termination conditions for the simulation.

    Returns an array of the path of the photon in the u, phi plane.
    """
    u, phi, w = u0, phi0, w0
    path = []
    path.append([1/u, phi])
    for i in range(numsteps):
        u, w = RK4_2(u, w, h, acc)
        phi += h
        path.append([1/u, phi])
        
        # Check stopping conditions
        if 1/u < R_min:      break # collision with black hole
        if 1/u > R_max:    break # collision with outer sphere
        
    return np.array(path)

Calculating photon trajectories

In order to calculate the deflection of a photon around a black hole we need to calculate the path of the photon. To simulate what a black hole would look like we need to define a camera that is looking at the black hole. This camera is placed at a distance $r_{0}$ from the black hole. We now simulate the path of a photon that is emitted from the camera at a certain angle $α$ . To do this we need to get the initial conditions of the photon $(u_{0}, φ_{0}, w_{0})$ .

\begin{aligned} u_{0} & = \frac{1}{r_{0}} \\ φ_{0} & = 0 \\ w_{0} & = \frac{u_{0}}{\tan (α)} \end{aligned}

We call the angle by which a photons' path gets changed around a black hole the deflection angle $β$ . It can be obtained by calculating the angle between the photons initial direction and its final direction. Let $\vec{a}$ be the initial photon direction and $\vec{b}$ the final direction then

β = ∠ (\vec{a}, \vec{b}) = \arccos (\frac{\vec{a} \cdot \vec{b}}{| | \vec{a} | | \cdot | | \vec{b} | |})

Note that while the actual deflection angle may be arbitrary high, this function will return deflection angles in the range of $β \in [0, π]$ .

Because the black hole system is rotation symmetric we can calculate these angles once and then use them multiple times. In order to do that efficiently we can use a cache to store the deflection angles.

Calculate the transformation matrix

We want to create a matrix that contains all the information on how certain pixels have to be shifted in order to account for the distortion caused by a black hole. The transformation matrix will be based on equirectangular images. An equirectangular image is an image that, like the most common maps of planet earth, contains the complete surface of a sphere on a rectangular 2d plane. This means that the $(y, x)$ and $(l a t, l o n)$ values are dependent on each other and can easily be converted.

(y, x) \hat{=} (l a t, l o n)

Let the output image also be equirectangular. Then, because we defined the (0,0) coordinate to be at the center of the input and output image, we can easily obtain the viewing direction α as the distance from any position to the center.

α = \arccos (\cos (l a t) \cdot \cos (l o n))

For every $α$ we can calculate the deflection angle $β$ and store it in a cache.

Then we can calculate the angle $α$ for every pixel in the output image and use it to obtain the deflection angle $β$ from the cache.

Now we need the rotation angle $γ$ that represents the rotation of our 2D reference system inside the 3D simulation. If we limit our camera to be on the $z$ axis then $γ$ is the same as the angle between the y-axis up vector $e_{y}$ and the $(x, y, 0)$ values of the $(l a t, l o n)$ transferred to Cartesian coordinates.

γ = \arccos (\frac{e_{y} \cdot (x, y, 0)}{| | e_{y} | | \cdot | | (x, y, 0) | |})

Note that when converting spherical to cartesian coordinates the fact that we defined the center of the image to be (0,0) and not the top results in a rotated system. The conversion from spherical to cartesian coordinates is therefore done by:

spheric2cart (l a t, l o n) = [\begin{array}{c} x = \sin (l a t) \cdot \cos (l o n) \\ y = \sin (l a t) \cdot \sin (l o n) \\ z = \cos (l a t) \end{array}]

[\begin{array}{c} z \\ x \\ y \end{array}] = spheric2cart (l a t + \frac{π}{2}, l o n)

For a given $(l a t, l o n)$ we will now have the deflection angle $β$ and the rotation angle $γ$ . By using these values we can calculate the new lat and lon values $(l a t^{'}, l o n^{'})$ of the input image that are deflected by the black hole to the output image position $(l a t, l o n)$ .

\begin{aligned} l a t^{'} & = l a t + β \cdot \cos (γ) \\ l o n^{'} & = l o n - β \cdot \sin (γ) \cdot sign (l o n) \end{aligned}

Next we have to calculate the pixel value of the input image for $(l a t^{'}, l o n^{'})$ .

python

def angle_to_pixel_rad(lat, lon, img):
    x = (lon + np.pi) / (2 * np.pi) * img.shape[1]
    y = (lat + np.pi/2) / np.pi * img.shape[0]
    
    # change to int
    x = np.array(x, dtype=int)
    y = np.array(y, dtype=int)
    
    # set boundary to be repeating
    x = x % img.shape[1]
    y = y % img.shape[0]
    return y, x

For a complete implementation of this algorithm see the GitHub repository.

If we run this simulation for the milky way galaxy using the following parameters we will get a distorted image of the milky way galaxy.

python

R_S = 1         #   schwarzschild radius of black hole
r_0 = 6         #   distance of camera 
h = 0.01        #   integration precision
numsteps = 2000 #   max integration steps

# calculate deflection angles
num_alphas = 1024 * 12 # number of calculated alphas (view directions)

# transformation params
size_x, size_y = 6000, 3000 # transformation matrix size
zoom = 1                    # allows to only compute center of the image. must be <= 1

We can now change the zoom parameter to only compute the center of the image and get a faster result.

python

zoom = 0.5