Higher order gravitational-wave modes with likelihood reweighting

TL;DR

This paper presents an efficient two-step Bayesian method to incorporate higher-order gravitational-wave modes into astrophysical inference, improving parameter estimation without high computational costs, demonstrated on GWTC-1 data.

Contribution

The authors introduce a novel likelihood reweighting approach that allows inclusion of higher-order modes in gravitational-wave analysis efficiently.

Findings

Higher-order modes provide additional information but are computationally challenging to include.

The proposed method successfully incorporates higher modes into Bayesian inference.

Weak preference for higher modes was found in GWTC-1 events.

Abstract

The gravitational waveform of a merging stellar-mass binary is described at leading order by a quadrupolar mode. However, the complete waveform includes higher-order modes, which encode valuable information not accessible from the leading-order mode alone. Despite this, the majority of astrophysical inferences so far obtained with observations of gravitational waves employ only the leading order mode because calculations with higher-order modes are often computationally challenging. We show how to efficiently incorporate higher-order modes into astrophysical inference calculations with a two step procedure. First, we carry out Bayesian parameter estimation using a computationally cheap leading-order-mode waveform, which provides an initial estimate of binary parameters. Second, we weight the initial estimate using higher-order mode waveforms in order to fold in the extra information from the full waveform. We use mock data to demonstrate the effectiveness of this method. We apply the method to each binary black hole event in the first gravitational-wave transient catalog GWTC-1 to obtain posterior distributions and Bayesian evidence with higher-order modes. Performing Bayesian model selection on the events in GWTC-1, we find only a weak preference for waveforms with higher order modes. We discuss how this method can be generalized to a variety of other applications.

Tables

Table 1. Table 1: The log Bayes factor and the “efficiency,” equal to the number of effective samples n eff subscript 𝑛 eff n_{\text{eff}} divided by the initial number of samples n 𝑛 n .

event	$\ln \text{BF}$	$n_{\text{eff}}/n$
simulated	9.01	$2.2\times {10}^{-3}$
GW150914	-0.20	0.31
GW151012	0.40	0.36
GW151226	0.01	0.82
GW170104	0.01	0.57
GW170608	-0.16	0.73
GW170729	1.23	0.13
GW170809	-0.09	0.67
GW170814	0.19	0.48
GW170818	0.31	0.74
GW170823	-0.21	0.69
GWTC-1	1.49	N/A

Equations

$$p(\theta |d)=\frac{\mathcal{L}(d|\theta )\pi (\theta )}{\mathcal{Z}},$$

$$\mathcal{Z}=\int d\theta \mathcal{L}(d|\theta )\pi (\theta ),$$

$$p(\theta |d)=$$

$$=$$

$$w(d|\theta )\equiv \frac{\mathcal{L}(d|\theta )}{{\mathcal{L}}_{\text{O}}(d|\theta )},$$

$$\mathcal{Z}=$$

$$=$$

$$p_{\text{O}}(\theta |d)=\frac{{\mathcal{L}}_{\text{O}}(d|\theta )\pi (\theta )}{{\mathcal{Z}}_{\text{O}}},$$

$$\mathcal{O}\equiv \underset{t,\varphi }{\max }\frac{⟨h_{\text{O}}^{+},h^{+}⟩+⟨h_{\text{O}}^{\times },h^{\times }⟩}{\sqrt{\left(⟨h_{\text{O}}^{+},h_{\text{O}}^{+}⟩+⟨h_{\text{O}}^{\times },h_{\text{O}}^{\times }⟩\right)\left(⟨h^{+},h^{+}⟩+⟨h^{\times },h^{\times }⟩\right)}}.$$

$$n_{\text{eff}}=\frac{{\left({\sum }_k{w}_k\right)}^2}{{\sum }_k{w}_k^2}.$$

Full text

I Introduction

