Black Hole rendering with SageMath

Florentin Jaffredo


This notebook is a step-by-step implementation of a basic rendering engine in curved spacetime. The objective is to obtain a somewhat realistic image of the accretion disk around a black hole.

The technique consists in launching lightlike geodesics toward the past from a single point (the virtual camera), using the geodesic integrator of SageMath. To reduce computation time, the spacetime is assumed be spherical symmetric; this reduces the number of required geodesics to produce an image of $n_x\times n_y$ pixels from about $O\left(n_x n_y\right)$ to $O\left(\sqrt{n_x^2+n_y^2}\right)$.

This work relies heavily on the SageManifolds Project. Advanced SageMath notions will alo be used throughout this notebook, like Cython compilation and multithreading.

This notebook requires a version of SageMath at least equal to 8.5:

In [1]:
'SageMath version 8.7, Release Date: 2019-03-23'


The code is separated into 9 parts.

  • Declaring the spacetime
  • Launching a geodesic
  • Launching a lot of geodesics!
  • Figuring out where it intersects with the accretion disk
  • Adding thickness to the disk
  • Using black-body radiation and convert spectra to rgb
  • First relativistic effect: Doppler effect
  • Second relativistic effect: aberration (forward focalisation)
  • Conclusion


This notebook can be quite ressource hungry to run. For that reason different configurations options are provided. It is recommended to start with the lowest one to check that everything works properly. You can of course adapt the number of CPUs to your needs.

First configuration: will run in less than a minute on a 4-core laptop. Produces tiny images with no details (no secondary image).

In [2]:
# n_cpu = 4 # 4 Go Ram minimum
# n_geod = 100
# nx, ny = 180, 90

Second configuration: will run in about 5 minutes on a workstation, produces a reasonably sized image:

In [3]:
n_cpu = 8 # 8 Go Ram minimum
n_geod = 1000
nx, ny = 720, 360

Third configuration: will run in 30 minutes on the Google Cloud Compute Engine. Produces a 4K image showing tiny details on the secondary disk images.

In [4]:
# n_cpu = 36 # 144 Go Ram minimum
# n_geod = 30000
# nx, ny = 4000, 2000
In [5]:
%display latex

Declaring the spacetime

Let's start slow by declaring the spacetime we'll use for rendering: it is the Schwarzschild spacetime.

It is important to use a coordinate system that is regular at the horizon. Here we use the Eddington-Finkelstein coordinates.

Let $m$ be the masse of the black hole (that we'll take equal to 2 later).

We also add a restriction to ensure that nothing touches the central singularity, and we set the metric $g$.

In [6]:
M = Manifold(4, 'M', structure='Lorentzian')
In [7]:
C.<t, r, th, ph> = M.chart(r't r:(1,+oo) th:(0,pi):\theta ph:\phi')
In [8]:
m = var('m')
In [9]:
g = M.metric()
g[0,0] = -(1-2*m/r)
g[0,1] = 2*m/r
g[1,1] = 1+2*m/r
g[2,2] = r^2
g[3,3] = (r*sin(th))^2
In [10]:

We also define a 3-dimensional Euclidean space $E$ to plot some results, using a map $\phi: M \rightarrow E$:

In [11]:
E.<x, y, z> = EuclideanSpace()
phi = M.diff_map(E, [r*sin(th)*cos(ph), r*sin(th)*sin(ph), r*cos(th)])

Launching a geodesic

Geodesic integration was first implemented in SageMath in 2017 and perfected in 2018 to support fast integration and event handling (used to detect the singularity in our case).

To introduce the method, let's plot an orbit around a black hole.

To do that, we need to find a starting point $p$ as well as an inital velocity vector $v$. It can be quite troublesome to find a suitable one, but here is a free one:

In [12]:
p = M((0, 14.98, pi/2, 0))
Tp = M.tangent_space(p)
v = Tp((2, 0, 0.005, 0.05))
v = v / sqrt(, v))

