Gravitational waveforms for compact binaries from second-order self-force theory

Over the holidays we released our letter with the first post-adiabatic waveform from self-force theory. The comparison with numerical relativity (NR) is particularly exciting as at a mass ratio of \( q=10 \) the agreement is excellent up until close to merger — see the figure below from the letter.


There has been growing evidence that self-force results, once extended to second-order in the mass ratio, will provide accurate models of intermediate mass-ratio inspirals (IMRIs). This is the first time we have confirmed this via direct calculation. This is particular exciting as the LIGO-Virgo Collaboration are already seeing binaries with IMRI mass ratios, but currently the models they use for search and parameter estimation (PE) are not calibrated for mass ratios greater than 10:1.

The calculation of the second-order self-force needed to compute the above waveform takes many days on a computer with 40 cores (we could write more efficient code to speed this up but it will also be slow). This is an offline step though, and once the self-force is computed for a range of orbital radii this can be efficiently interpolated and then the above waveform is computed in milliseconds. That means it is almost ready-to-go for use in data analysis. The only caveat using it for IMRI search and PE is that we don’t have the transition to plunge, which contributes a very significant portion of the signal in ground-based detectors. One of our next goals is to add that piece of the waveform.

The waveform letter is in collaboration with Barry Wardell, Adam Pound, Jeremy Miller, Leanne Durkan, and Alexandre Le Tiec. It builds upon two papers published in Physical Review Letters that calculate the binding energy, and gravitational wave flux, at second-order in the mass ratio.

Gravitational-wave energy flux for compact binaries through second order in the mass ratio

We recently put out a key letter from our second-order self-force framework. In this work we compute the gravitational wave (GW) flux through second-order in the mass ratio for the first time. We do for quasi-circular orbits around an initially non-rotating black hole, but at second-order we can also account for slow rotating of the primary (large) black hole. Our flux results show remarkable agreement with numerical relativity simulations from the SXS collaboration for 10:1 mass ratios or smaller. Second-order self-force calculations are crucial to enable fundamental physics tests with extreme mass ratio inspirals (EMRIs) and our work shows they will also be useful to modelling intermediate mass ratio inspirals (IMRIs) as well. Our result are summarized by the following plot for a binary with mass ratio 10:1 \( (q=10)\).


Screenshot 2021-07-27 at 23.37.47

The figure shows the Newtonian-normalized flux for the the (2,2)-mode as a function of the (inverse) orbital separation, \(\bar{x}\) which is computed from the frequency, \(\varpi\) extracted from the waveform. Weak field orbits are to the left of the plot and strong field orbits are on the right. The location of the innermost stable circular orbit (ISCO) for geodesic motion is shown by the dashed vertical line. On the figure we plot the (2,2)-mode flux computed using various methods.

The gold standard reference is the numerical relativity result (NR) which is taken from SXS:BBH:1107. We estimate the error bars in the NR data using the different extrapolation orders they use to compute the waveform at null infinity (shown by the shaded region between the two blue curves). The oscillations in the NR result is likely due to a combination of residual eccentricity and motion of the centre of mass of the binary influencing the mode decomposition. The plot also shows the PN result which is valid in the weak field but diverges strongly from the NR result in the strong field.

The dashed, green curve is the first-order self-force result that could have been plotted since the 60s when people first solve the Regge-Wheeler equation. Our new second-order (2SF) result is the solid, red curve with triangular markers. Over much of the frequency range where we have NR data our 2SF result is within the error bars of the NR result. Only towards the merger do we see a small discrepancy before our result breaks down as the ISCO is approached. Future work will improve this by attaching a transition to plunge.

The letter goes on to show 6 more figures that show (i) the agreement between NR and 2SF at \(q=1\) is still, unexpectedly, very good, (ii) the agreement  is less good for subdominant modes, but we can resum the results to drastically improve the comparison, (iii) we also see good agreement for spin aligned binaries where the angular momentum of each black hole is small. For the secondary, so long as the mass ratio is small, this corresponds to any value of the Kerr spin parameter and for the primary this corresponds to small values of the Kerr spin parameter, and (iv) we also see agreement with all the known PN terms up to 3.5PN.