The precise morphology of a gravitational waveform encodes a wealth of information about the binary that produced it. Merging stellar-mass binaries are typically characterized by fifteen parameters: two mass parameters, six spin parameters, and seven extrinsic parameters, which describe the location and orientation of the binary with respect to the detector 111Binary neutron stars can be characterized by additional tidal parameters and all compact binaries can be characterized with additional eccentricity parameters, depending on their formation scenario.. Extracting binary parameters from gravitational-wave measurements enables tremendous science including sky maps for electromagnetic follow-up Abbott et al. (2019, 2017a), measurement of the neutron star equation of state Abbott et al. (2018a), measurement of cosmological parameters Schutz (1986); Abbott et al. (2017b), probing the fate of massive stars Fishbach and Holz (2017); Talbot and Thrane (2018); Abbott et al. (2018b), understanding the formation mechanisms of compact binaries Vitale et al. (2017); Stevenson et al. (2017); Talbot and Thrane (2017); Gerosa and Berti (2017); Farr et al. (2017); Wysocki et al. (2018); Lower et al. (2018), and testing general relativity Abbott et al. (2016, 2017c).

The parameters of compact binaries are estimated using Bayesian inference software Veitch et al. (2015); Biwer et al. (2019); Ashton et al. (2019); Pankow et al. (2015); Lange et al. (2018). The software employs nested sampling Skilling (2004), Markov Chain Monte Carlo Metropolis et al. (1953); Hastings (1970); Hogg and Foreman-Mackey (2018), or adaptive mesh refinement Pankow et al. (2015); Lange et al. (2018) in order to construct posterior distributions for binary parameters and/or to calculate the Bayesian evidence. Bayesian inference calculations in gravitational-wave astronomy are computationally demanding, and so significant research has been carried out in order to bring down the wall time of calculations, thereby enabling new science; see, e.g., Smith et al. (2016); Pürrer (2014); Canizares et al. (2013); Talbot et al. (2019); Wysocki et al. (2019).

The computational demands of inference have created a premium for fast approximate gravitational waveforms or “approximants” Hannam et al. (2014); Khan et al. (2016a). Fast approximants have enabled breakthrough science. However, the speed can come at a cost. The current approximants most commonly used in gravitational-wave inference are constructed using only the leading order, $\ell=2$ modes in the spin-weighted spherical harmonic decomposition, although see London et al. (2018). While these leading-order approximants provide reasonably good estimates of binary parameters, they do not incorporate all of the information in a gravitational waveform, and therefore provide an incomplete picture. Inference with higher-order modes can provide tighter constraints than those obtained with leading-order waveforms alone. In particular, higher-order modes are useful breaking degeneracy between binary parameters. For example, the $\ell=|m|=2$ waveform is totally invariant under a transformation in which the polarization angle and phase of coalescence advance by $\pi/2$. The ability to break this degeneracy is key to detecting gravitational-wave memory Lasky et al. (2016). This is just one example highlighting the scientific potential of inference with higher-order modes.

Astrophysical inference with higher-order mode waveforms was first demonstrated in Kumar et al. (2018); Lange et al. (2018), which applied a numerical relativity surrogate model Blackman et al. (2017) to produce posterior distributions for GW150914, GW170104, GW170608, and GW170814. More recently, adaptive mesh methods have been employed in order to derive posterior distributions and Bayesian evidence for GW170729 and other events in GWTC-1 Abbott et al. (2018c) using a variety of higher-order-mode waveforms Chatziioannou et al. (2019).

In this paper we demonstrate a fast and effective method to calculate posterior distributions and Bayesian evidence for gravitational-wave signals with higher-order-modes using a technique known as “importance sampling”; see, e.g., Robert and Casella (2004); Liu (2004). First, we carry out Bayesian parameter estimation using a low-cost, $\ell=|m|=2$ waveform, which yields an approximate answer on which we can improve. In the second step, we calculate a weight factor for each posterior sample, which incorporates information from the more expensive higher-order-mode waveform. Using the weights, we obtain the posterior and evidence, which we would have obtained if we had carried out the entire calculation using the more expensive higher-order mode waveform.

The remainder of this paper is organized as follows. In Section II, we describe the reweighting formalism. Then, in Section III, we show results obtained for both simulated data and for events in the LIGO/Virgo catalog, GWTC-1 Abbott et al. (2018c); Vallisneri et al. (2014). Posterior samples and Bayesian evidence for every event in GWTC-1 are available on a companion web page Payne et al. (2019a), along with the code used in our analysis. We provide a table summarizing the evidence of higher-order modes in GWTC-1. We show that, while higher-order modes produce tighter constraints than $\ell=2$ analyses, there is not yet a strong signature of higher-order modes in published LIGO/Virgo detections. In Section IV, we discuss other possible directions for future research including novel applications of this reweighting technique, which may be useful for a variety of problems in astrophysics. Finally, we include an appendix, Section B, with the rest of the results from GWTC-1.

