Intro

In flat spacetime, light moves on straight lines. In curved spacetime, light moves on geodesics, the straightest-possible lines as prescribed by the geometry. Around a black hole, spacetime is so curved that light can go into orbit! Below you can play around with those orbits and see what they look like.

You have control over two numbers: the spin $0\le a \le 1$, and the radius of the orbit (the limits depend on the spin).

You can also ask the code to try to match some frequency ratio (it might be impossible). Try some small-integer ratios, and drag the spin until a solution is found (the r slider will be disabled).

The details behind the math and physics appear below.

Drag/zoom around the 3d visualization to get a feeling for these orbits.

Details

Everything is in geometric units, $G=1=c$, and furthermore the total mass has been scaled out, $M=1$. The $a$ which appears here is $\tfrac{a}{M}$, and the $r$ is $\tfrac{r}{M}$. Everything is in Boyer-Lindquist coordinates.

Controls and visualization

In the visualization you will see a red curve animating around, a gray ellipsoid, and a blue arrow. The gray ellipsoid represents the event horizon (a constant-r surface in B-L coords), and the blue arrow represents the direction of the spin vector. The red curve is a chunk of the photon’s trajectory. Some history of the trajectory is saved.

• Play/pause: no explanation needed.
• $a$ slider: this controls the spin of the black hole, from $0 \le a \le 1$.
• $r$ slider: this controls the radius of the orbit. The minimum and maximum radius depend on the spin. Disabled when “Try to match” is checked.
• History length slider: how much history of the photon’s trajectory is plotted. This is in terms of an arbitrary affine parameter.
• $\tfrac{\Delta\phi}{2\pi}$: this shows the ratio of frequencies between the azimuthal motion and the polar motion.
• “Try to match” checkbox and field: See below in resonances. With the box checked, the r value will be automatically determined from the a value; so move the a slider around, because the r slider will be disabled.

Correction 2019-01-19: Many thanks to Charles Gammie for spotting a typo in the equation for $u_0^2$ (this typo did not affect the code).

Update 2016-03-05: Boyer-Lindquist coordinates are converted into Cartesian coordinates by treating them as oblate spheroidal coordinates, with the appropriate oblateness determined by the spin (thanks to Steve Drasco for this suggestion).

Physics

Disclaimer: since this is visualizing the trajectories of photons, this is not something anybody could “see” with their own eyes–because the only photons you see are those which make it to your eyeballs. Moreover while the photon trajectories themselves are good geometric objects, I’m plotting their Boyer-Lindquist coordinates turned into standard Cartesian coordinates. It’s very coordinate-dependent, and remember that in general relativity, coordinates don’t mean anything!

Now on to the physics. Circular photon orbits are usually ignored, because they are unstable: any tiny perturbation and the photon will end up shooting out to infinity, or plunging into the black hole. Still, they’re useful to understand (see e.g. reference 5), and they’re just fun to visualize.

Finding photon orbits firstly requires finding null geodesic curves, i.e. geodesics where the tangent vector is lightlike. We also need it to be null, $u^a u_a = 0$ where the tangent vector is $u^a \equiv dx^a/d\lambda$ for some coordinate trajectory $x^a(\lambda)$ with $\lambda$ an affine parameter.

Finding the system of ordinary differential equations is pretty straightforward (see any of the references). For example, from Teo’s paper: \begin{align} \Delta \Sigma \dot{t} &= ((r^2+a^2)^2-\Delta a^2\sin^2\theta)E - 2Mr a L_z \\ \Sigma^2 \dot{r}^2 &= E^2 r^4 + (a^2E^2-L_z^2-\mathcal{Q})r^2 + 2M[(aE-L_z)^2+\mathcal{Q}]r - a^2 \mathcal{Q} \\ \Sigma^2\dot{\theta}^2 &= \mathcal{Q} - \left( \frac{L_z^2}{\sin^2\theta}-E^2a^2 \right)\cos^2\theta \\ \Delta\Sigma\dot{\phi} &= 2MraE + (\Sigma-2Mr)\frac{L_z}{\sin^2\theta} \end{align} where overdots are derivatives with respect to an affine parameter, and with the usual Kerr quantities $\Delta = r^2-2Mr+a^2$, $\Sigma = r^2+a^2\cos^2\theta$, and where $\mathcal{Q}$ is the so-called Carter constant.

Parameters

