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.

FEW_waveform

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.

Efnzer-XgAEFt5N

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.

FEW_timing

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)

FEW_accuracy

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.

BHPToolkitWorkshopAttendeeMap

 

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.

Second-order self-force calculation of gravitational binding energy

It is with great pleasure that we present today the first second-order (in the mass-ratio) self-force result. From the monopole piece of the second-order metric perturbation, and specializing to circular orbits in Schwarzschild spacetime, we compute the binding energy of the binary. We have compared the results of our calculation with the result from the first law of binary mechanics (FLBM) and find great agreement — see the figure below.

Screenshot 2019-08-21 09.34.25

In comparing against the FLBM we find a slight discrepancy in the very strong field. This discrepancy is not unexpected as the FLBM only applies to a conservative spacetime (no inspiral) whereas we include the effects of dissipation in our calculation.

This work builds upon almost 7 years of effort since the foundational second-order papers. It’s been a long road with 10+ papers already out and I can think of at least another 5 in preparation.

This first concrete calculation really marks the end of the beginning. As we compute higher modes we will get the binding angular momentum, the radiated flux, gauge invariants, local force, and then the first complete post-adiabatic waveform. Watch this space!

RailJourneysMap-small

Traveling by train

I think my favourite way to travel is by train so when the opportunity arises, if possible, I will try to take the train there, back or both. A good example of this is that I was fortunate enough to be invited to lecture at a summer school in Beijing in 2016. Whilst I flew out there the only logical choice was to take the train back (or at least as far as St. Petersburg).

One thing I’ve been thinking of doing for a while now is making map of all the rail journey’s I’ve been on. This has resulted in my Rail Journey Map which shows all the intercity rail journeys I have been on totalling ~24,000 km of travel.

Compiling this map was rather challenging and it took me some time to research how to do it. Below is the method I used. It might not be the best, but it worked for me.

Creating the rail journey map

To create the map made use of a piece of open-source software called OpenRailRouting. The instructions are pretty clear on GitHub page about how to install and run it but a key piece of information is missing, at least for those not used to using OpenStreetMap (OSM) data.

To get the OSM data I had to download the 44GB planet PBF file and then use the Osmium tool to extract all the railway information. The command for this is:

osmium tags-filter -o planet-rail.osm.pbf planet.osm.pbf nw/railway

This results in a much more manageable ~350MB file. With this in hand you can run OpenRailRouting as described on the software’s webpage. This will setup a GraphHopper server locally which you can reach at localhost:8989. With this GUI you can select start and end points on a rail journey and the software will provide the route in between. You can also set intermediate points in case the software returns another route you didn’t take. Once you’ve settled on a route, you can export the GPX file by clicking the little ‘gpx’ button.

Finally, to plot all the routes together I used Mathematica as this made it very easy to change the map projection.

One useful thing to note is that the default configuration for OpenRailRouting only routes on standard gauge (1435mm) tracks. Many countries use this rail gauge but some, such as Ireland, Russia and Mongolia, use different gauges. You can modify the config file to add additional gauges easily.