II Methodology

We show how to carry out Bayesian inference using an approximate “proposal” likelihood ${\cal L}_{\O}(d|\theta)$ to obtain initial posterior samples, which are then reweighted using a more computationally expensive “target” likelihood ${\cal L}(d|\theta)$. Here, $d$ is the data and $\theta$ represents the parameters. The proposal likelihood is an approximation for the target likelihood. In order for reweighting to be efficient, the proposal likelihood should be similar to the target likelihood, so that the two likelihoods overlap significantly. For demonstration purposes, we use as our proposal waveform IMRPhenomD Khan et al. (2016b), an aligned-spin, $\ell=|m|=2$ approximant, which is widely used in astrophysical inference thanks to its reliability and speed. For our target waveform, we use NRHybSur3dq8 Varma et al. (2019), a numerical relativity surrogate model, which includes higher-order modes up to $\ell=4$ excepting $(4,\pm 1)$ and $(4,0)$, but including $(5,\pm 5)$, aligned spin, and mass ratios $m_{2}/m_{1}\geq 0.125$.

Our goal is to derive expressions for the target posterior

[TABLE]

and the target Bayesian evidence

[TABLE]

written in terms of a fast-to-calculate, proposal likelihood.

The proposal quantities are linked to the target quantities by a weight factor. Multiplying the target posterior by unity, we obtain

[TABLE]

Here,

[TABLE]

is the weight function. Multiplying by unity again, we obtain the following expression for the evidence

[TABLE]

The second line replaces the integral with a discrete sum over $n$ proposal posterior samples; see Thrane and Talbot (2019).

Carrying out Bayesian inference with the proposal likelihood, we obtain “proposal posterior samples” for the distribution

[TABLE]

where ${\cal Z_{\O}}$ is the proposal evidence. We generate our proposal samples using the Bilby Ashton et al. (2019) implementation of CPNest Pozzo and Veitch (2019). We compute the weights using likelihoods marginalized over the coalescence time and reference phase (see, e.g., Thrane and Talbot (2019)) in order to avoid differences in the definition of these parameters between waveform models.

In a previous iteration of this work Payne et al. (2019b), we employed a slightly different approach. Instead of marginalizing over time and phase, we chose the values of these parameters that maximized the overlap ${\cal O}$ of the proposal waveform and target waveform:

[TABLE]

Here, $h_{\O}^{+,\times}$ are the plus and cross components of the proposal waveform while $h^{+,\times}$ are the plus and cross components of the target waveform. The angled brackets denote noise-weighted inner products. This approach, used subsequently in other publications, e.g., Romero-Shaw et al. (2019), produces qualitatively similar results to marginalizing over ($t$,$\phi$), albeit with reduced computational efficiency. Moreover, the marginalization approach is better justified from a statistical point of view.

Weighting each sample by $w(d|\theta)$, and renormalizing, we convert the proposal posterior samples into target posterior samples. The application of weight factors has the effect of reducing the effective number of samples Kish (1995); Víctor Elvira (2018)

[TABLE]

After reweighting, it is therefore prudent to calculate $n_{\text{eff}}$ in order to determine that there are a suitably large number of samples. It is straightforward to generate more weighted posterior samples by simply combining the results from multiple proposal analyses run in parallel.

The method of likelihood reweighting outlined here is similar to the procedure of “recycling” commonly used to study the population properties of compact objects; see, e.g., Abbott et al. (2018b); Thrane and Talbot (2019). Previous applications of recycling have, in effect, carried out reweighting to change the prior in post-processing. The principle here is the same, except we change the likelihood. Our formalism can be straightforwardly extended to simultaneously alter the prior (informed by astrophysics) and likelihood (using more sophisticated waveforms). Reweighting is an application of importance sampling. For additional details, readers are referred to Robert and Casella (2004); Liu (2004).

III Results

