We start from the Schwarzschild metric(Lambourne 2010) for a simple black hole with zero angular momentum and zero charge,\[d s^{2} = \left( 1 - \frac{r_{s}}{r} \right) d t^{2} - \left( 1 - \frac{r_{s}}{r}\right)^{-1} d r^{2} - r^{2}d\Omega^{2},\] using the standard convention \(G = c = 1\). Taking advantage of the spherical symmetry of the problem we can always choose a coordinate reference in the orbital plane \(\left(\theta = \frac{\pi}{2}\right)\) of our test particle. Doing so simplifies our problem greatly, transforming our line element into, \[d s^{2} = \left( 1 - \frac{r_{s}}{r} \right) d t^{2} - \left( 1 - \frac{r_{s}}{r}\right)^{-1} d r^{2} - r^{2} d \phi^{2}.\] It is useful at this time to define the constants \(l\) and \(e\), derived directly from our metric \(g_{\mu\nu}\), which represent the pseudo angular momentum and pseudo energy respectively. These come from the *Killing Vectors*(Antonelli 2015) and are represented as, \[r^{2}\frac{d\phi}{d\tau} = l,\] \[\left(1 - \frac{r_{s}}{r}\right)\frac{dt}{d\tau} = e.\] By substituting (3) and (4) into (2) and dividing through by \(d\phi^{2}\) it can be shown through minimal algebra that (2) can be written of the form, \[\left(\frac{d r}{d\phi}\right)^{2} = \frac{r^{4}}{b^{2}} - \left( 1 - \frac{r_{s}}{r}\right)\left[\frac{r^{4}}{a^{2}} + r^{2}\right].\] The constants \(a = l\) and \(b = \frac{l}{e}\) are defined in order to give this equation the familiar form of an elliptical orbit. We make one more substitution of \(u = \frac{1}{r}\), with \(d r = - r^{2} d u\) to arrive at the convenient form of (5), \[\left(\frac{d u}{d\phi}\right)^{2} = \frac{1}{b^{2}} - \left(1 - r_{s}u\right)\left[a^{-2} + u^{2}\right].\] By taking the limit as \(a\rightarrow\infty\) to represent the massless particle we can derive an equation of motion by taking the derivative of (6) with respect to \(\phi\). For convenience we will represent \(\frac{d u}{d\phi}\) as \(\dot{u}\) with the dot operator being in respect to \(\frac{d}{d\phi}\). \[\begin{gathered}
\frac{d}{d\phi}\left(\dot{u}^{2} = \frac{1}{b^{2}} - \left( 1 - r_{s}u \right)u^{2}\right),\\
2\dot{u}\ddot{u} = - 2 \dot{u} u + 3 r_{s} \dot{u} u^{2},\\\end{gathered}\] \[\ddot{u}\left(\phi\right) = -u\left(1-\frac{3r_{s}}{2}u\right).\]

If we look carefully at (7) we can recognize that for motion in the orbital plane this equation exactly satisfies the Binet Equation(Antonelli 2015) for a massive particle moving in a central potential, \[\begin{gathered} \frac{d^{2}}{d t^{2}}\vec{x} = \frac{1}{m}F\left(r\right),\\\end{gathered}\] \[\ddot{u} + u = \frac{-1}{m h^{2} u^{2}}F(u).\] We can finally plug equation (7) into (8) and resubstitute \(r=\frac{1}{u}\) to obtain our effective central force on a massless particle in orbit near a black hole, \[\vec{F}(r) = - \frac{3r_{s}}{2}h^{2}\frac{\hat{r}}{r^{4}}.\] We have set \(m=1\) for simplicity. To arrive at (9) we make the step that by evolving a particle under a central force we arrive at the Binet equation, then the Binet equation may be used as a tool to solve for our central force(Antonelli 2015).

To solve our path of light due to the effective central force in (9) we use a Runge-Kutta 4 algorithm(Press 2007) to step the equation forward. We use the standard definition for a force that depends on a single variable, \[\begin{gathered} \dot{x} = v, \\ \ddot{x} = \vec{F}(x). \\\end{gathered}\]

The RK4 algorithm can then easily be summarized as the set of four smaller steps combined to perform a single step of size \(h\) as follows,

