Kerr orbit visualizer

In March I wrote a code to compute the frequencies and constants of motion associated with generic, bound, timelike Kerr geodesics. Back then I said I wanted to write an online tool for visualizing the associated orbits. A few weeks back I did just that and I finally have time to share it and write a little about it now. Rather than trying to embed it in the WordPress layout the Kerr timelike orbit visualizer can be found here.

The tool is basic but it will plot most orbits. Note it wont tell you if you enter parameters that do not correspond to Kerr a bound geodesic, you’ll just get a blank plot. Also, special cases like \(a=0, e=0\) and \(\theta_\text{inc}=0\) are not implemented (but you can set a very close to zero value to get it to work).

You can plot the orbit in Boyer-Lindquist coordinates and also in co-rotating coordinates [as defined in Eq. (3.2) of arXiv:0904.3810] where the orbit usually looks much simpler. A nice little feature, which works in most modern browsers, is you can `play’ the orbital frequencies*. I like really like this feature as it allows you to hear the character of each orbit. If you find a set of parameters near a resonance you can hear the beating between the fundamental frequencies.

This visualizer is pretty basic but I’ve had discussions with Leo Stein and Scott Hughes about improving it. In particular, it would be nice to make it more user friendly, provide more information as the orbits are plotted, give a nice set of example orbits, plot the black hole, allow animation of the orbit, etc. It might also be possible to speed up the calculating using fast Fourier transform methods.

The current code uses plotly.js to visualize the orbit. The orbit is computed by numerically integrating the geodesic equations with a fixed-step RK4 integrator found in JSXgraph. Much of the inspiration to create this online visualizer came from Leo Stein’s visualizer for bound, spherical null geodesics in Kerr. Documenting the equations the code solves, as Leo does so nicely, is something that also needs to be done.

* the frequencies do not correspond to any astrophysical extreme mass-ratio binary as such systems have frequencies in the milli-Hertz regime which cannot be heard by humans. Instead I just increase all the frequencies by a constant multiple to make them audible.


Invited lecturer at Spring School in 北京 (Beijing)

I have been invited to lecture on numerical approaches to black hole perturbation theory and the self-force problem at a spring summer school in Beijing. The two-week long school is focusing on numerical relativity and gravitational waves. The school will take place from 15th-25th of May at the Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing and is aimed at PhD students and postdocs. More information about the school can be found here:

Other lectures include: David Hilditch (Friedrich Schiller University of Jena), Koutarou Kyutoku (KEK, IPNS), Matthias Hanauske (Goethe University Frankfurt), Andrea Taracchini (Max Planck Institute for Gravitational Physics), and David Weir (University of Helsinki). They will cover topics such as post-Newtonian (PN) theory, effective-one-body theory (EOB), numerical relativity (NR), and gravitational wave sources such as neutron star (NS) binaries, black hole (BH)-NS binaries, and exotics sources.

My own lectures will concentrate on black hole perturbation as applied to small-mass ratio binaries. I’m still sketching the content of the individual lectures but once that is done I will update this page.

This will be my second visit to Beijing. I had a great time first time I went (as part of China-Japan trip in 2013) and am looking forward to heading out there again in May.


Bound timelike Kerr geodesics: constants of motion and frequencies

\(\theta_\text{inc}:\) deg


Brief Explanation

Any timelike geodesic about a Kerr black hole has associated with it three constants of motion: energy \(\mathcal{E}\), angular momentum \(\mathcal{L}\) and the Carter constant \(\mathcal{Q}\). In addition, generic bound orbits have three frequencies associated with the orbital motion. With respect to Boyer-Lindquist coordinate time these are the rate of azimuthal angle accumulation \(\Omega_\varphi\), radial libration \(\Omega_r\) and polar libration \(\Omega_\theta\). These frequencies can also be defined with respect to Mino time. The Mino time frequencies, \(\Upsilon_\alpha\), can be calculated via \(\Upsilon_\alpha = \Gamma \Omega_\alpha \). The above interactive Javascript code computes these constants and frequencies given the orbital parameters.

Method, Code and Limitations

The equations to compute the constants of motion come from W. Schmidt’s Celestial mechanics in Kerr spacetime.

A method to compute the orbital frequencies can also be found in that work but that method requires numerically integrating various equations. Instead I use equations in terms of elliptic integrals from Fujita and Hikida’s Analytical solutions of bound timelike geodesic orbits in Kerr spacetime.

To compute elliptic integrals in Javascript I use Leo Stein‘s library which can be found on GitHub.

Limitations: Not all parameters will lead to bound stable orbits. The code will not warn you of this, in part because calculating when a set of parameters leads to a stable orbit is not straight forward for generic orbits. Instead you will get NaN’s (Not a Number). Also the code gives an error for the Schwarzschild \(a=0\) case. It’s on my todo list to fix that (note you can put in very small values of \(a\)).


I also wrote, and am making freely available, code to compute the orbital constants and frequencies in:

I also have C and C++ (arbitrary precision using Arb) implementations but they need a little tidying up before public release.

Future work

Plot the associated orbits in an interactive manner similar to Leo’s nice visualization of photon orbits. These orbits have a very rich structure, e.g.,


This orbit has parameters \((a,p,e,\theta_\text{min}) = (0.998,3,0.7,\pi/4) \). The apparent complexity of the orbit can be made simpler by moving to a co-rotating coordinate system. This resulting orbit is quite pretty:


I already have Mathematica code to make the calculation (that code produced the above two images). The Mathematica code could easily be converted to Python but Javascript would be more difficult due to a lack of good numerical libraries.


X-ray reverberation

In 2016 I spent a little time exploring x-ray reverberation with colleagues from the University of Southampton. We’re interested in modeling flares in accreting black hole systems. Initially we’re considering a simplified lamppost model where the flare takes place along the symmetry axis of the system. Our goal is to calculate the difference in the time of arrival of the photons that travel directly from the flare to the observer (us on Earth) and the photons that reverberate around the black hole + accretion disk system (this includes photons that are bent around the black hole and photons that ‘reflect’ off the accretion disk).

The project, which is still in its infancy, involves ray-tracing a black hole spacetime with an accretion disk. This involves solving the null geodesic equations, working out where and when the photons interact with the accretion disk and calculating the time of arrival. As well as the time of arrival we’re interested in the intensity of the emission and its spectrum. Currently the first of these (time of arrival) is implemented for a Schwarzschild black hole with a simple accretion disk.

The code is written in Mathematica and uses semi-analytic methods to calculate the path of the light. Using the code I made a few short animations which can be found here: The page describes what each animation depicts.


Research trip to Rio de Janeiro

At the start of 2016, myself and Prof. Marc Casals applied to the Brazil-US Exchange program funded by the Sociedade Brasileira de Física (SBF) and the American Physical Society (APS). Our proposal was accepted and I was awarded up to $3000 to travel to Rio to work with Marc for two weeks, within one year of the grant being awarded.

We completed our two weeks of intensive research in the first two weeks of October 2016. Whilst I was there we worked on implementing new methods for computing analytic series solutions to the Teukolsky radial equation. This was my second trip to Rio to visit Marc. I always find the second time you visit a city to be very enjoyable; the first time you run around doing all the touristy stuff, but on the second trip you can relax into the city and experience a little more of what it would be like to live there.

Rio is a really good looking city. Nestled between mountains and on the Atlantic coast it just looks good, particularly if you can get up anywhere high. The people I met at the institute – Centro Brasileiro de Pesquisas Físicas (CBPF) – and elsewhere were friendly and welcoming. It was a very pleasant place to spend two weeks working.

If you are interested in apply for this program yourself more details can be found on the APS Brazil-U.S. Exchange Program webpage. They also have a program for collabration between the US and India in case that interests you (which it should as India is another awesome place to visit).

Screenshot 2016-08-25 01.09.43

Gargantua at CQG+

I wrote an entry for CQG+, the companion blog to the journal Classical and Quantum Gravity. The article is aimed at someone with a general science background and includes some nice sounds comparing a standard chirp waveform with a (non-chirping) Gargantua waveform. Head over to CQG+ to check it out.

A nice upshot of the CQG+ posting is that the published version of the paper is open access for the next month. After that time has passed the article can still be found, free to read, on

Screenshot 2016-06-15 14.18.15

Second gravitational-wave transient: GW151226

The LIGO Scientific Collaboration and Virgo Collaboration just announced the second detection of gravitational waves. The detection paper can be found in Physical Review Letters. As with the first event the source was the merger of two black holes. This time the mass of the two components was 14.2 and 7.5 solar masses and the signal was observed for 55 cycles lasting ~1 second.

There will be a lot written about this event in the coming days. There is a nice summary at