We demonstrate likelihood reweighting using a simulated binary black hole merger signal injected into Gaussian noise. We assume a two-detector LIGO network operating at design sensitivity Abadie et al. (2015). Using NRHybSur3dq8, we inject a binary black hole waveform. The binary, located at a luminosity distance $d_{L}=400\,\mathrm{Mpc}$, has chirp mass $\mathcal{M}=30M_{\odot}$ and mass ratio $q=m_{2}/m_{1}=0.8$. The dimensionless aligned spins are $\chi_{1}=0.4,\chi_{2}=0.3$. The signal has a network matched-filter signal-to-noise ratio of $\rho_{\text{mf}}=55$. In Fig. 1 we provide a corner plot showing the posterior distribution and credible intervals obtained for this simulated event. The blue shades indicate the posterior derived using our proposal IMRPhenomD waveform while the green shades indicates the posterior obtained after reweighting with our target NRHybSur3dq8 waveform. The darkness of the contours indicate credible intervals at $1\sigma,2\sigma,3\sigma$. The true values of each parameter are indicated by orange markers.

While both posteriors includes the true parameter values, the blue IMRPhenomD posterior is broad in comparison to the green NRHybSur3dq8 posterior. For some parameters, the posterior shrinks dramatically when we add information from higher-order modes. Higher-order modes break the degeneracy between distance and inclination while improving our ability to measure the mass ratio. This, in turn, improves our estimation of the spins.

We calculate the Bayes factor comparing the hypothesis that the data are best fit by NRHybSur3dq8 to the hypothesis that they are best fit by IMRPhenomD. This is a measure of degree to which the data prefer a model including higher-order modes 222To be clear, we are not proposing $\ell|m|=22$ gravitational waves as a serious alternative to general relativity. Rather, we are using this Bayes factor in order to ascertain if we can resolve the higher-order mode content predicted by general relativity.. The log Bayes factor for our simulated event are reported in the first row of Tab. 1. We include also the “efficiency,” the number of effective samples used in each calculation, normalized by the number of proposal samples.

For the simulated event, the signal-to-noise ratio is sufficiently large to “detect” the presence of higher-order modes with a high significance $\ln\text{BF}=9.01$. Since this simulated event has a high signal-to-noise ratio, the ratio of the effective number of samples to the number of proposal samples is small, $\approx 0.2\%$. Thus, a large number of proposal runs ($\approx 1000$) are required in order to produce a well-sampled target posterior. Fortunately, these proposal runs are embarrassingly parallel, which means that the wall time is no longer than a single proposal inference run, provided sufficient computational resources are available. It should also be noted that we are able to get posterior samples and evidence for all of the events in GWTC-1 with just a few parallel runs. We further note that it may be possible to improve efficiency by restricting the proposal prior based on early returns from the target posterior. The events in GWTC-1 have lower signal-to-noise ratio, and so the reweighting procedure is much more efficient.

We now turn toward real data in the GWTC-1 catalog. For each event, we produce: a set of posterior samples with weights, the proposal Bayesian evidence (obtained with IMRPhenomD), the target Bayesian evidence (obtained with NRHybSur3dq8), the NRHybSur3dq8/IMRPhenomD Bayes factor, and corner plots with credible intervals. The full results are available here: Payne et al. (2019a). Key summary statistics are provided in Tab. 1. For illustrative purposes, we also include in Fig. 2 the corner plot for GW170729, the event with the greatest support for higher-order modes ($\ln\text{BF}=1.23$). Our target posterior for GW170729 is qualitatively similar to previous results from Chatziioannou et al. (2019). We observe increased support for non-unity mass ratio $q=m_{2}/m_{1}$, slightly more support for zero-spin, and changes to the posterior distributions of the extrinsic parameters. An additional validation of the results for GW170729 are included in Appendix A.

IV Discussion

As the gravitational-wave catalog grows, higher-order modes will provide a more precise view of the properties of compact objects. Higher-order modes become increasingly important for systems: with large mass ratios $q\lesssim 0.5$, with significant spin, and with large total mass Varma and Ajith (2017a). The method described here will be particularly useful for exploring such events. In order to make use of likelihood reweighting, it is necessary to employ a suitable proposal model, capable of producing posterior samples, which provide a reasonable approximation of the true distribution. For inference with higher-order modes, this appears to be achievable using quadrupolar approximants.