$v$ is defined as a member of the tangent space at $p$. The last line is used to normalize $v$ as a unit timelike vector.

Next is the definition of the geodesic. We need to pass a symbolic variable for the proper time (which will not be used). The starting point is deduced from the velocity vector (as the point where the velocity vector is defined).

In [13]:
tau = var('tau')
curve = M.integrated_geodesic(g, (tau, 0, 3000), v)

The integration should be very fast. Don't forget to give some numerical value to $m$ here.

In [14]:
sol = curve.solve(step = 1, method="ode_int", parameters_values={m: 2})
# sol = curve.solve(step = 1, parameters_values={m: 2})

Plotting the solution requires an interpolation. This is automatically done in the next line.

In [15]:
interp = curve.interpolate()

The following cell plots the result using the mapping we provided previously. We also add a grey sphere at $r_s = 2m = 4$ (the event horizon) to give a scale.

In [16]:
P = curve.plot_integrated(mapping=phi, color="red", thickness=2, plot_points=3000)
P = P + sage.plot.plot3d.shapes.Sphere(4, color='grey')[1, 1, 1], viewer='threejs', online=True)

You can see that it look nothing like an ellipse, as we are used to in classical celestial mechanics. At this step, you can try adding an angular momentum to the black hole, in other words going from Schwarzschild to Kerr by setting an non-zero angular momentum in the definition of the manifold ($J=1$ works fine). When this is the case, the orbits are not even included in a plane. Don't forget de revert back your changes before proceeding to the next part.

Launching a lot of geodesics!

Of course one geodesic is not enough for us, we'll need at least a few hundred of them.

Because we don't need to compute the equation again each time, what is done is simply copying the previous declaration of the geodesic while changing the initial point and velocity.

It will be usefull here to introduce the Python module multiprocessing and progress bars as widgets:

In [17]:
import multiprocessing
from ipywidgets import FloatProgress
from IPython.display import display

It wouldn't be a great idea to set "1 job = 1 geodesic integration". Indeed, that would mean copying the geodesic declaration a few hundred times, which would be quite slow. What is done instead is seperating geodesics in batches using the following function:

In [18]:
def chunks(l, n):
    """Yield successive n-sized chunks from l."""
    for i in range(0, len(l), n):
        yield l[i:i + n]

The number of batches per CPU in not very important. If set to 1, some CPUs may run faster than other ones and stay idle at the end. If too high, too much time will be spend copying the curve setting. I found 3 to be a good value.

In [19]:
n_batches_per_cpu = 3

We also redefine the previous geodesic to our new needs: fewer steps and the ability to check for chart boundaries when integrating. The $v$ used here will not be used, it will always be overwritten before starting any integration.

In [20]:
curve = M.integrated_geodesic(g, (tau, 0, 200), v, across_charts=True)

When using multiprocessing, functions can only accept a single argument. To overcome this limitation, each argument will be a tuple (curve, start index, number of curve to integrate).

In [21]:
args = []
start_index = 0

for chunk in chunks(range(n_geod), n_geod//(n_batches_per_cpu*n_cpu)):
    args += [(loads(curve.dumps()), start_index, len(chunk))]
    start_index += len(chunk)

The next line prints the list of arguments. We can check that each of the 100 geodesic is correctly set. Our little trick allowed us to only define 13 geodesics (about 3 per core, as we wanted)

In [22]:
(Integrated geodesic in the 4-dimensional Lorentzian manifold M, 984, 16)

Now comes a question: which vector can be used as the starting 4-velocity?

We need a past-oriented lightlike vector pointing toward the center but with a linearly increasing angle. The 3 space components are already imposed. The time component must then be chosen so that the total vector is lightlike.

Let $p$ be the initial point and $v$ the initial 4-velociy, with an unknown time coordinate $dt$ ($y$ depends on the angle, it is a known quantity).

In [23]:
dt, y, r0 = var('dt, y, r0')
In [24]:
p = M((0, r0, pi/2, 0))
Tp = M.tangent_space(p)
v = Tp((dt, -1, 0, y))

The norm of $v$ is currently given by:

In [25]:, v)

