# Dirichlet Process Gaussian-mixture model: An application to localizing coalescing binary neutron stars with gravitational-wave observations

Where do gravitational waves like GW170817 come from? Using our network of detectors, we cannot pinpoint a source, but we can make a good estimate—the amplitude of the signal tells us about the distance; the time delay between the signal arriving at different detectors, and relative amplitudes of the signal in different detectors tells us about the sky position (see the excellent video by Leo Singer below).

In this paper we look at full three-dimensional localization of gravitational-wave sources; we important a (rather cunning) technique from computer vision to construct a probability distribution for the source’s location, and then explore how well we could localise a set of simulated binary neutron stars. Knowing the source location enables lots of cool science. First, it aids direct follow-up observations with non-gravitational-wave observatories, searching for electromagnetic or neutrino counterparts. It’s especially helpful if you can cross-reference with galaxy catalogues, to find the most probable source locations (this technique was used to find the kilonova associated with GW170817). Even without finding a counterpart, knowing the most probable host galaxy helps us figure out how the source formed (have lots of stars been born recently, or are all the stars old?), and allows us measure the expansion of the Universe. Having a reliable technique to reconstruct source locations is useful!

This was a fun paper to write [bonus note]. I’m sure it will be valuable, both for showing how to perform this type of reconstruction of a multi-dimensional probability density, and for its implications for source localization and follow-up of gravitational-wave signals. I go into details of both below, first discussing our statistical model (this is a bit technical), then looking at our results for a set of binary neutron stars (which have implications for hunting for counterparts) .

### Dirichlet process Gaussian mixture model

When we analyse gravitational-wave data to infer the source properties (location, masses, etc.), we map out parameter space with a set of samples: a list of points in the parameter space, with there being more around more probable locations and fewer in less probable locations. These samples encode everything about the probability distribution for the different parameters, we just need to extract it…

For our application, we want a nice smooth probability density. How do we convert a bunch of discrete samples to a smooth distribution? The simplest thing is to bin the samples. However, picking the right bin size is difficult, and becomes much harder in higher dimensions. Another popular option is to use kernel density estimation. This is better at ensuring smooth results, but you now have to worry about the size of your kernels.

Our approach is in essence to use a kernel density estimate, but to learn the size and position of the kernels (as well as the number) from the data as an extra layer of inference. The “Gaussian mixture model” part of the name refers to the kernels—we use several different Gaussians. The “Dirichlet process” part refers to how we assign their properties (their means and standard deviations). What I really like about this technique, as opposed to the usual rule-of-thumb approaches used for kernel density estimation,  is that it is well justified from a theoretical point of view.

I hadn’t come across a Dirchlet process before. Section 2 of the paper is a walkthrough of how I built up an understanding of this mathematical object, and it contains lots of helpful references if you’d like to dig deeper.

In our application, you can think of the Dirichlet process as being a probability distribution for probability distributions. We want a probability distribution describing the source location. Given our samples, we infer what this looks like. We could put all the probability into one big Gaussian, or we could put it into lots of little Gaussians. The Gaussians could be wide or narrow or a mix. The Dirichlet distribution allows us to assign probabilities to each configuration of Gaussians; for example, if our samples are all in the northern hemisphere, we probably want Gaussians centred around there, rather than in the southern hemisphere.

With the resulting probability distribution for the source location, we can quickly evaluate it at a single point. This means we can rapidly produce a list of most probable source galaxies—extremely handy if you need to know where to point a telescope before a kilonova fades away (or someone else finds it).

### Gravitational-wave localization

To verify our technique works, and develop an intuition for three-dimensional localizations, we used studied a set of simulated binary neutron star signals created for the First 2 Years trilogy of papers. This data set is well studied now, it illustrates performance it what we anticipated to be the first two observing runs of the advanced detectors, which turned out to be not too far from the truth. We have previously looked at three-dimensional localizations for these signals using a super rapid approximation.

The plots below show how well we could localise the sources of our binary neutron star sources. Specifically, the plots show the size of the volume which has a 90% probability of containing the source verses the signal-to-noise ratio (the loudness) of the signal. Typically, volumes are $10^4$$10^5~\mathrm{Mpc}^3$, which is about $10^{68}$$10^{69}$ Olympic swimming pools. Such a volume would contain something like $100$$1000$ galaxies.

Localization volume as a function of signal-to-noise ratio. The top panel shows results for two-detector observations: the LIGO-Hanford and LIGO-Livingston (HL) network similar to in the first observing run, and the LIGO and Virgo (HLV) network similar to the second observing run. The bottom panel shows all observations for the HLV network including those with all three detectors which are colour coded by the fraction of the total signal-to-noise ratio from Virgo. In both panels, there are fiducial lines scaling inversely with the sixth power of the signal-to-noise ratio. Adapted from Fig. 4 of Del Pozzo et al. (2018).