\[\begin{gathered} k^{0}_{1} = \dot{x_{n}}, && k^{1}_{1} = \vec{F}\left(x_{n}\right), \\ k^{0}_{2} = \dot{x_{n}} + \frac{h}{2}k^{1}_{1}, && k^{1}_{2} = \vec{F}\left(x_{n} + \frac{h}{2}k^{0}_{1}\right), \\ k^{0}_{3} = \dot{x_{n}} + \frac{h}{2}k^{1}_{2}, && k^{1}_{3} = \vec{F}\left(x_{n} + \frac{h}{2}k^{0}_{2}\right), \\ k^{0}_{4} = \dot{x_{n}} + h k^{1}_{3}, && k^{1}_{4} = \vec{F}\left(x_{n} + h k^{0}_{3}\right), \\ x_{n+1} = x_{n} + \frac{h}{6}\left(k^{0}_{1} + 2 k^{0}_{2} + 2 k^{0}_{3} + k^{0}_{4}\right), && \dot{x}_{n+1} = \dot{x}_{n} + \frac{h}{6}\left(k^{1}_{1} + 2 k^{1}_{2} + 2 k^{1}_{3} + k^{1}_{4}\right). \\\end{gathered}\]

For \(l\) light rays we evolve an initial condition of \(\left(x_{\small{0}},v_{\small{0}}\right)\) along a total path length of \(s\). Our step size \(h\) is given by \(\small{h} = \frac{s}{n}\) where n is the total number of points along our path or, our resolution when calculating the RK4 algorithm. A listing of the implementation is presented in appendix A.

We have two light sources used to generate initial conditions for \(l\) rays. The first is a set of parallel rays coming in from \(r_{\small{0}} \small{\gg} r_{s}\). These have an initial impact parameter, \(b\), and a distribution of \(w\) with equal spacing of all rays, in units of \(r_{s}\). The second is an angular distribution of light rays all originating from the same point \(r_{\small{0}}\). Their angular distribution is again represented by \(w\), with \(\small{0} \small{\leq} w \small{\leq} \small{2\pi}\). All \(l\) rays are evenly distributed through this angle \(w\). Implementation of these sources are presented in appendix B.

```
// Vector2_2 is of the form (pos,vel) -> ((x,y),(xdot,ydot))
private Vector2_2 RK4F(Vector2_2 r, float h2)
{
Vector2_2 R = new Vector2_2(r.Vel, Vector2.zero);
Vector2 v = r.Pos;
// r^-4 potential that approximates path of light near black hole
v = - (1.5f * h2 / Mathf.Pow(v.magnitude, 4f)) * v.normalized;
R.Vel = v;
return R;
}
private Vector2_2 RK4(Vector2_2 r, float step, float h2)
{
// Standard Runge-Kutta 4
Vector2_2 k1 = RK4F(r, h2);
Vector2_2 k2 = RK4F(r + 0.5f * step * k1, h2);
Vector2_2 k3 = RK4F(r + 0.5f * step * k2, h2);
Vector2_2 k4 = RK4F(r + step * k3, h2);
Vector2_2 k = (float)(step / 6f) * (k1 + 2f * k2 + 2f * k3 + k4);
return r + k;
}
private void DrawLightPath(Vector2_2 r, int k)
{
// step size is total path length (s) divided by total steps (n)
float step = (float)(s / n);
// square of angular momentum in the orbital plane (r x v)
float h2 = Mathf.Pow( r.Pos.x * r.Vel.y - r.Pos.y * r.Vel.x, 2f );
for (int i = 0; i < n; i++)
{
// Runge-Kutta 4 to increment our light point r a single step
r = RK4(r, step, h2);
// Process each intermediate value of r here...
}
// Process final value of r here...
}
private struct Vector2_2 {
Vector2 pos;
Vector2 vel;
public Vector2_2(Vector2 p, Vector2 v)
{
pos = p;
vel = v;
}
public Vector2 Pos {
get{ return pos; }
set{ pos = value; }
}
public Vector2 Vel {
get{ return vel; }
set{ vel = value; }
}
public static Vector2_2 operator *(float f, Vector2_2 r)
{
return new Vector2_2(f * r.Pos, f * r.Vel);
}
public static Vector2_2 operator +(Vector2_2 r1, Vector2_2 r2)
{
return new Vector2_2(r1.Pos + r2.Pos, r1.Vel + r2.Vel);
}
public static Vector2_2 operator -(Vector2_2 r1, Vector2_2 r2)
{
return new Vector2_2(r1.Pos - r2.Pos, r1.Vel - r2.Vel);
}
};
```

```
/*
i - ray number
b - impact parameter (z position)
l - number of lines
w - distribution/spread
x - x position
*/
private static Vector2_2 ParallelLightSource(int i, int n, float b, int s, int l, float w, float x)
{
return new Vector2_2(new Vector2(x, b + (float)i * w / (float)l), Vector2.left);
}
private static Vector2_2 ConeLightSource(int i, int n, float b, int s, int l, float w, float x)
{
Vector3 temp = new Vector3( -1f, 0f, 0f);
Vector3 target = new Vector3(1f, 0f, 0f);
temp = Vector3.RotateTowards(temp,target, (i - (float)(l-1f)/2f) * w / (float)l, 0f);
return new Vector2_2(new Vector2(x, b), new Vector2(temp.x, temp.z));
}
```

## Share on Social Media