The remarkable agreement that we observe between NR and 2SF results at \(q=10\) strongly suggests that our 2SF results will be excellent models for binaries with, say, \(q=50\). Very recently we have been able the waveforms associated with the binary and at \( q=10 \) the comparison with NR looks amazing. Watch this space for more details soon.

FastEMRIWaveforms: New tools for millihertz gravitational-wave data analysis

Today we released a follow up paper to our Physical Review Letter that presented the first fully relativistic EMRI waveform model fast enough for performing parameter estimation. I outlined the model in a previous blog post. This paper provides a lot more detail and shows the utility of the new model by performing the first Markov Chain Monte Carlo-based EMRI parameter estimation with a relativistic model. We find two key results from these studies: (i) we can significantly reduce the number of harmonic modes in the model without introduction biases in the signal extraction. This helps to speed up the model even further so that, on a GPU, you can compute an EMRI waveform in under 50 milliseconds. And (ii) using semi-relativistic, ‘kludge’ waveform amplitudes introduces significant biases in the posterior distribution for some parts of the parameter space. The below corner plot shows some of the biases that can be introduced when using a Schwarzschild version of the Augmented Analytic Kludge to generate the waveform amplitudes. We gave both inspirals have exactly the same phase from an adiabatic inspiral so the differences are just due to the waveform generation step.


Screenshot 2021-04-13 at 22.05.49

Our present model is only for non-rotating black holes but it is available from a common Python interface, with all the code available from the Black Hole Perturbation Toolkit. In the LISA Data Challenges it is useful to models that cover the full parameter space of generic Kerr orbits (even if they are not phase accurate with the true waveform). For that reason we also provide an updated Augmented Analytic Kludge model via the same common interface. Both waveform models can be GPU accelerated so they generate waveforms in under a second and the common interface means we can seamless replace the kludge model with the fully relativistic model in the future.

There are a lot of nice results in the paper so I recommend you take a look (if EMRI waveforms are your thing). We are (virtual) running the BritGrav meeting at UCD this week and I spent 3 hours chair a session today so I will leave it there for now and get some sleep!

Adiabatic waveforms for extreme mass-ratio inspirals via multivoice decomposition in time and frequency

Before we get to our new paper, our previous letter was published today in Physical Review Letters. In the letter we demonstrated a method to rapidly compute fully relativistic EMRI waveforms. As part of that work we created the Fast EMRI Waveforms (FEW). Since we put out the letter the framework has been extended to have a generic waveform interface. With this new interface you can now access generic Kerr kludge waveforms computed using the Augmented Analytic Kludge (AAK). Via the generic interface you can also calculate our fast, fully relativistic Schwarzschild waveform in the detector frame. The plan is that the FEW framework will be used in future LISA Data Challenges (LDCs), and as such the older EMRI Kludge Suite is now discontinued.


Ok, on to our new paper out today. In this work we look at computing adiabatic waveforms for EMRIs in both the time and frequency domain. There are a few nice result in the paper but the main two are (i) direct calculation of EMRI waveforms in the frequency domain and (ii) the calculation of waveforms for generic inspirals into Kerr black hole. Let’s look at each of these.


Frequency domain waveforms: Often we think of EMRI waveforms as a sum of many “voices”, where each voice corresponds to a particular multipolar and frequency mode of the radiation. Each voice has a simple chirping behaviour that is similar to the well known waveforms for quasi-circular binaries. If you like you can think of an EMRI waveform as being equivalent to thousands of simultaneous chirping binaries (with a precise phase and amplitude relation between them). There is one key difference though between the quasi-circular case and the voices of the EMRI: the latter can both chirp (\(\dot{f} >0\)) and backwards chirp (\(\dot{f} < 0\)). The figure below shows the amplitude (top panel) and frequency (bottom panel) of a voice for an eccentric, equatorial Kerr inspiral with \(a=0.9, p=12, e=0.7\). In the bottom panel you can clearly see that \(\dot{f} = 0\) near the end of the inspiral