Looking at the results in detail, we can learn a number of things

1. The localization volume is roughly inversely proportional to the sixth power of the signal-to-noise ratio [bonus note]. Loud signals are localized much better than quieter ones!
2. The localization dramatically improves when we have three-detector observations. The extra detector improves the sky localization, which reduces the localization volume.
3. To get the benefit of the extra detector, the source needs to be close enough that all the detectors could get a decent amount of the signal-to-noise ratio. In our case, Virgo is the least sensitive, and we see the the best localizations are when it has a fair share of the signal-to-noise ratio.
4. Considering the cases where we only have two detectors, localization volumes get bigger at a given signal-to-noise ration as the detectors get more sensitive. This is because we can detect sources at greater distances.

Putting all these bits together, I think in the future, when we have lots of detections, it would make most sense to prioritise following up the loudest signals. These are the best localised, and will also be the brightest since they are the closest, meaning there’s the greatest potential for actually finding a counterpart. As the sensitivity of the detectors improves, its only going to get more difficult to find a counterpart to a typical gravitational-wave signal, as sources will be further away and less well localized. However, having more sensitive detectors also means that we are more likely to have a really loud signal, which should be really well localized.

Left: Localization (yellow) with a network of two low-sensitivity detectors. The sky location is uncertain, but we know the source must be nearby. Right: Localization (green) with a network of three high-sensitivity detectors. We have good constraints on the source location, but it could now be at a much greater range of distances. Not to scale.

Using our localization volumes as a guide, you would only need to search one galaxy to find the true source in about 7% of cases with a three-detector network similar to at the end of our second observing run. Similarly, only ten would need to be searched in 23% of cases. It might be possible to get even better performance by considering which galaxies are most probable because they are the biggest or the most likely to produce merging binary neutron stars. This is definitely a good approach to follow.

Galaxies within the 90% credible volume of an example simulated source, colour coded by probability. The galaxies are from the GLADE Catalog; incompleteness in the plane of the Milky Way causes the missing wedge of galaxies. The true source location is marked by a cross [bonus note]. Part of Figure 5 of Del Pozzo et al. (2018).

arXiv: 1801.08009 [astro-ph.IM]
Journal: Monthly Notices of the Royal Astronomical Society; 479(1):601–614; 2018
Code: 3d_volume
Buzzword bingo: Interdisciplinary (we worked with computer scientist Tom Haines); machine learning (the inference involving our Dirichlet process Gaussian mixture model); multimessenger astronomy (as our results are useful for following up gravitational-wave signals in the search for counterparts)

### Bonus notes

#### Writing

We started writing this paper back before the first observing run of Advanced LIGO. We had a pretty complete draft on Friday 11 September 2015. We just needed to gather together a few extra numbers and polish up the figures and we’d be done! At 10:50 am on Monday 14 September 2015, we made our first detection of gravitational waves. The paper was put on hold. The pace of discoveries over the coming years meant we never quite found enough time to get it together—I’ve rewritten the introduction a dozen times. It’s extremely satisfying to have it done. This is a shame, as it meant that this study came out much later than our other three-dimensional localization study. The delay has the advantage of justifying one of my favourite acknowledgement sections.

#### Sixth power

We find that the localization volume $\Delta V$ is inversely proportional to the sixth power of the signal-to-noise ration $\varrho$. This is what you would expect. The localization volume depends upon the angular uncertainty on the sky $\Delta \Omega$, the distance to the source $D$, and the distance uncertainty $\Delta D$,

$\Delta V \sim D^2 \Delta \Omega \Delta D$.

Typically, the uncertainty on a parameter (like the masses) scales inversely with the signal-to-noise ratio. This is the case for the logarithm of the distance, which means

$\displaystyle \frac{\Delta D}{D} \propto \varrho^{-1}$.

The uncertainty in the sky location (being two dimensional) scales inversely with the square of the signal-to-noise ration,

$\Delta \Omega \propto \varrho^{-2}$.

