# Light Bending 2D

## Light Path Derivation

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).

## Light Path Computation

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.

## Sources

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.

# Appendicies

## Appendex A: Listings

// Vector2_2 is of the form (pos,vel) -&gt; ((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 &lt; 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);
}
}; 

## Appendex B: Ray Sources

/*
i - ray number
b - impact parameter (z position)
l - number of lines
}