Chapter 3: Ray Tracing in the Schwarzschild Geometry

The Schwarzschild black hole is a spherically symmetric, non-rotating, charge-less black hole that was discovered by Karl Schwarzschild in 1915. For these reasons it is the simplest type of black hole. In this chapter we will derive the line metric of the empty space around this black hole’s body, followed by the equations of motion for photon’s trajectories around it. We will then write up a piece of code made in python that contains an algorithm that produces distorted images of these black holes.

3.1 The Schwarzschild Metric

The trajectory of photons around Schwarzschild black holes are governed by its line metric. This line metric is derived below using solutions of Einstein’s field equations. The work below has been based on Hobson’s book ^[2], chapter $9$, Sections $1$ and $2$.

Because we are dealing with a spherically symmetric geometry, we use the spherical coordinate line metric form of Minkowski space, derived in the previous chapter. It now also consists of multiplying each of the terms by arbitrary coefficient functions of $r$ and $t$ to make this line metric unique to this Schwarzschild geometry, which are our metric functions, $g_{\alpha\beta}$. So the line metric becomes $$ds^2 = A(r,t)dt^2-B(r,t)dr^2-C(r,t)r^2(d\theta^2+\sin^2 (\theta) d\phi^2) - D(r,t)dtdr.\ \label{2.1}$$ The spherically symmetric geometry means that the only spatial coordinate that is affected by the gravitational field force is the radial coordinate, $r$, so by taking a new coordinate system of $r^2$, it is found that the function $C(r,t) =1$. The line metric needs to be invariant under changing the time component from positive to negative, $t \to -t$, because the geometry is at rest. This causes the coefficient multiplied to $drdt$ to equate to zero, hence we have $D(r,t) = 0$. This also requires our $g_{\alpha\beta}$ functions to be independent of the time-like coordinate, $t$, turning the line metric into

$$ds^2 = A(r)dt^2-B(r)dr^2-r^2(d\theta^2+\sin^2 (\theta) d\phi^2) .\ \label{2.2}$$ The line metric now satisfies the form of a spacetime metric that is static spatially isotropic because all of its terms are invariant under rotation of its space-like coordinates and all of its terms, except for the time-like coordinate, are independent of time.