After you study the radial equation, you learn that the only bound photon trajectories—that is, orbits!—are those for which $r=\mathrm{const}$ in Boyer-Lindquist coordinates. This is why these photon orbits are sometimes called “circular” or “spherical.”

Now to make life easier: all photons, regardless of energy, follow the same trajectories. This means we can just set $E=1$ in the above equations. This is exactly the same as doing an affine tranformation on the trajectory. This allows us to parameterize orbits in terms of two parameters instead of three: instead of $(E,L_z, \mathcal{Q})$, we use just $(\Phi,Q)$ where $\Phi \equiv L_z/E$ and $Q \equiv \mathcal{Q}/E^2$.

Not all choices of parameters give solutions. I recommend following along in Teo’s paper to understand which parts of parameter space are allowed.

In the end, you see that at each $a$, there is a one-parameter family of trajectories given by the radius $r$, which must be between the two limits $r_1(a) \le r \le r_2(a)$. The innermost photon orbit is a prograde circle lying in the equatorial plane, and the outermost orbit is a retrograde circle lying in the equatorial plane. You saw those functions previously on my Kerr calculator. You can switch between $(r,a)$ and $(\Phi,Q)$ (in the appropriate part of parameter space) via \begin{align} \Phi &= -\frac{r^3-3Mr^2+a^2 r+a^2 M}{a(r-M)} \\ Q &= - \frac{r^3(r^3-6Mr^2+9M^2r-4a^2M)}{a^2(r-M)^2} \end{align}

For any radius between the minimum and maximum, the photon will oscillate in latitude between a minimum and maximum, which you can read off of the angular potential. This minimum/maximum is given by: \[ u_0^2 = \cos^2\theta_\mathrm{min\ or\ max} = \frac{1}{2a^2} \left[ (a^2-Q-\Phi^2) + \sqrt{(a^2-Q-\Phi^2)^2+4a^2Q} \right] \] This function is plotted below:

The photon will hit every latitude value between the two blue curves at fixed $(a,r)$, making something like a truncated sphere (justifying the name “spherical” photon orbits)—unless it’s resonant!

Reparameterizing latitudinal motion

By construction, the radius is constant. Keeping track of the Boyer-Lindquist coordinate time doesn’t really add information. Therefore the only thing I integrate here is the $(\theta,\phi)$ motion.

The differential equation for $\phi$ is straightforward. But look back at the $\theta$ equation: it is quadratic in $\dot{\theta}$, so there are two roots: one positive, and one negative. These correspond to the photon ascending/descending in latitude, so if you were to use this equation you’d have to switch which root you’re using at the turning points ($\cos\theta = \pm u_0$).

That’s pretty annoying, plus you might worry about the Lipschitz condition. All of this can be avoided by reparameterizing the latitudinal motion with a periodic function of a monotonically increasing/decreasing phase angle. I learned this trick from Scott Hughes (see e.g. reference 4), though I’m not sure where it originally came from.

First, instead of using $\theta$, we switch variables to $u = \cos\theta$. The differential equation for $u$ is \[ \Sigma^2 \dot{u}^2 = a^2 (u_0^2-u^2)(u^2-u_1^2) \] where $u_1^2$ is the other root of the potential on the right hand side: \[ u_1^2 = \frac{1}{2a^2} \left[ (a^2-Q-\Phi^2) - \sqrt{(a^2-Q-\Phi^2)^2+4a^2Q} \right] \] Note that my $u_1$ is not the same as Teo’s $u_1$. Also keep track of the fact that $u_1^2 \le 0$, but I’ll stick with the name.

Now we write $u = u_0 \sin\chi(\lambda)$ with a new monotonic phase function $\chi(\lambda)$, and do a tiny bit of algebra to get \[ \dot{\chi} = \pm \frac{a}{\Sigma} \sqrt{u^2 - u_1^2} = \pm \frac{a}{r^2 + a^2 u_0^2 \sin^2 \chi} \sqrt{u_0^2 \sin^2 \chi - u_1^2} \] Because of the cancellation of the root at the turning point, there is clearly no Lipschitz problem. And, since the phase angle is monotonic, we just pick the positive root once and for all.

Resonances

For “most” choices of $(a,r)$, the frequency $\Omega_\phi$ of the azimuthal motion and the frequency $\Omega_\theta$ of the longitudinal motion are incommensurate. That is, the ratio $\Omega_\phi/\Omega_\theta$ is an irrational number “almost everywhere” (in the measure-theoretic sense). That means that most of these photon orbits will completely fill their respectively allowed truncated spheres.