Screenshot 2021-02-05 at 16.57.08

Ok, so what? Well, as EMRIs evolve slowly we can appeal to the stationary phase approximation (SPA) to compute the Fourier transform of the waveform. Unfortunately the standard SPA breaks down when \(f=0\) and so in this work we extended the SPA to handle this case. When then apply our ‘extended SPA’ to real EMRI data and compute the frequency-domain waveform. As a check on our calculation we compare the results with the Discrete Fourier Transform (DFT) of the equivalent time-domain waveform.

Screenshot 2021-02-05 at 17.00.42

Generic inspirals into a Kerr black hole:  The second main result this paper is the calculation of waveforms (in both the time and frequency domain) for generic (eccentric and inclined) inspirals into a Kerr black hole. This involves computing the fluxes and mode amplitudes across a region of the 3D parameter space. From the fluxes we can compute the (phase space) inspiral and we can stitch the amplitudes together along this inspiral to compute the waveform in the time-domain, or frequency domain using the extended SPA discussed above. There are some additional phase corrections required in this calculation that have been discussed elsewhere but here we implement them for the first time. We verify our ‘stitched together’ waveforms against waveforms computed directly using a 2+1 time-domain Teukolsky code. The result is beautiful generic Kerr EMRI waveforms (as well as spherical and eccentric, equatorial orbits).

Screenshot 2021-02-05 at 17.08.55

What Next?  The present work can calculate frequency domain and generic Kerr waveform (driven by adiabatic fluxes). Our next goal is to incorporate this work in to the Fast EMRI Waveforms (FEW) framework. There are some new challenges here: the parameter space is now 3D (rather than 2D as in Schwarzschild) and we need methods to rapidly compute or interpolation the spheroidal harmonics. Implementing extended SPA in an efficient way will also require some thought. These are tasks for another day.


Rapid generation of fully relativistic extreme-mass-ratio-inspiral waveform templates for LISA data analysis

Our paper on rapidly computing fully relativistic EMRI waveforms is out today at arXiv:2008.06071. This is the first time that waveforms with full harmonic content can be computed on timescales useful for LISA data analysis. With our new model you can compute a year-long EMRI waveform in 10s of seconds on a CPU and in < 1s on a GPU. This is the same timing as the kludge models but without making any weak-field approximations. All the code is available in the BHPToolkit. It’s written in Python, with C++ and CUDA backends, which makes it very easy to use — see the screenshot below.


At the moment the model computes adiabatic inspirals into a Schwarzschild black hole, but the techniques used and the modular code is extensible to post-adiabatic inspirals into a Kerr black hole (watch this space). The key observation with post-adiabatic EMRI waveforms (i.e., those that include self-force corrections) is that although the waveform phase needs to be known to post-adiabatic order, the instantaneous waveform amplitudes only need to be known at leading order. Self-force corrections also add a short orbital timescale to the equations of motion. This can greatly slowdown the calculation of the inspiral but we know how to over this with either near-identity transformations or a two-timescale framework.


The main bottleneck to rapidly computing EMRI waveforms is thus interpolating the mode amplitudes and summing over thousands of modes at each time step (and each year-long waveform will have millions of time steps at LISAs sampling rate). For the interpolation we used order-reduction methods and a neural network. Neural networks are not great at high precision interpolation but we don’t need that as the error in the amplitudes does not accumulate over the inspiral. Both the neural network for the waveform amplitudes and the splines we use to sample the waveform can be evaluated extremely efficiently on a GPU. This results in year-long waveforms that can be computed in 100s of milliseconds.


The waveforms are also very accurate, with a mismatch of less than 5e-4 (relative to slow to compute high precision waveforms) across the strong-field Schwarzschild inspiral parameter space (e <= 0.7)