We need to find $dt$ so that this expression is equal to 0 (lightlike condition). this is easy:

In [26]:
sol =, v).solve(dt)

As expected, there are two solutions: one past-oriented and one future-oriented. In fact, in our case it does not matter, given that the Schwartzschild spacetime is static.

Next cell defines the function that will be called by multiprocessing. It starts by unpacking the arguments, setting an empty dictionnary as the result, and defining the starting position.

The initial velocity is then overwritten using the formula above, the integration is performed, and the result is added to the dictionnary.

In [27]:
def calc_some_geodesics(args):
    Compute nb geodesics starting at index n0
    global sol, M
    curve, n0, nb = args
    res = {}
    r = 100
    posi = [0, r, pi/2, 0]
    p = M(posi)
    Tp = M.tangent_space(p)
    for i in range(n0, n0+nb):
        # starting vector
        dy = i*0.006/n_geod
        v = Tp([sol[0].rhs()(r0=r, y=dy, m=2).n(), -1, 0, dy])
        # overwrite the starting vector
        curve._initial_tangent_vector = v
        # integration with m=2
        curve.solve_across_charts(step = 0.2, parameters_values={m:2})
        # copy and clear solution
        res[i] = (p.coord(), curve._solutions.copy())
    return res

geo will keep the numerical solutions. I like to see pool as a hole in which I can throw some jobs. multiprocessing will then magically do them for me using every ressource available on the computer.

In [28]:
geo = {}
pool = multiprocessing.Pool(n_cpu)

# progress bar display
f = FloatProgress(min=0, max=n_geod)

for i, some_res in enumerate(pool.imap_unordered(calc_some_geodesics, args)): # do and wait
    # progress bar update
    f.value += len(some_res)
    # update result

# clean exit

If, for any reason, you don't want to use parallel computing, you can replace the previous cell with this one:

In [29]:
# geo = calc_some_geodesics((c, 0, n_geod))

We can now try to visualize those geodesics. Next cell will plot 20 of them.

In [30]:
# add the sphere
P = sage.plot.plot3d.shapes.Sphere(4, color='grey')

# cycle through the solutions
for i in range(0, n_geod, 5*n_geod/100):    
    # set solution
    curve._solutions = geo[i][1]
    # do interpolation
    interp = curve.interpolate()
    # plot the curve
    P += curve.plot_integrated(mapping=phi, color=["red"], thickness=2, plot_points=150, 
                               label_axes=False, across_charts=True)

# show the result[1, 1, 1], viewer='threejs', online=True)

We can see that some falls inside the black hole toward the singularity. That's not an issue because the integration is automaticaly stopped when the geodesic leaves the chart domain defined in part 1.

intersection with the accretion disk

Time to transform those simulated light-rays into an image. To do this, we need to first compute the intersection between each geodesic and the accretion disk.

For this exemple, the disk span from $r=8$ to $r=50$, and is tilted by an angle $\alpha = - \frac{\pi}{20}$.

In [31]:
disk_min = 12
disk_max = 50
alpha = -pi/20

Let's plot the disk on top of the last figure.

(We cheat a little bit here and use a flatten torus).

In [32]:
D = sage.plot.plot3d.shapes.Torus((disk_min+disk_max)/2,
In [33]:
(P+D).show(aspect_ratio=[1, 1, 1], viewer='threejs', online=True)

The same but tilted on the X-axis by an angle $\beta=\frac{\pi}{3}$. As explained earlyer, the final image will be obtained by computing for each pixel :

  • Which geodisic best describe the light-ray
  • Which angle $\beta$ shoulde the disk be tilted
  • The intersection between the disk and that geodesic
In [34]:
(P+D.rotateX(pi/3)).show(aspect_ratio=[1, 1, 1], viewer='threejs', online=True)