To find expressions of the functions $A(r) \text{ and } B(r)$ a unique solution of the Einstein’s equations needs to be determined. Starting with the Ricci tensor, we find that it vanishes in empty space because of the absence of matter, which only contains a gravitational field, so we have the Ricci tensor as $$\begin{aligned} R_{ab} &= 0.\\ \end{aligned} \label{2.2a}$$ Equating the Riemann curvature equation with the Ricci tensor equation we obtain $$\begin{aligned} R_{ab} &\equiv \partial_b \Gamma^{c}_{ac} -\partial_c \Gamma^{c}_{ab} + \Gamma^{e}_{ac} \Gamma^{c}_{eb} - \Gamma^{e}_{ab} \Gamma^{c}_{ec},\\ \end{aligned} \label{2.2b}$$ with our Christoffel symbol, $\Gamma^{c}_{ab}$, as before having a relation with the metric, $g_{ab}$, equation as $$\begin{aligned} \Gamma^{\eta}_{\beta\zeta} &= \frac{1}{2} g^{\eta\alpha} (\partial_{\zeta}g_{\alpha\beta}+ \partial_{\beta}g_{\alpha\zeta}- \partial_{\alpha}g_{\beta\zeta}).\\ \end{aligned} \label{2.2c}$$ The calculus tensor method of solving the Einstein equations is done by directly calculating $\Gamma^{c}_{ab}$ from our equation, using the values of $g^{ab}$, where we have our expressions of $g_{ab}$ and $g^{ab}$ taken directly from the line metric equation to give $$\begin{aligned} g_{00} &= A(r), & g_{11} &= -B(r) , & g_{22} &= -r^2, &g_{33} &= -r^2\sin^2\theta, \\ g^{00} &=1/A(r), &g^{11} &= -1/B(r) ,&g^{22} &=-1/r^2, &g^{33} &= -1/r^2\sin^2\theta.\\ \end{aligned} \label{2.2d}$$ Using these values of $g^{\alpha\beta}$ and substituting them into equation gives the first non-zero Christoffel symbol, $\Gamma^{0}_{01}$, as $$\begin{aligned} \Gamma^{0}_{01} &= \frac{1}{2} g^{0\alpha} (\partial_{1}g_{\alpha0}+ \partial_{0}g_{\alpha1}- \partial_{\alpha}g_{01})=\frac{1}{2} g^{00} (\partial_{1}g_{00})= \frac{1}{2A(r)} \partial_1A(r) =\frac{A'}{2A}.\\ \end{aligned} \label{2.2da}$$ Repeating this for all of the Christoffel symbols, we find that there are only 9 non-zero elements of the 40 independent connection coefficients, with them being $$\begin{aligned} \Gamma^{0}_{01} &= A'/(2A), &\Gamma^{1}_{00} &= A'/(2B), & \Gamma^{1}_{11} &= B'/(2B),\\ \Gamma^{1}_{22} &= -r/(B), &\Gamma^{1}_{33} &= -r^2 \sin^2 \theta/(B),& \Gamma^{2}_{12} &= 1/r,\\ \Gamma^{2}_{33} &= -\sin\theta\cos\theta,& \Gamma^{3}_{13} &= 1/r, & \Gamma^{3}_{23} &= \cot\theta.\\ \end{aligned} \label{2.2e}$$ The Ricci tensors are then calculated tediously by substituting these Christoffel symbols into equation. Substituting in the relevant expressions gives the first Ricci tensor as $$\begin{aligned} R_{00} &= \partial_0 \Gamma^{\alpha}_{0\alpha} -\partial_\alpha \Gamma^{\alpha}_{00} + \Gamma^{\beta}_{0\alpha} \Gamma^{\alpha}_{\beta0} - \Gamma^{\beta}_{00} \Gamma^{\alpha}_{\beta\alpha} ,\\ &= -\partial_1 \Gamma^{1}_{00} + \Gamma^{0}_{01} \Gamma^{1}_{00} - \Gamma^{1}_{00}( \Gamma^{1}_{11}+ \Gamma^{2}_{12} + \Gamma^{3}_{13}),\\ & = \frac{-2A''B + A'B'}{4B^2} +\frac{A'}{2B'} +\frac{A'}{2A}\frac{A'}{2B} - \frac{A'}{2B}\left( \frac{B'}{2B}+ \frac{1}{r}+ \frac{1}{r}\right),\\ &= -\frac{A''}{2B} +\frac{A'}{4B} \left(\frac{A'}{A} + \frac{B'}{B}\right) - \frac{A'}{rB}.\\ \end{aligned} \label{2.2ea}$$ All of the non-diagonal components of the Ricci tensors equations equate to zero, $R_{\alpha\beta}=0$ when $\alpha \neq \beta$, giving the diagonal terms as $$\begin{aligned} R_{00} &= -\frac{A''}{2B} +\frac{A'}{4B} \left(\frac{A'}{A} + \frac{B'}{B}\right) - \frac{A'}{rB},\\ R_{11} &= \frac{A''}{2A} -\frac{A'}{4A} \left(\frac{A'}{A} + \frac{B'}{B}\right) - \frac{B'}{rB},\\ R_{22} &= \frac{1}{B} -1 +\frac{r}{2B} \left(\frac{A'}{A} - \frac{B'}{B}\right),\\ R_{33} &= R_{22} \sin^2\theta.\\ \end{aligned} \label{2.2f}$$ Equating the Ricci tensor equations to zero as they are in an empty space and simplifying the terms, by adding $\frac{B}{A} (R_{00}) + R_{11}$ gives the following relation $$\begin{aligned} A'B + AB' &= 0,\\ \end{aligned} \label{2.2g}$$ that proves that $AB$ is constant. So we then set $AB= \alpha$ where $\alpha$ is an arbitrary constant. Substituting $B=\frac{\alpha}{A}$ into our expression for $R_{22}$ in the equation above gives $$\begin{aligned} \frac{A}{\alpha} -1 +\frac{rA}{2\alpha} \left(\frac{A'}{A} - \frac{A}{A'}\right) &= \frac{A}{\alpha} -1 +\frac{rA}{2\alpha} \left(\frac{2A'}{A} \right) =0,\\ \alpha &= A+rA',\\ \Rightarrow \frac{d(rA)}{dr} &= \alpha.\\ \end{aligned} \label{2.2h}$$ By integrating the last expression in the equation above we then get $$\begin{aligned} rA&=\alpha(r+k),\\ \end{aligned} \label{2.2ha}$$ where $k$ is a constant of integration. So we finally arrive at equations of the functions $A$ and $B$ as $$\begin{aligned} A(r) &= \alpha \left( 1 + \frac{k}{r}\right) ,\\ B(r) &= \left( 1 + \frac{k}{r}\right)^{-1}.\\ \end{aligned} \label{2.2i}$$

Comparing equations above to that of equations of a weak-field, we see that the constant $k$ needs to represents the mass of an object that is causing the gravitational field, so we substitute $k$ with the Newtonian gravitational potential, $\Phi$, so that $k=2\Phi$ and then we require

$$\begin{aligned} \frac{A(r)}{c^2} &\to 1 + \frac{2\Phi}{c^2}.\\ \end{aligned} \label{2.2j}$$ For a spherically symmetric mass, $M$, and radial distance, $r$, from the centre of the mass we have $\Phi = -GM/r$, hence the constant $k = -2GM/c^2$ and constant $\alpha= c^2$. Setting the usual constants to one with $G=c=1$ the functions of $r$ then become

$$\begin{aligned} A(r) &= 1-\frac{2GM}{c^2 r} = 1-\frac{2M}{ r},\\ B(r) &= c^2 A(r)^{-1} = \frac{c^2}{1-\frac{2GM}{c^2 r}} = \frac{1}{1-\frac{2M}{ r}}.\\ \end{aligned} \label{2.3}$$ The geometry very far away from the black hole needs to tend towards the Minkowski line geometry, because at infinite distances away, the gravitational force of the black hole will tend to zero, hence as $r \to \infty$ then $ds^2 \to \eta_{\alpha\beta}(x) dx^\alpha dx^\beta$. This forces the functions $A(r) \text{ and } B(r)$ to tend to $1$, with $\lim_{r \to \infty} A(r) = c^2 =1$ and $\lim_{r \to \infty} B(r) = 1$. Our expressions for $A(r) \text{ and } B(r)$ satisfy these properties and hence agree with the Minkowski spacetime arguments at infinite distances from the black hole.

Substituting these functions of $r$ into our line metric equation yields the completed Schwarzschild line metric for the empty spacetime outside the spherically symmetric body as

$$ds^2 = \left(1-\frac{2M}{ r}\right)dt^2-\frac{1}{1-\frac{2M}{ r}}dr^2- r^2(d\theta^2+\sin^2 (\theta) d\phi^2),\ \label{2.5}$$ with the metric tensor being $$g_{\alpha,\beta} = \begin{pmatrix} (1-\frac{2M}{ r}) & 0 & 0 & 0 \\ 0 & -(1-\frac{2M}{ r})^{-1} & 0 & 0 \\ 0 & 0 & -r^2& 0 \\ 0 & 0 & 0 & -r^2\sin^2\theta\\ \end{pmatrix}.$$

Now that the line metric has been calculated, a system of first order differential equations can be derived that serve as the equations of motion for photons along null geodesics in the Schwarzschild geometry.

3.2 The Equations of Motion

To derive the equations of motion, the line metric needs to be simplified. To do this one of the spatial dimensions of the geometry is removed, forcing the photon paths to be ray traced in a two dimensional equatorial plane through the black hole with the three coordinates, $(t, r, \phi)$. Hence the $\theta$ component is considered as being removed from our Schwarzschild line metric equation , so that our model now traces light in the equatorial plane only. The spherical symmetry in this geometry allows us to remove $\theta$ without loss of generality. Work in this section has been based on Shultz’s book ^[1], Chapter 11.1.

By removing the $\theta$ component from the Schwarzschild line metric equation and equating $ds^2=0$, because photons do not feel real time and hence they do not travel any real distance, gives

$$0 = -A(r) \dot{t}^2 + \frac{1}{A(r)} \dot{r}^2 + r^2 \dot{\phi}^2, \\$$ where $A(r)$ is our function defined as $A(r) = (1-2M/ r)$. The origin of this system is at the centre of the black hole and the photons have polar coordinates, $\{r, \phi\}$.

The energy of photons in this geometry, along each trajectory, are constant and can be set equal to $E=-\rho_0$, since time is independent in this metric. Spherical symmetry in this metric causes the angular momentum of photons along trajectories also to be constant and set equal to $L= p_\phi$. Because of these two conserved quantities it makes the momentum components of the photon become $$\begin{aligned} p^0 &= \left(1-\frac{2M}{r}\right) ^{-1} E ,\\ p^r &= \frac{dr}{d\lambda} = \dot{r}, \\ p^{\phi} &= \frac{d\phi}{d\lambda} = \dot\phi= \frac{L}{r^2}. \\ \end{aligned} \label{2.6}$$ The affine parameter, $\lambda$, in the above equation is defined by the photon’s $p^r$. Substituting these momentums into the line metric of the photon and setting $ds^2 =0$ produces $$\begin{aligned} 0 &= -E^2 \left(1-\frac{2M}{r}\right) ^{-1} + \left(1-\frac{2M}{r}\right) ^{-1}\left(\frac{dr}{d\lambda}\right)^2 + \frac{L^2}{r^2} .\\ \end{aligned} \label{2.7}$$ Rearranging the equation results in $$\begin{aligned} \left(\frac{dr}{d\lambda}\right)^2 &= E^2 - \left(1-\frac{2M}{r}\right) \frac{L^2}{r^2} = E^2 -V_{eff} \\ \end{aligned} \label{2.8}$$ where $V_{eff}$ is defined as the effective potential of the photon. This potential will be discussed in more detail in the next chapter. Differentiating this equation with respect to $\lambda$ yields

$$\begin{aligned} 2\left(\frac{dr}{d\lambda}\right) \left(\frac{d^2r}{d\lambda^2}\right) &= - \frac{d}{dr}V_{eff}(r)\left(\frac{dr}{d\lambda}\right) . \\ \end{aligned} \label{2.8b}$$

with

$$\begin{aligned} \left(\frac{d^2r}{d\lambda^2}\right) &= - \frac{1}{2}\frac{d}{dr}V_{eff}(r). \\ \end{aligned} \label{2.8c}$$

This can then be expressed as

$$\begin{aligned} \ddot{r} &= -\frac{L^2}{2} \left(-\frac{2}{r^3} + \frac{6M}{r^4}\right), \\ \ddot{r} &= -\frac{L^2(3M-r)}{r^4}.\\ \end{aligned} \label{2.9}$$

To make this second order equation into a system of first order equations, a variable $s$ is defined as $s\equiv \dot{r}$. Substituting this into the above equation then gives the following system of first order equations:

$$\begin{aligned} \dot{\phi} &= L/r^2\,,\\[1ex] \dot{r} &= s\,,\\[1ex] \dot{s} &= -\frac{L^2(3M-r)}{r^4}\,. \label{2.10} \end{aligned}$$ The equations above are those that govern the motion of a photon in the equatorial plane of the Schwarzschild geometry. They can be analytically solved without much difficulty, but the equations of motion become a lot more complicated for less symmetric cases of black holes, such as rotating ones, which won’t be covered as they are beyond the scope of this report. For rotating black holes, the geometries are no longer symmetric, hence they cannot be generalised to two dimensions.

Now that these equations of motion have been derived we can use them to ray trace the photon’s paths around Schwarzschild black holes.

3.3 Ray Tracing in the Schwarzschild Geometry

To produce images of Schwarzschild black holes our code will need to ray trace the photons paths. The images will be what appears from a camera that is looking directly at a black hole with a celestial sphere star field at an infinite distance behind it. This is shown in Figure 3.1 where a pinhole camera is on the left, the black hole is the black disc in the centre and a star field on a celestial sphere is on the right. In this figure we see a photon path, the green line, starting from a point on the celestial sphere being diverted around the black hole and still ending in the pinhole camera. So from the pinhole camera’s perspective the image of the point on the celestial sphere looks like it is coming from a different direction than it actually is.

Pinhole cameras work by having an infinitesimally small hole in a light-proof box as shown in Figure 3.2. When light from various angles pass through this aperture it forms an image that has been inverted both vertically and horizontally on a photographic plate on the opposite side of the box. The width of the lightproof box is known as it focal length, $f$, which is a constant we set before running the code. The focal length affects how much light passes into the box, to determine the size and view of the image. The code produces the images by modelling what appears on the photographic plate of a pinhole camera in Figure 3.2.

To ray trace these paths an integrating function in python solves the three first order differential equations for the Schwarzschild Geometry, by starting the photons at the aperture of the pinhole camera and iterates over each time step, by integrating the equations to produce an array, ($\phi,r,s$), of the photon’s location at each time step. This function uses the initial values of ($\phi,r,s$) are used to find solutions of the derivatives for each time step.It then integrates over the functions to produce the array of the photon’s locations as $z= [[\phi_0,r_0,s_0],…,[\phi_t,r_t,s_t]]$ for each time step from $0$ to $t$. Setting the amount of time steps to a very large integer, gives an accurate location of the photon’s location when $t \to \infty$. Plotting the array of locations of the photon shows the paths that the light travels along in two dimensions around and near the black hole. Using a range of values of angular momentum, $L$, for the photons we see the different paths that light travels along.

t  = np.linspace(0,1000,200) # time steps
yo = [2, 3]                  # constants [M, L]
zo = [np.pi, 20, -0.13]      # initial [phi0, r0, s0]

def function(z, t, y):      # ODE functions
    phi, r, s = z
    M, L = y

    dpdt = L/r**2
    drdt = s
    dsdt = -(L**2)*(3.0*M-r)/(r**4)
    return [dpdt, drdt, dsdt]

z = scipy.integrate.odeint(function,zo,t, args=(yo,))  # solve ODE
phi, r, s = z[:,0], z[:,1], z[:,2]

Ray tracing these light paths for each pixel on the photographic plate will take the code large periods of time to run and produce good quality images, so it would be an ineffective method. To speed up the code, we find a relation between the curvature of the light paths, which generalises for all paths near the black hole. Because this is a symmetric problem we can find a correlation between two different angles in the two dimensional plane. Figure 3.3 shows the first angle being $\alpha'$, which is the angle under which a star’s light would classically enter the camera’s aperture in the absence of the black hole, which is the star’s natural angle as it would be in Newtonian physics. The second angle, $\alpha$, in Figure 3.3 is the angle at which the light from the same star actually enters the camera’s aperture due to the bending caused by the black hole’s gravitational field. Both of these angles will be measured with respect to the horizontal in the black hole’s equatorial plane.

These angles are calculated simply by comparing the different locations of photons and using their change in positions. To get the real angle of the star, $\alpha'$, the end location of the photon is recorded after many time iterations to get the angle $\phi$ which is what $\alpha'$ tends to when you hit the celestial sphere because it is at an infinite distance away. The angle $\alpha$ is calculated using trigonometric equations using the change in location of the photon between the first two time steps. After getting a list of various angles of $\alpha$ and $\alpha'$ by using many different values of angular momentum, $L$, an interpolating function in python produces a map from $\alpha \to \alpha'$ by comparing the relations in values of the two angles. This map acts as a function that after an angle $\alpha$ is entered into it, it finds the corresponding value of $\alpha'$.

This map between $\alpha$ and $\alpha'$ has been plotted, shown in Figure 3.4, where as you would expect for small angles of $\alpha$ the plot has no corresponding $\alpha'$ values as the photons have been absorbed by the black hole and cannot escape. For this reason there do not exist any $\alpha$ values lower than this one, as there are limitations to the amount of gravitational lensing that can occur. This is when light travels around the black hole in different directions and intersect in the future, which will be explained in more detail in the next chapter. The photon that has the maximum curvature around the black hole is plotted with the smallest value of $\alpha$, hence no plots exists below this. Increasing $\alpha$ from this point the plot is steep with a large gradient and corresponds to negative $\alpha'$ values, which has been caused by photons traveling from the lower side of the celestial sphere, curving up and around the top of black hole and still making it to the pinhole camera. As $\alpha$ gets larger, and closer to $\pi/2$, the graph becomes less steep and the values of $\alpha$ and $\alpha'$ tend to $\alpha=\alpha'$, represented by the straight blue dashed line. This occurs because when the photons are traveling further away from the black hole, they are less curved by it, causing the map of the angles to have a direct correlation. For negative $\alpha$ values the plot would appear as the plot in Figure 3.4 reflected in the $\alpha = \alpha'$ line.

Another angle, $\beta$, is introduced to generalise, by rotational symmetry, our model to a three dimensional problem as $\theta$ was removed from the line metric equation to simplify the ray tracing in the two dimensional equatorial plane. The angle $\beta$ will be defined as the angle between the photon and the direct vertical, which we were using as the equatorial plane. In this model light paths will only be traced in two dimensions, so that the angle $\beta$ will not be affected by the curvature of spacetime. The only change to $\beta$ will be from the pinhole camera. So after the photon has left the pinhole camera, the angle is just inverted horizontally and vertically as we pass the black hole because the Schwarzschild geometry is symmetric and because of the properties of pinhole cameras. So this relation is defined as $\beta' =- \beta$.

Now that all of the crucial steps of ray tracing have been covered, a backwards ray tracing algorithm can now be derived to produce the black hole images.

3.4 The Backwards Ray Tracing Algorithm

To produce these images of Schwarzschild black holes, a code with an algorithm that will be outlined below is used. The following steps need to be taken for the code to produce an image of a black hole:

• Set the initial location of the pinhole camera’s aperture, $\{ r, s,\phi\}$ at $\tau =0$, the focal length, $f$, and the mass of the black hole, $M$.

• Start at a pixel, $\{i,j\}$, on the photographic plate of the pinhole camera and calculate the angles $\alpha$ and $\beta$ of this photon. This is done by using the change in distances of this pixel from the centre of the photographic plate as shown in the left image of Figure 3.5. It is showing the photographic plate of width $x$ and height $y$ in pixels, the change in vertical and horizontal distance, $dx,dy$, and the hypotenuse of the triangle being $dxy$. These equations are defined as follows: $$\begin{aligned} (dx,dy) &= \left(\frac{x}{2} -i,\frac{y}{2} -j\right), \\ dxy &= \sqrt{dx^2+dy^2}.\\ \end{aligned} \label{2.11a}$$

  

  The angle $\alpha$ is calculated using the focal length, $f$, and the distance of the photon from the centre of the photographic plate as

  $$\begin{aligned} \alpha &= \tan^{-1}\left(\frac{dxy}{f}\right).\\ \end{aligned} \label{2.11}$$ The angle $\beta$ is calculated by using the position of the pixel compared to the vertical, which changes depending on which side of the photographic plate it is on. This angle is defined by $$\beta = \begin{cases} 0, & \text{if}\ dxy=0, \\ \cos^{-1}(dy/dxy), & \text{if} \ i\geq x/2,\\ -\cos^{-1}(dy/dxy), & \text{if} \ i<x/2.\\ \end{cases} \label{2.12}$$

• The interpolating function discussed in the previous section is used to map $\alpha \to \alpha'$, . This gives $\alpha' = \alpha$ when $\alpha$ is larger than the greatest $\alpha$ value in the plot in Figure 3.4 because the plot tends towards the blue dashed line. For the rest of the $\alpha$ values, the interpolating map is used to evaluate $\alpha'$. The angle $\beta$ is just inverted horizontally and vertically as we pass the black hole so therefore we get $\beta' = - \beta$.

• Use $\alpha' \text{ and } \beta'$ to calculate the location of where the photon would have originally come from on the celestial sphere. Calculate this pixel on the celestial sphere, $\{i',j'\}$, by first calculating the change in distances on the celestial sphere image on the right image of Figure 3.5, with $$\begin{aligned} dxy' &= f (\tan(\alpha')),\\ (dx',dy') &= (dxy' \sin(\beta'),dxy' \cos(\beta')), \\ (i',j') &=\left(\frac{x}{2} + dx', \frac{y}{2} + dy'\right).\\ \end{aligned} \label{2.13}$$

• The red, green and blue colour scheme is used where each pixel on the celestial sphere is assigned three integer values between $0$ and $255$ depending on its colour. Then, the colour of the pixel, $\{i',j'\}$, is extracted from the celestial sphere, so that the corresponding pixel, $\{i,j\}$, on the photographic plate takes this colour.

• If the photon from a pixel does not make it to the celestial sphere and is instead absorbed by the black hole you set the colour of that pixel as (0,0,0), which represents the colour black, as no light has been absorbed from this direction. This is if the angle of $\alpha$ is lower than that of the lowest point plotted in our Figure 3.4.

• Repeating this for all pixels, $\{i,j\}$ for $i\in \{x\}$, $j\in \{y\}$, on the photographic plate produces a deformed image of a Schwarzschild black hole.

Now that the photon’s equations of motions in the Schwarzschild geometry have been derived and used in our piece of code that contains the algorithm to ray trace and produce images of these black holes, we can move onto analysing the images of the black holes.