Finally, let me emphasize again the modular framework of the code. This should make it easy to add post-adiabatic corrections as well as environmental effects or beyond-GR/standard model physics that effect the inspiral phasing.

The (virtual) Black Hole Perturbation Toolkit Workshop

Before the lockdown we applied for and had been awarded GWVerse COST funding to host a workshop for the Black Hole Perturbation Toolkit in Prague in the Czech Republic. Unfortunately we had to put that face-to-face meeting on hold – hoping to host the workshop in March 2021 instead – but this opened a new opportunity: a virtual workshop. For the original workshop we had expected 25-30 participants but by the time the virtual workshop came around we had almost 200 participants registered! We had a range of great talks from contributors to the Toolkit and many of them also prepared awesome exercises, notebooks and extra material help people get to grips with various pieces of the Toolkit. It was amazing to see ~110 people in the Zoom room and ~20-25 people watching live on YouTube on the first day. All the talks were recorded on YouTube and can be watched again on the BHPToolkit YouTube channel and the talk material (slides, notebooks) can be found on the workshop website.

At the time of writing the videos have already amassed over 700 views. The first two days of the workshop were devoted to introducing the tools to new users with an emphasis on giving worked examples and exercises showing how the code works and how to use the code in your own projects. We also had short contributed talks from 8-9 people who spoke for 5 mins each on a topic related to, or already using, the Toolkit. It was awesome to meet new people using our code to make interesting calculations. The third day was the developer day and we set up a Slack channel to facilitate interactions between developers and potential new developers. This was a great success and already new code has been added to the Toolkit and I can see a lot more coming by looking at the spike of activity in the development branches.

The virtual workshop format worked far better than I could have ever imagined. A huge thank-you to our “local” organizing committee (Vojtech Witzany and Georgios Lukes-Gerakopoulos) for all their work behind the scenes; a huge thank you to all the speakers who gave consistently excellent talks that they had clearly put a lot of effort into preparing; and finally a huge thank you for all the attendees who turns up on mass and engaged throughout. It’s clear that virtual workshops give an opportunity to reach a large audience: we had 84 attendees from Europe, 37 from North America, 16 from Asia, and 5 from South America (rough stats based on an analysis of the email addresses people used to sign up – usually their host institution addresses but I had to remove all the .coms etc so not a perfect analysis). The map below shows the countries, based on my rough analysis, that attendees came from.



Given the huge success of the virtual workshop I think that even once the lockdowns are all over we will give serious consideration to running a virtual workshop again.

(virtual) Talk at the Gravity Exploration Institute at Cardiff University

I gave my second virtual talk of the lockdown today at the Gravity Exploration Institute of Cardiff University. Screenshot 2020-05-19 19.01.42


The talk was on three techniques for rapidly computing accurate EMRI waveforms. By combining these techniques we will be able to compute year-long EMRI waveforms in under a second. The three techniques are:

  • Non-linear black hole perturbation theory. This is needed in order to compute the post-adiabatic corrections to the waveform phase that will leave a residual error in the phase that scales with the mass ratio (so it the error in the waveform phase should be tiny for an EMRI). Our first paper on this topic presenting a result for a binary was recently published in Physical Review Letters. In this talk I also presented some exciting new (preliminary) results where we have computed the second-order (in the mass-ratio) flux for the first time. There are still many checks to go (and for some reason the l=2,m=2 mode is not working yet) but it feels great to be reaching a major goal with this long term project.
  • Near-identity (averaging) transformations. Once the oscillatory pieces of the self-force are included in the inspiral models the phase space trajectory and the orbital phases oscillate (with the number of oscillations scaling with the mass ratio). Numerically computing the phase space trajectory and the phases thus becomes very slow as the integrator has to resolve tens to hundreds of thousands of small oscillations. This can be overcome using so-called near-identity (averaging) transformations, or NITs. These average out the short timescale physics to produce equations of motion that capture the correct long-term behaviour of the system but that do not oscillate on the orbital timescale. The resulting NIT’ed equations of motion can be solved in milliseconds rather than minutes. Maarten van de Meent and I published a paper on this in Classical and Quantum Gravity, which looked at the general transformation and applied to to inspirals in Schwarzschild spacetimes. Since then one of my PhD students, Philip Lynch, has explicitly derived the NIT for generic inspiral into a Kerr black hole and has a numerical model implemented for equatorial orbits in Kerr spacetime.
  • Fast EMRI waveforms. To compute an EMRI waveform you need to sum over many thousands of harmonics at each time step. A straightforward implementation of this approach is quite slow, resulting in year-long waveforms that take 10s-100s of second to produce. In collaboration with Alvin Chua, Michael L. Katz, and Scott Hughes, I presented our new method for computing year-long EMRI waveforms in ~200ms. This involves a computation of a reduced-order model (ROM), neural networks and acceleration using graphics processor units (GPUs).

