Thoughts about light bending

ebvalaim 2015-06-17 0 Comments

I had a sudden moment of clarity recently when thinking about how to remove a graphical artifact from the Black Hole Simulator.

The problem

The current simulator has one ugly aspect of the rendered image. Exactly 90 degrees from the direction to the black hole a graphical artifact appears - a strip of smudged, incorrectly calculated pixels (see the picture to the right). The reason is hidden deeply in the rendering mechanism.

In short, it looks like this - you can't efficiently calculate the color of every pixel by raytracing. Taking advantage of the symmetry of the Schwarzschild black hole, I created a table of angles of light deflection. Thanks to the symmetry, I can describe each ray of light by only one parameter - simply put, the minimal distance from the black hole along its path (actually, the impact parameter). To calculate the deflection, I need one more thing, and that is the distance from the black hole at which I send/receive it - the greater part of the whole path the ray has to travel, the greater the deflection.

Such a table was being sent to the graphics card. Then, during rendering, a direction of the ray was calculated for each pixel, and this was converted into the impact parameter. The distance was known independently. Appropriate deflection was being read from the table and this was used to calculate the color of the pixel.

In theory everything is fine, but one problem appeared - the light rays sent in directions close to 90 degrees from the black hole have very similar impact parameters. This gives nearly identical deflection angles, which makes many pixels the same color. An ugly strip appears.

The solution

The whole problem arises because of the usage of the impact parameter for description of the light ray. This makes it easier to generate the deflection table, but as you can see, it causes big problems. What if we abandoned this idea? Let's try to describe each ray with the distance from the black hole and the angle between its direction and the direction to the black hole.

First we need a way of generating the four-velocity of a ray being given an angle. It would be done most convieniently using two unit, orthogonal spatial vectors (with which we would generate a circle) and one timelike vector, which we would add with an appropriate coefficient, such that we will get a null vector (null vectors describe light rays).

In order to be able to describe the interior of the black hole as well as its exterior, the simulator uses the Eddington-Finkelstein coordinates $(u, r, \theta, \phi)$ (I think in the literature what I denote by $u$ is denoted by $v$ instead - I'll leave it as $u$, because I use this notation everywhere and a change would cause mistakes). In these coordinates $\partial_u$ is the same as $\partial_t$ in Schwarzschild coordinates and $\partial_r$ is a null vector directed to the past.

$\partial_r$ being directed to the past is actually a good thing - the observer doesn't send light rays, but receives them. The calculations are being started at the observer though, so we need to follow the paths of the rays into the past. The initial direction should thus point to the past.

The $\partial_r$ vector is then a great candidate for one of our parametrized rays - I'll want to get it as a result for the angle $\alpha=\pi$ ($\alpha$ will be the angle from the direction towards the black hole). Let's now find a good ray for the direction $\alpha=0$.

For now we will only consider the radial direction and time, which means that we don't care for $\theta$ nor $\phi$. This leaves us with $v = a\partial_r + b\partial_u$.
It would be nice for this vector to be "normalized": $g(v, \partial_r) = 1$.
And, of course, it should be a null vector: $g(v, v) = 0$.

In E-F coordinates we have $g_{uu} = \left(1-\frac{2M}{r}\right)$, $g_{ur} = -1$, $g_{rr} = 0$. This yields:
$1 = -b$
$0 = -2ab + b^2\left(1-\frac{2M}{r}\right)$

The solution:
$b = -1$
$a = \frac{M}{r} - \frac{1}{2}$

From this we get $v = \left(\frac{M}{r} - \frac{1}{2}\right) \partial_r - \partial_u$.

Having two null vectors, one pointing towards and one away from the black hole, it is easy to find a timelike vector - we just add them together $v+\partial_r$ and normalize the result. This yields:
$T = \frac{1}{\sqrt{2}}\left[ \left( \frac{M}{r} + \frac{1}{2} \right) \partial_r - \partial_u \right]$

Now, one of the spatial vectors will be a combination of $\partial_r$ and $T$, hence we get:
$A = \sqrt{2}\partial_r - T$
The second one will be just in the direction of $\partial_\phi$:
$B = \frac{1}{r}\partial_\phi$

We now have a complete orthonormal basis for the rays. A ray sent in the direction at an angle $\alpha$ to the direction towards the black hole will have its four-velocity equal to:
$U(\alpha) = T - A\cos\alpha + B\sin\alpha$
($A$ is negated so that $\alpha=0$ means the direction to the black hole, not away from it).

Because our spatial vectors are orthogonal to $T$, $-T$ is the four-velocity of an observer relative to which we generate the circle. As it turns out, this observer moves relative to the black hole, which means that $\alpha$ can't be directly converted into a pixel of the background, but fortunately it's not a problem. $\alpha$ will only be needed as a parametrization of the table, but the final point reached by the ray will be defined by the $\phi$ coordinate of the ray's position, and it is uniformly distributed in the resting frame. Small problems may arise when picking appropriate $\alpha$ angles for the entries in the table, but I think I'll solve this in an entirely different way.

The pixel calculation algorithm will therefore look like this: calculate the direction of the ray for each pixel; calculate $\alpha$ from the direction (by $\cos\alpha = -g(U, A)$); read the deflection angle for given $\alpha$ and $r$; find the proper pixel.

What is left to decide is what the values in the table should actually be. If the simulator is to render the black hole also from great distances, we will need something that extrapolates nicely or can be calculated approximately when the observer is far from the black hole. Since possible distances are infinite, we won't be able to generate table values for all of them. A nice candidate is the difference from the final angle that we would get in a flat space.

The final $\phi$ in a flat space is easily found: we need to find $U(\alpha)$, substitute $M=0$, find the coordinates in the directions $\partial_t$ and $\partial_R$ (remembering that $\partial_t = \partial_u$ and $\partial_r = \partial_R - \partial_t$) in ordinary spherical coordinates and from that calculate $\cos \alpha'$. We get:
$\cos\alpha' = \frac{3\cos\alpha - 1}{3-\cos\alpha}$

Assuming that the observer is at $\phi=\pi$ and $\alpha=0$ is in the direction of the black hole, $\alpha'$ will be exactly the final $\phi$ of the ray with no deflections. We can subtract that from the raytraced coordinate and thus find the difference from the flat case, which should be approximately $\frac{4M}{r_{min}}$ for large distances ($r_{min}$ being the minimal distance from the black hole along the path of the ray, related to the impact parameter $b$ by $\frac{1}{b^2} = \left(1-\frac{2M}{r_{min}}\right)\frac{1}{r_{min}^2}$).

For further improvements, we can try to approximate the deflections by an analytic function of the angle for each distance. It is something I intend to explore a bit more. If it can be done, what will be left is putting the values of the coefficients defining the function in the table and then interpolating them between predefined distances. This would be probably the best way. One important thing to remember: for a range of the angles there won't be any values of the deflection. That's because the rays sent in some directions will hit the black hole. The critical impact parameter will be given by $r_{min} = 3M$ (the radius of the photon sphere), which gives $b = 3\sqrt{3}M$. Calculatig critical $\alpha$ from that will be relatively simple.

My plan for the next days is therefore this: try to generate the deflection angles for different values of $r$, draw charts and check if some reasonable function can be chosen to approximate them. Of course, I will describe the progress here :)

'Till the next time!