One must exercise a degree of caution when choosing the proposal likelihood to make sure it is suitable for the problem at hand. For example, binary mergers with very low mass ratios $q\lesssim 1/6$ (not yet observed by LIGO/Virgo) can produce so much higher-order mode content that the quadrupolar approximants yield biased parameter estimates Calderón Bustillo et al. (2016); Varma and Ajith (2017b). There are two telltale signs, which indicate that the proposal likelihood is poorly suited to the problem. First, the efficiency $n_{\text{eff}}/n$ will be low even for modest signal-to-noise ratio because the proposal samples are all in the tail of the target likelihood distribution. The sample(s) closest to the target likelihood peak will be weighted as far more important than the other samples. Second, the reweighted posterior distributions will be peak at the edges of the sampled region of parameter space. By comparing the target posterior to the proposal, it is straightforward to determine if the target posterior is railing against the edge of the proposal samples 333Of course, if the target posterior is multi-modal, there is an additional risk that the proposal posterior might miss a mode, which would not be apparent through visual inspection. We have not observed such multi-modal behavior in the examples considered in this manuscript, but one can guard against this failure mode through $pp$ tests as described in the appendix.. If these two tests indicate a poorly chosen proposal likelihood, one should devise a better one.

We foresee a number of useful applications for likelihood reweighting including: inferences about gravitational-wave memory Lasky et al. (2016), noise models that add complexity beyond the usual Gaussian assumption, treatment of calibration errors, and inference with computationally expensive waveforms that include tidal effects or eccentricity 444After this paper appeared on the arxiv, our method was used in Romero-Shaw et al. (2019)..

V Acknowledgements

We thank Rory Smith and Greg Ashton for sharing a beta version of their parallelized nested sampling algorithm, used for benchmarking in the appendix. We thank Rory Smith for helping to debug our parallel inference runs. We thank Katerina Chatziioannou, Richard O’Shaughnessy, and Vijay Varma for helpful comments. We thank Moritz Hübner, and the Bilby team for support. This is document LIGO-P1900128. EP, ET, and CT are supported by ARC CE170100004. ET is supported by ARC FT150100281. This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. LIGO is funded by the U.S. National Science Foundation. Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by Polish and Hungarian institutes.

Appendix A Validation

In order to test the reweighting procedure, we calculate the posterior for GW170729 using direct sampling and compare the results obtained with reweighting. The direct sampling method is computationally expensive. However, we are able to make use of a recent work by Smith & Ashton Smith and Ashton (2019), who have created a parallelized implementation of the nested sampler, dynesty Speagle (2019). The calculation ran for 25 hours on 830 cores, which suggests it would have taken more than two years to complete on a single core. This is in contrast to our calculation, which took approximately 5 hours on 25 cores to generate the target samples and 9 hours on 450 cores to calculate weights (66 core days in total). In Fig. 3, we plot the posterior for GW170729 obtained using reweighting (green) alongside the posterior obtained using the parallel inference code from Smith et al. (2016) (purple). The results are in approximate agreement, with small differences that are typical of the differences that arise when a posterior is calculated with two different implementations of nested sampling, for example, dynesty and CPNest Pozzo and Veitch (2019). It is difficult to ascertain, which result is a more accurate representation of the true posterior distribution, but the level of disagreement should be taken as indicative of the systematic error from reweighting.

In Fig. 4, we include a “$pp$ plot,” which provides another means of validating the reweighting method. The $pp$ plot shows the fraction of 100 simulated events (vertical axis), for which the true parameter value falls within a given confidence interval (horizontal axis). Ideally, the colored traces (representing different astrophysical parameters) should be consistent with a slope=1 diagonal line. The gray region indicates the typical fluctuations (90% credible interval) for a single parameter due to the finite sample size. The distribution of colored traces about the slope=1 line indicates that the reweighting method is producing reliable posteriors to within the precision we can test.

Appendix B Results from GWTC-1

In this appendix, we present posterior distributions and credible intervals for the rest of the events in GWTC-1. In each figure, the blue distributions show the posteriors obtained using the $\ell=|m|=2$ approximant IMRPhenomD Khan et al. (2016b). The green distribution shows the posteriors obtained using the approximant NRHybSur3dq8 Varma et al. (2019), which includes higher-order modes $\ell\leq 4$ excepting $(4,\pm 1)$ and $(4,0)$, but including $(5,\pm 5)$. The intrinsic binary parameters (total mass, mass ratio, and effective aligned spin) are on the left. The extrinsic parameters, luminosity distance and inclination angle, are on the right.