Ray-tracing a Black Hole Image

This is a project that I completed as a part of my Physics minor coursework (PHY 427 Computation Physics) at Washington University in St. Louis.

The Boyer-Lindquist method used here provides a metric for the coordinates of a Schwarzschild black hole, and is the same method used to create the black holes in Interstellar:

Most of the heavy computation is performed in C++. This part of the project has no dependencies outside of the standard library, so running it is pretty trivial.

The rendering for the Einstein rings and unstable spherical orbits is done within a Jupyter notebook, but the rendering of the main black hole image is done somewhat naively with PPM files (See Section 4.2).

Note

What follows are my original submission notes for the project, which are probably not very helpful without any of the context.

1) Methods

1.1) Boyer-Lindquist Metric

The Boyer-Lindquist metric was initially implemented using derivatives which were obtained via hand- computation. While this provided successful results for some of the validation examples, namely the Einstein Ring and the stable spherical orbit A, all other calculations were wildly incorrect. Some or all of this error could be caused by the reuse of previously computed values in the calculations of derivatives whenever possible.

Unable to find the error in the manually obtained derivatives, I implemented numeric differentiation, however this solution, while correct, resulted in significant error. Reducing the finite difference $h$ in an effort to obtain a tolerable level of error significantly increased computation time. This numeric solution is implemented in boyer_lindquist_metric.h.

Unsatisfied with the large error and slow computation time provided by the numerical differentiation, I found Maxima, a free and open-source alternative to Mathematica (to which I do not have access), was able to provide exact symbolic representations of the necessary derivatives required by the Boyer- Lindquist metric. Maxima also conveniently provides a method to export symbolic expressions to formats compatible with a variety of programming languages, including C++, and these are included directly in boyer_lindquist_metric.h.

1.2) Ordinary Differential Equation Solver

An adaptive 4th order Runge-Kutta ODE solver was selected for its ease of use, high-accuracy, and efficient computation. Both absolute and relative tolerance were set to $10^{-14}$. The initial step size was set to $0.01$ for all problems. This is implemented in rk4_adaptive.h.

2) Einstein Ring

This validation test was immediately successful (see Figure 1 below) even under when the incorrect analytical derivatives were used. As this example has $u_\theta=\pi/2$ and kept constant, the problematic derivatives were not revealed until later in testing.

Nevertheless, the path of the photon accurately orbited the black hole and returned to the original starting position without significant deviation from the starting position. The stop conditions used for this test were $|r-r_0|<0.1$ and $u_r>0$ to ensure the photon was moving away from the black hole and the stop condition was not raised at $t=0$.

File einsteinRing.cpp generates a CSV file which can be rendered into the images below using the Jupyter Notebook einsteinRing.ipynb.

Figure 1: Einstein Ring. Left: example with total $\phi=2\pi$. Right: example with total $\phi=4\pi$

3) Unstable Spherical Orbits

Figure 2: Orbit C. Right: Three-dimensional plot of photon path. Color indicates the value of $t$. Right: The error term $|r-r_0|/r_0$ over time.

The spherical orbits serve as a great test of the numerical precision and accuracy of the underlying Boyer-Lindquist model. The various combinations of initial conditions validate different parts of the metric. See Table 1 below for the results obtained from spherical orbits.

Orbit	Result	Duration of Stability ($t$)
A	Fell into black hole	~2000
B	Ejected	~175
C	Ejected	~155
D	Ejected	~150
E	Fell into black hole	~200

Table 1: Results of Unstable Spherical Orbits

The duration of stability is the approximate value of $t$ before the error term $|r-r_0|/r_0$ began increasing exponentially.

The calculations for the spherical orbits are implemented in sphericalOrbit.cpp and the 3D render and error plot are created with the sphericalOrbits.ipynb Jupyter Notebook.

4) Ray-traced Black Hole Image

The main event!

4.1) Ray-tracing

The initial conditions are constructed from the relative location of the pixels being rendered onto the image. The position of the virtual “camera” is the center of the screen. The screen then has a fixed size, and resolution is increased instead by increasing the density of pixels per unit of area.

Figure 3 : A selection of photon paths corresponding to a column of pixels in a low-res rendering bending around the black hole (hole enlarged for visibility)

Once the initial conditions are calculated based upon the pixel offset from the camera location, the differential equations are repeatedly stepped through until one of the ending conditions is met:

• The photon goes off to infinity (see Figure 3 above), approximated by $r>D*1.1$, where $D$ is the distance between the black hole and the camera.
• The photon falls into the black hole, determined by $r<r_H^+$, where $r_H^+$ is the event horizon.
• The photon hits the accretion disk.

The first two cases result in a black pixel. The third results in a non-black pixel and is described below.

The complete ray-tracing and rendering procedure is implemented in blackhole.cpp.

4.1.1) The Accretion Disk

The accretion disk is approximated as a disk with inner radius $r_{in}=5M$, outer radius $r_{out}=20M$, and thickness $0.1M$. The disk can be hit head on, in the case of photons approaching the front of the black hole, or by photons which are bent over the top or under the bottom of the black hole.

The brightness of the displayed pixel, $I_{obs}$, is determined by the redshift factor:

$$1+z=\frac{1.0+\Omega\frac{u_\phi}{u_0}}{(-g_{00}-\Omega^2g_{\phi\phi}-2\Omega g_{t\phi})^\frac{1}{2}}$$

$$I_{obs}=\frac{1}{(1+z)^3}$$

4.2) Rendering

The “Portable Pixel Map” or PPM file format is an intermediate photo format intended for the simple creation of images and is ideal for a rendering program such as this. The output file can be created with a simple std::ofstream by outputting RGB triplets for each pixel in ASCII text. This enables the image of the black hole to be created as the calculations are completed.

Although most operating system image viewer programs can open PPM files directly, the inefficient data storage method results in large file sizes, so it is best to convert them to standard image formats. This can be accomplished with the ImageMagick CLI tool:

convert blackhole.ppm blackhole.png

4.3) Observations & Comments

The overage average brightness of the accretion disk seems to correlate directly with. This makes sense, as a faster spinning black hole will result in the orbiting matter to also be spinning faster, and generating more heat, and thus be brighter.

Figure 5: Left: Low-res (128x72) top-down view ($\theta=10^\circ$) of black hole. Right: Low-spin ($a=0.1$) rendering

When viewed from above, the true geometry of the accretion disk is revealed. Simply a flat disk, the upper and lower portions which appear when viewing from the side are “projected” by the bending of light.

Figure 6: Black hole render with $\phi_0=0$

Regardless of the angle $\phi$ of the camera, the position of the inner photon ring doesn’t change. Comparing Figure 4 where $\phi_0=90^\circ$ with Figure 6 where $\phi_0=0$, the photon ring is offset to the right in both cases.

The solid white line which appears in the lower resolution renderings (see Figures 5 and 6) is an artifact of the large step size used. In these cases, for the column of pixels at $x_{sc}=0$, the angle $\beta=0$, and the initial $\phi$-component of momentum, which $\propto\sin{\beta}$ is also zero. This results in the brightness calculation (see 4.1.1) evaluating to zero for these pixels. I have elected to overwrite the color of these pixels to pure white, rather than keep them at black. This issue does not appear in Figure 4, which uses a step size which results in $x$ never being set to $0$.