However, when $\Omega_\phi/\Omega_\theta$ is the ratio of small integers, the orbit will close/recur, and the photon on that orbit will visit only a small part of the truncated sphere.

Let’s calculate $\Omega_\phi/\Omega_\theta$, which is actually the same as the extra azimuthal phase (over $2\pi$) for a complete longitudinal cycle. That is, we want to be able to compute \[ \Delta\phi = \oint_{\chi-\mathrm{cycle}} d\phi = \int_0^{2\pi} \frac{d\phi}{d\chi} d\chi = \int_0^{2\pi} \frac{d\phi/d\lambda}{d\chi/d\lambda} d\chi \] After plugging in the above equations of motion and doing a bit of algebra, this quantity can be written in terms of complete elliptic integrals (see what conventions I use for complete elliptic integrals): \[ \Delta\phi = \frac{4}{\sqrt{-u_1^2}} \left[ \frac{2r-a\Phi}{\Delta} K\left(\frac{u_0^2}{u_1^2}\right) + \frac{\Phi}{a} \Pi\left(u_0^2, \frac{u_0^2}{u_1^2}\right) \right] . \] This expression is somewhat simpler than Teo’s.

Teo’s paper goes into a bit more depth on this function, but for now let’s just examine a plot of $\Delta\phi/2\pi$:

Yes, that jump is real. That’s where the orbits switch from being prograde (on the left, positive $\Delta\phi$) to being retrograde (on the right). It’s the same point where $\Phi=0$, which is where $u_0^2=1$, which is a polar orbit. There is a jump of exactly 2 between the left and right branches of $\Delta\phi/2\pi$, and the exactly polar orbit is half way in between the two branches.

A perhaps more informative way to plot this is as contours of $\Delta\phi/2\pi$ in the $r-a$ plane, which I’ve only done statically in Mathematica (maybe at some point in the future I’ll do it in javascript). Here’s what this function looks like:

This should give you an idea of some of the easier-to-find resonances in parameter space.

If you want to see what one of these resonances looks like, check the “Try to match” checkbox and type in the desired ratio (in decimal) in the field at right. Then drag the slider for a into the appropriate range, and the simulation should solve for the appropriate value of r. The r slider will be disabled while the checkbox is checked.

If you move a outside the range where a certain resonance exists, the simulation will probably rail to $r_1$ or $r_2$.

If you choose a target value for $\Delta\phi/2\pi$ that only exists on the $\Phi=0$ line (the discontinuity), you’re gonna have a bad time—the simulation is pretty wonky along the $\Phi=0$ line, for technical reasons.

Code

This page was made using the JSXGraph toolkit, three.js, threestrap, and javascript that I wrote by hand. The latter includes code for complete elliptic integrals in javascript (which is only used for the frequency ratio calculations).

The numerical integrator from JSXGraph that I use is a fixed step-size RK4 integrator. It is sadly not adaptive! This is a problem for orbits close to the spin axis.

Limitations

• The major limitation is that this is in javascript, which is not ideal for numerical work. There are very few numerical libraries.
• JSXGraph only has fixed-step-size integrators, no adaptive ones. I picked what seemed like a reasonable step size in most of parameter space.
• Orbits that pass very close to the spin axis experience rapid variation in $\phi$ in a small amount of affine parameter. These numerics will be especially bad.
• For exactly polar orbits, the situation is even worse. In principle, $\phi$ goes through a jump discontinuity of $\pi$ when going across the pole. This isn’t possible with the formulation I’ve got here. Maybe there’s a better way.
• The iterative scheme for elliptic integrals may be bad at converging for certain inputs.
• The code is a bit of spaghetti. I have not made much effort to untangle it. Maybe in the future.
• Memory leaks. There are probably memory leaks, even though javascript is a garbage-collected language… it needs its hand held to be convinced that I really don’t need this earlier data.

References

1. Any of a number of standard textbooks on GR, e.g. MTW, Wald, or Carroll
2. Chandra’s book
3. E. Teo, Spherical photon orbits around a Kerr black hole, General Relativity and Gravitation, 2003, Volume 35, Issue 11, pp 1909-1926.
4. Some papers by Scott Hughes, where I learned the nicer parameterization for the latitudinal motion; for example, this paper.
5. Some of these calculations also appeared in Yang et al., Quasinormal-mode spectrum of Kerr black holes and its geometric interpretation, Phys. Rev. D 86, 104006 (2012) [arXiv:1207.4253]

Future

Timelike orbits? We shall see!