The signal-to-noise ratio itself is inversely proportional to the distance to the source (sources further way are quieter. Therefore, putting everything together gives

$\Delta V \propto \varrho^{-6}$.

#### Treasure

We all know that treasure is marked by a cross. In the case of a binary neutron star merger, dense material ejected from the neutron stars will decay to heavy elements like gold and platinum, so there is definitely a lot of treasure at the source location.

# GW170814—Enter Virgo

On 14 August 2017 a gravitational wave signal (GW170814), originating from the coalescence of a binary black hole system, was observed by the global gravitational-wave observatory network of the two Advanced LIGO detectors and Advanced Virgo.  That’s right, Virgo is in the game!

Very few things excite me like unlocking a new character in Smash Bros. A new gravitational wave observatory might come close.

Advanced Virgo joined O2, the second observing run of the advanced detector era, on 1 August. This was a huge achievement. It has not been an easy route commissioning the new detector—it never ceases to amaze me how sensitive these machines are. Together, Advanced Virgo (near Pisa) and the two Advanced LIGO detectors (in Livingston and Hanford in the US) would take data until the end of O2 on 25 August.

On 14 August, we found a signal. A signal that was observable in all three detectors [bonus note]. Virgo is less sensitive than the LIGO instruments, so there is no impressive plot that shows something clearly popping out, but the Virgo data do complement the LIGO observations, indicating a consistent signal in all three detectors [bonus note].

A cartoon of three different ways to visualise GW170814 in the three detectors. These take a bit of explaining. The top panel shows the signal-to-noise ratio the search template that matched GW170814. They peak at the time corresponding to the merger. The peaks are clear in Hanford and Livingston. The peak in Virgo is less exceptional, but it matches the expected time delay and amplitude for the signal. The middle panels show time–frequency plots. The upward sweeping chirp is visible in Hanford and Livingston, but less so in Virgo as it is less sensitive. The plot is zoomed in so that its possible to pick out the detail in Virgo, but the chirp is visible for a longer stretch of time than plotted in Livingston. The bottom panel shows whitened and band-passed strain data, together with the 90% region of the binary black hole templates used to infer the parameters of the source (the narrow dark band), and an unmodelled, coherent reconstruction of the signal (the wider light band) . The agreement between the templates and the reconstruction is a check that the gravitational waves match our expectations for binary black holes. The whitening of the data mirrors how we do the analysis, by weighting noise at different frequency by an estimate of their typical fluctuations. The signal does certainly look like the inspiral, merger and ringdown of a binary black hole. Figure 1 of the GW170814 Paper.

The signal originated from the coalescence of two black holes. GW170814 is thus added to the growing family of GW150914, LVT151012, GW151226 and GW170104.

GW170814 most closely resembles GW150914 and GW170104 (perhaps there’s something about ending with a 4). If we compare the masses of the two component black holes of the binary ($m_1$ and $m_2$), and the black hole they merge to form ($M_\mathrm{f}$), they are all quite similar

• GW150914: $m_1 = 36.2^{+5.2}_{-3.8} M_\odot$, $m_2 = 29.1^{+3.7}_{-4.4} M_\odot$, $M_\mathrm{f} = 62.3^{+3.7}_{-3.1} M_\odot$;
• GW170104: $m_1 = 31.2^{+5.4}_{-6.0} M_\odot$, $m_2 = 19.4^{+5.3}_{-5.9} M_\odot$, $M_\mathrm{f} = 48.7^{+5.7}_{-4.6} M_\odot$;
• GW170814: $m_1 = 30.5^{+5.7}_{-3.0} M_\odot$, $m_2 = 25.3^{+2.8}_{-4.2} M_\odot$, $M_\mathrm{f} = 53.2^{+3.2}_{-2.5} M_\odot$.

GW170814’s source is another high-mass black hole system. It’s not too surprising (now we know that these systems exist) that we observe lots of these, as more massive black holes produce louder gravitational wave signals.

GW170814 is also comparable in terms of black holes spins. Spins are more difficult to measure than masses, so we’ll just look at the effective inspiral spin $\chi_\mathrm{eff}$, a particular combination of the two component spins that influences how they inspiral together, and the spin of the final black hole $a_\mathrm{f}$

• GW150914: $\chi_\mathrm{eff} = -0.06^{+0.14}_{-0.14}$, $a_\mathrm{f} = 0.70^{+0.07}_{-0.05}$;
• GW170104:$\chi_\mathrm{eff} = -0.12^{+0.21}_{-0.30}$, $a_\mathrm{f} = 0.64^{+0.09}_{-0.20}$;
• GW170814:$\chi_\mathrm{eff} = 0.06^{+0.12}_{-0.12}$, $a_\mathrm{f} = 0.70^{+0.07}_{-0.05}$.

There’s some spread, but the effective inspiral spins are all consistent with being close to zero. Small values occur when the individual spins are small, if the spins are misaligned with each other, or some combination of the two. I’m starting to ponder if high-mass black holes might have small spins. We don’t have enough information to tease these apart yet, but this new system is consistent with the story so far.

One of the things Virgo helps a lot with is localizing the source on the sky. Most of the information about the source location comes from the difference in arrival times at the detectors (since we know that gravitational waves should travel at the speed of light). With two detectors, the time delay constrains the source to a ring on the sky; with three detectors, time delays can narrow the possible locations down to a couple of blobs. Folding in the amplitude of the signal as measured by the different detectors adds extra information, since detectors are not equally sensitive to all points on the sky (they are most sensitive to sources over head or underneath). This can even help when you don’t observe the signal in all detectors, as you know the source must be in a direction that detector isn’t too sensitive too. GW170814 arrived at LIGO Livingston first (although it’s not a competition), then ~8 ms later at LIGO Hanford, and finally ~14 ms later at Virgo.  If we only had the two LIGO detectors, we’d have an uncertainty on the source’s sky position of over 1000 square degrees, but adding in Virgo, we get this down to 60 square degrees. That’s still pretty large by astronomical standards (the full Moon is about a quarter of a square degree), but a fantastic improvement [bonus note]!

90% probability localizations for GW170814. The large banana shaped (and banana coloured, but not banana flavoured) curve uses just the two LIGO detectors, the area is 1160 square degrees. The green shows the improvement adding Virgo, the area is just 100 square degrees. Both of these are calculated using BAYESTAR, a rapid localization algorithm. The purple map is the final localization from our full parameter estimation analysis (LALInference), its area is just 60 square degrees! Whereas BAYESTAR only uses the best matching template from the search, the full parameter estimation analysis is free to explore a range of different templates. Part of Figure 3 of the GW170814 Paper.

Having additional detectors can help improve gravitational wave measurements in other ways too. One of the predictions of general relativity is that gravitational waves come in two polarizations. These polarizations describe the pattern of stretching and squashing as the wave passes, and are illustrated below.

The two polarizations of gravitational waves: plus (left) and cross (right). Here, the wave is travelling into or out of the screen. Animations adapted from those by MOBle on Wikipedia.

These two polarizations are the two tensor polarizations, but other patterns of squeezing could be present in modified theories of gravity. If we could detect any of these we would immediately know that general relativity is wrong. The two LIGO detectors are almost exactly aligned, so its difficult to get any information on other polarizations. (We tried with GW150914 and couldn’t say anything either way). With Virgo, we get a little more information. As a first illustration of what we may be able to do, we compared how well the observed pattern of radiation at the detectors matched different polarizations, to see how general relativity’s tensor polarizations compared to a signal of entirely vector or scalar radiation. The tensor polarizations are clearly preferred, so general relativity lives another day. This isn’t too surprising, as most modified theories of gravity with other polarizations predict mixtures of the different polarizations (rather than all of one). To be able to constrain all the  mixtures with these short signals we really need a network of five detectors, so we’ll have to wait for KAGRA and LIGO-India to come on-line.

The six polarizations of a metric theory of gravity. The wave is travelling in the $z$ direction. (a) and (b) are the plus and cross tensor polarizations of general relativity. (c) and (d) are the scalar breathing and longitudinal modes, and (e) and (f) are the vector $x$ and $y$ polarizations. The tensor polarizations (in red) are transverse, the vector and longitudinal scalar mode (in green) are longitudinal. The scalar breathing mode (in blue) is an isotropic expansion and contraction, so its a bit of a mix of transverse and longitudinal. Figure 10 from (the excellent) Will (2014).

We’ll be presenting a more detailed analysis of GW170814 later, in papers summarising our O2 results, so stay tuned for more.

Title: GW170814: A three-detector observation of gravitational waves from a binary black hole coalescence
arXiv: 1709.09660 [gr-qc]
Journal: Physical Review Letters; 119(14):141101(16) [bonus note]
Data release: LIGO Open Science Center
Science summary: GW170814: A three-detector observation of gravitational waves from a binary black hole coalescence

If you’re looking for the most up-to-date results regarding GW170814, check out the O2 Catalogue Paper.

### Bonus notes

#### Signs of paranoia

Those of you who have been following the story of gravitational waves for a while may remember the case of the Big Dog. This was a blind injection of a signal during the initial detector era. One of the things that made it an interesting signal to analyse, was that it had been injected with an inconsistent sign in Virgo compared to the two LIGO instruments (basically it was upside down). Making this type of sign error is easy, and we were a little worried that we might make this sort of mistake when analysing the real data. The Virgo calibration team were extremely careful about this, and confident in their results. Of course, we’re quite paranoid, so during the preliminary analysis of GW170814, we tried some parameter estimation runs with the data from Virgo flipped. This was clearly disfavoured compared to the right sign, so we all breathed easily.

I am starting to believe that God may be a detector commissioner. At the start of O1, we didn’t have the hardware injection systems operational, but GW150914 showed that things were working properly. Now, with a third detector on-line, GW170814 shows that the network is functioning properly. Astrophysical injections are definitely the best way to confirm things are working!

#### Signal hunting

Our usual way to search for binary black hole signals is compare the data to a bank of waveform templates. Since Virgo is less sensitive the the two LIGO detectors, and would only be running for a short amount of time, these main searches weren’t extended to use data from all three detectors. This seemed like a sensible plan, we were confident that this wouldn’t cause us to miss anything, and we can detect GW170814 with high significance using just data from Livingston and Hanford—the false alarm rate is estimated to be less than 1 in 27000 years (meaning that if the detectors were left running in the same state, we’d expect random noise to make something this signal-like less than once every 27000 years). However, we realised that we wanted to be able to show that Virgo had indeed seen something, and the search wasn’t set up for this.

Therefore, for the paper, we list three different checks to show that Virgo did indeed see the signal.

1. In a similar spirit to the main searches, we took the best fitting template (it doesn’t matter in terms of results if this is the best matching template found by the search algorithms, or the maximum likelihood waveform from parameter estimation), and compared this to a stretch of data. We then calculated the probability of seeing a peak in the signal-to-noise ratio (as shown in the top row of Figure 1) at least as large as identified for GW170814, within the time window expected for a real signal. Little blips of noise can cause peaks in the signal-to-noise ratio, for example, there’s a glitch about 50 ms after GW170814 which shows up. We find that there’s a 0.3% probability of getting a signal-to-ratio peak as large as GW170814. That’s pretty solid evidence for Virgo having seen the signal, but perhaps not overwhelming.
2. Binary black hole coalescences can also be detected (if the signals are short) by our searches for unmodelled signals. This was the case for GW170814. These searches were using data from all three detectors, so we can compare results with and without Virgo. Using just the two LIGO detectors, we calculate a false alarm rate of 1 per 300 years. This is good enough to claim a detection. Adding in Virgo, the false alarm rate drops to 1 per 5900 years! We see adding in Virgo improves the significance by almost a factor of 20.
3. Using our parameter estimation analysis, we calculate the evidence (marginal likelihood) for (i) there being a coherent signal in Livingston and Hanford, and Gaussian noise in Virgo, and (ii) there being a coherent signal in all three detectors. We then take the ratio to calculate the Bayes factor. We find that a coherent signal in all three detectors is preferred by a factor of over 1600. This is a variant of a test proposed in Veitch & Vecchio (2010); it could be fooled if the noise in Virgo is non-Gaussian (if there is a glitch), but together with the above we think that the simplest explanation for Virgo’s data is that there is a signal.

In conclusion: Virgo works. Probably.

#### Follow-up observations

Adding Virgo to the network greatly improves localization of the source, which is a huge advantage when searching for counterparts. For a binary black hole, as we have here, we don’t expect a counterpart (which would make finding one even more exciting). So far, no counterpart has been reported.

i

#### Announcement

This is the first observation we’ve announced before being published. The draft made public at time at announcement was accepted, pending fixing up some minor points raised by the referees (who were fantastically quick in reporting back). I guess that binary black holes are now familiar enough that we are on solid ground claiming them. I’d be interested to know if people think that it would be good if we didn’t always wait for the rubber stamp of peer review, or whether they would prefer to for detections to be externally vetted? Sharing papers before publication would mean that we get more chance for feedback from the community, which is would be good, but perhaps the Collaboration should be seen to do things properly?

One reason that the draft paper is being shared early is because of an opportunity to present to the G7 Science Ministers Meeting in Italy. I think any excuse to remind politicians that international collaboration is a good thing™ is worth taking. Although I would have liked the paper to be a little more polished [bonus advice]. The opportunity to present here only popped up recently, which is one reason why things aren’t as perfect as usual.

I also suspect that Virgo were keen to demonstrate that they had detected something prior to any Nobel Prize announcement. There’s a big difference between stories being written about LIGO and Virgo’s discoveries, and having as an afterthought that Virgo also ran in August.

The main reason, however, was to get this paper out before the announcement of GW170817. The identification of GW170817’s counterpart relied on us being able to localize the source. In that case, there wasn’t a clear signal in Virgo (the lack of a signal tells us the source wan’t in a direction wasn’t particularly sensitive to). People agreed that we really need to demonstrate that Virgo can detect gravitational waves in order to be convincing that not seeing a signal is useful information. We needed to demonstrate that Virgo does work so that our case for GW170817 was watertight and bulletproof (it’s important to be prepared).