For those who are interested here’s a link to the slides from my talk.

The show must go on

Recently I was invited to give a seminar as part of the Maynooth Theoretical Physics Seminar Series and with the lockdown it went online, so earlier today I gave my first virtual seminar. On the whole it worked pretty well. My talk was titled ‘Gravitational Waves from Galactic Cores’ and reviewed modelling of small mass ratio binaries. Towards the end I highlighted our recent result on second-order calculations (now published in Physical Review Letters). I also included material highlighting our work from last spring on gravitational waves from objects orbiting Sgr A*, which unfortunately I didn’t get time to talk about during the seminar.

Screenshot 2020-04-17 18.14.04

For those who are interested here’s a link to the slides from the talk.

Talk at Caltech

I recently visited Caltech to work with Alvin Chua (and Michael Katz, who was also visiting) on a new fast method to compute full relativistic EMRI waveforms that can incorporate all known self-force corrections to the phase. We made great progress during the week and are now writing up that work. I’ll post on that again when the paper is out.

Whilst at Caltech I gave a talk on our recent second-order in the mass-ratio result. That talk was recorded in case you are interested:

One week, two papers

It’s been a busy end to the year and this week I had two papers out on the arXiv. I briefly review them below.

The location of the last stable orbit in Kerr spacetime

The first was with Leo Stein. One of the key differences between Newtonian and Einstein gravity is the presence of unstable and plunging orbits around compact objects. The surface in parameter space which separates stable and unstable orbits is called the ‘separatrix’. The location of the separatrix is non-trivial and in this work we present a new approach to finding it based upon finding roots of a high-order polynomial. This allows for analytic study of the separatrix. This allows us to easily derive many limits in previously calculated in the literature (such as the ISCO, ISSO, IBSO, etc). The polynomial is of order 12 in ‘p’ (the semi-latus rectum) so there is are closed form solutions in general but we provide a robust numerical scheme for finding the roots.


Screenshot 2019-12-18 11.13.02

The above figure shows the an extended view in (p,e,x) space of the solutions of the separatrix polynomial. The physical separatrix is one of the sheets (restricted to 0 <= e <=1).


Dissipation in extreme-mass ratio binaries with a spinning secondary

In this work we examine the gravitational-wave flux from a compact binary where the secondary is spinning. We derive the balance law in this case, which relates the flux to the local change in the secondary’s orbital constants of motion (energy, angular momentum). We then perform an explicit calculation for the case of a spin-aligned binary where the primary is a non-rotating (Schwarzschild) black hole. We make this calculation using the Teukolsky formalism reconstructing in the radiation gauge for the local calculation (we make a post-Newtonian and numerical calculation using this method), and a numerical Lorenz-gauge calculation. We find perfect agreement between the approaches and show explicitly that our balance law holds.

Screenshot 2019-12-20 08.23.53Screenshot 2019-12-20 08.24.09

The above two plots show the comparison of the flux to infinity between the PN and numerical results, and the Lorenz-gauge metric perturbation for the (2,2)-mode for a spinning body. All the PN and numerical results can be found in the Black Hole Perturbation Toolkit.