Testing general relativity using gravitational wave signals from the inspiral, merger and ringdown of binary black holes

TL;DR

This paper explores how gravitational wave observations from binary black hole mergers can be used to test general relativity, especially by combining multiple modest SNR events to detect potential deviations.

Contribution

It introduces a method to combine data from many low SNR gravitational wave events to improve constraints on deviations from general relativity.

Findings

Combining multiple events enhances sensitivity to deviations from GR.

The method can identify certain types of deviations if present in the waveforms.

The robustness of the test is confirmed against analysis parameter variations.

Abstract

Advanced LIGO's recent observations of gravitational waves (GWs) from merging binary black holes have opened up a unique laboratory to test general relativity (GR) in the highly relativistic regime. One of the tests used to establish the consistency of the first LIGO event with a binary black hole merger predicted by GR was the inspiral-merger-ringdown consistency test. This involves inferring the mass and spin of the remnant black hole from the inspiral (low-frequency) part of the observed signal and checking for the consistency of the inferred parameters with the same estimated from the post-inspiral (high-frequency) part of the signal. Based on the observed rate of binary black hole mergers, we expect the advanced GW observatories to observe hundreds of binary black hole mergers every year when operating at their design sensitivities, most of them with modest signal to noise ratios (SNRs). Anticipating such observations, this paper shows how constraints from a large number of events with modest SNRs can be combined to produce strong constraints on deviations from GR. Using kludge modified GR waveforms, we demonstrate how this test could identify certain types of deviations from GR if such deviations are present in the signal waveforms. We also study the robustness of this test against reasonable variations of a variety of different analysis parameters.

Tables

Table 1. Table 1: The table shows the parameters of some of the injected modified GR signals ( D L subscript 𝐷 L D_{\mathrm{L}} is the luminosity distance) and how well the PyCBC pipeline is able to detect these injected signals in a single detector search (LIGO Livingston). Optimal SNR is the optimal SNR available in the injection, PyCBC SNR is the matched filter SNR extracted by the PyCBC pipeline, and re-weighted SNR is the chi-squared weighted SNR defined in equation (6) of  [ 52 ] . The optimal SNR is computed using the detector’s noise power spectrum (it is defined by ⟨ h | h ⟩ 1 / 2 superscript inner-product ℎ ℎ 1 2 \langle h|h\rangle^{1/2} , for a signal h ℎ h ), while the matched filter SNR is a statistical variable computed using a specific noise realization, which is why the PyCBC SNR can be larger than the optimal SNR.

$m_1\left(M_{\odot }\right)$	$m_2\left(M_{\odot }\right)$	$D_L$ (Mpc)	Optimal SNR	PyCBC SNR	Reduced ${\chi }^2$	Re-weighted SNR
139.1	69.5	5625.8	8.1	8.6	0.76	8.6
88.1	14.7	803.5	14.3	14.0	0.78	14.0
40.8	46.8	787.0	21.6	21.7	0.81	21.7

Equations

$$P(\mathbf{\lambda }|d,\mathcal{H},I)=\frac{P(\mathbf{\lambda }|\mathcal{H},I)\mathcal{L}(d|\mathbf{\lambda },\mathcal{H},I)}{P(d|\mathcal{H},I)},$$

$$\mathcal{L}(d|\boldsymbol{\lambda},\mathcal{H},I)\propto\exp\big{[}-\frac{1}{2}\langle d-h(\boldsymbol{\lambda})\,|\,d-h(\boldsymbol{\lambda})\rangle\big{]},$$

$$⟨B|C⟩:=2{\int }_{f_{\mathrm{low}}}^{f_{\mathrm{cut}}}\frac{{\tilde{B}}^{\ast }(f)\tilde{C}(f)+\tilde{B}(f){\tilde{C}}^{\ast }(f)}{S_n(f)}\mathrm{d}f.$$

$$M_f=M_f(m_1,m_2,{\mathbf{a}}_1,{\mathbf{a}}_2),a_f=a_f(m_1,m_2,{\mathbf{a}}_1,{\mathbf{a}}_2),$$

$$ϵ:=\frac{\Delta M_f}{{\bar{M}}_f},\sigma :=\frac{\Delta a_f}{{\bar{a}}_f},$$

$$\Delta M_f:=M_f^{\text{textsc}i}-M_f^{\text{textsc}mr},\Delta a_f:=a_f^{\text{textsc}i}-a_f^{\text{textsc}mr},$$

$${\bar{M}}_f:=\frac{M_f^{\text{textsc}i}+M_f^{\text{textsc}mr}}{2},{\bar{a}}_f:=\frac{a_f^{\text{textsc}i}+a_f^{\text{textsc}mr}}{2}.$$

$$P(ϵ,\sigma |\{d_j\})=P(ϵ,\sigma )\prod _{j=1}^N\frac{P_j(ϵ,\sigma |d_j)}{P_j(ϵ,\sigma )},$$

$$\mathrm{credible}\mathrm{level}={\int }_{\{ϵ,\sigma \}\mathrm{where}P(ϵ,\sigma )>P({ϵ}_0,{\sigma }_0)}P(ϵ,\sigma )\mathrm{d}ϵ\mathrm{d}\sigma .$$

$$|\tilde{h}(f)|=\frac{A(t_f)}{\sqrt{\dot{F}(t_f)}}$$

$$M_f^{\text{textsc}i}$$

$$M_f^{\text{textsc}mr}$$

$$P(ϵ,\sigma |d)$$

Full text

\MakePerPage

footnote

1 Introduction

On September 14, 2015 the Advanced LIGO detectors [1, 2] made the first observation of gravitational waves (GWs) from a binary black hole system [3]. The signal, termed GW150914, was inferred to be produced by the merger of two black holes (with masses $\sim 36{M}_{\odot}$ and $29{M}_{\odot}$) at a distance of $\sim 400~{}\mathrm{Mpc}$ from the Earth [4, 5]. This observation was followed by a second one on December 26, 2015 [6] and a third one on January 4, 2017 [7]. The second signal, termed GW151226, was produced by the coalescence of two black holes with masses $\sim 14{M}_{\odot}$ and $8{M}_{\odot}$, while the third signal GW170104 was produced by a binary black hole system more like GW150914’s source (masses $\sim 31{M}_{\odot}$ and $19{M}_{\odot}$). These observations herald the beginning of a new branch of observational astronomy. Apart from providing the first direct evidence of the existence of GWs, these observations confirmed the existence of stellar mass black holes that are much more massive than previously thought ($\gtrsim 25{M}_{\odot}$) [8]. They also provided the first evidence of binary black holes that inspiral under GW emission and subsequently merge. These observations also enabled us to perform the first tests of general relativity (GR) in the highly relativistic and nonlinear regime of gravity [9, 5] — a regime with binary orbital velocities $\sim 0.5\,c$, which is inaccessible by other astronomical observations and laboratory tests.

A handful of tests have been performed to check the consistency of the observed signals with those predicted by GR. In particular, due to the large masses of the black holes in GW150914, the observed signal consists of signatures from the inspiral and merger of the two black holes and the subsequent ringdown of the final black hole. This allowed us to perform several consistency tests making use of one or more phases of the coalescence: The first involved testing the consistency between the mass and spin of the final remnant, determined from the low-frequency (inspiral) and high-frequency (post-inspiral) parts of the observed signal [10]. The second involved testing the consistency of the data after the peak of the observed signal (corresponding to the merger) with a quasi-normal mode (QNM) spectrum predicted by GR [11, 12, 13, 14]. The third involved bounding deviations from the GR predictions of the post-Newtonian (PN) coefficients describing the inspiral (and from phenomenological parameters describing the merger and ringdown) using parametrized waveform models [15, 16, 17]. The fourth involved constraining the amount of dispersion in the observed GW signal and converting it to a bound on the mass of the graviton [18]. In addition, the residuals after subtracting the best fit GR signal have been found to be consistent with detector noise. Within the statistical uncertainties, these investigations provided no evidence for deviations from GR [9]. Some of these tests were repeated by using the second and third LIGO events, and combining the results from multiple events allowed us to improve the constraints on certain departures from GR predictions [5, 7]. Of particular interest to this paper is the first test mentioned above that checks for the consistency of the properties of the remnant estimated from the inspiral and merger–ringdown parts of the observed signal, which we call the inspiral–merger–ringdown (IMR) consistency test [10]. This paper provides further details of the formulation and implementation of this test that was used to constrain certain departures from GR using the LIGO events GW150914 and GW170104 [7].

The Advanced LIGO detectors, are not yet operating at than their design sensitivities [1], but are expected to approach design sensitivity over the next few years [19]. The Advanced Virgo [20] detector in Italy joined the Advanced LIGO detectors towards the end of the second observing run, and commissioning should bring it to design sensitivity on a similar timescale as the Advanced LIGO detectors. Several observing runs will be conducted over this period. Based on the rate of GW events so far [7], and the anticipated improvement in the sensitivity of the instruments, we expect that GW observations will become routine. The binary black hole event rate is expected to approach several hundred detections per year when the LIGO detectors operate at their design sensitivity [21, 5]. The KAGRA detector [22] in Japan is expected to join the international network of advanced GW detectors in the next few years, and the LIGO-India [23] detector is planned to join the network around 2024 [19]. These additional detectors will dramatically improve the sky coverage, source localization accuracy and the polarization extraction capability of the international GW network.

Upcoming GW observations will allow us to make significant improvements to the present constraints on deviations from GR, or potentially to detect such deviations — either from a small number of particularly loud events or by combining constraints from a large number of moderate events. Thus, a careful characterization of the analysis pipelines and various systematic errors in the analysis is of vital importance, so that systematic errors in our analysis will not be mistaken for a true deviation from GR. In anticipation of the binary black hole events expected in the upcoming observational runs of Advanced LIGO and other ground-based observatories, we present here a detailed follow up of the analysis method of the IMR consistency test for binary black holes presented in our previous paper [10], which we will henceforth refer to as Paper I.

The IMR consistency test was developed based on the original idea by Hughes and Menou [24] proposed in the context of the space-based LISA observatory.111See also, Nakano, Tanaka and Nakamura [25] for a recent study in the context of a ground-based detector. The key idea is to infer the mass and spin of the remnant black hole, using two different parts of the observed signal, and then to compare these independent estimates. We first estimate the initial masses and spins from the inspiral (low-frequency) part of the signal, which allows us to infer the mass and spin of the final black hole making use of fitting formulas calibrated to numerical relativity (NR) simulations of binary black holes. Next, we estimate the same parameters independently from the merger–ringdown (high-frequency) part of the signal and then compare the two estimates. If the signal is correctly described by the merger of a (quasicircular) binary of Kerr black holes in GR, which is implicit in the waveform models that we use to estimate the parameters and the fits that went in, one should expect the two estimates to be consistent with each other. On the other hand, if there is a departure from GR, depending on the exact nature of the departure, it can manifest as a discrepancy between the two estimates. In particular, if the energy and angular momentum radiated during the merger regime (where the gravity is extremely strong and nonlinear) differ significantly from the GR predictions for these quantities, one could expect a discrepancy between the two estimates in the final mass and spin.

The original idea by Hughes and Menou [24] was to estimate the parameters purely from the early inspiral (well described by the PN approximation to GR) and the late ringdown (well described by a spectrum of QNMs). According to this, we would estimate the parameters of the binary from the early inspiral signal (where any departure from GR is presumably small) and use these estimates to predict the properties of the remnant black hole, which could be estimated independently from the QNM ringdown. However, such a test will be possible only using a small number of golden events, where both the early inspiral and late ringdown are observed with high SNRs. While such a test might be possible in the future using LISA or third generation ground-based instruments, such tests are unlikely to be possible using the current generation of GW observatories. Meanwhile, Advanced LIGO is expecting to observe several hundred binary black hole events of moderate SNR in the coming years. Hence, tests of GR in the near future are going to progress through building evidence by combining multiple events. Our formulation of the test is geared in this direction. Indeed, the transition from inspiral regime to merger regime does not happen at a precise time or frequency. We use some reasonable choice of a cutoff frequency to separate the signal to the low-frequency (inspiral) and high-frequency (merger–ringdown) parts. As part of our studies, we show that the test is robust against variations of this cutoff frequency that still give large enough SNRs in for parameter estimation in both the low- and high-frequency portions independently.

Section 2 contains a detailed description of the analysis method. We then show how a large number of GW signals from binary black holes can be combined to produce strong constraints on the deviations from GR (section 3), and how a true deviation from GR might show up in our analysis (section 4). We also demonstrate the robustness of our analysis against reasonable choices of the frequency that is used to demarcate the inspiral and post-inspiral parts of the observed signal (section 5.1), the choice of the waveform template (section 5.2), fitting formulas for the mass and spin of the final black hole (section 5.3), and the effects of precession and higher modes (section 5.4). We conclude and consider future work in section 6.

2 Testing the consistency between the inspiral, merger and ringdown

The GW signal from a spinning binary black hole with negligible orbital eccentricity is described by $15$ parameters (enumerated in, e.g., section II B of [26]). If the spins of the individual black holes are aligned (or antialigned) with the direction of the orbital angular momentum, then the number of parameters reduces to $11$. These parameters are typically estimated by means of Bayesian inference. Bayes’ theorem states that the posterior probability distribution of a parameter set $\boldsymbol{\lambda}:=\{\lambda_{i}\}$ of a model hypothesis $\mathcal{H}$, given data $d$ and any other information $I$, is:

[TABLE]

where $P(\boldsymbol{\lambda}|\mathcal{H},I)$ is the prior probability of the parameter set $\boldsymbol{\lambda}$ given $\mathcal{H}$ and $I$, while $\mathcal{L}(d|\boldsymbol{\lambda},\mathcal{H},I)$ is called the likelihood function, which is the probability of observing the data $d$, given $\boldsymbol{\lambda}$, $\mathcal{H}$ and $I$. What appears in the denominator is a normalization constant $P(d|\mathcal{H},I):=\int P(\boldsymbol{\lambda}|\mathcal{H},I)\,\mathcal{L}(d|\boldsymbol{\lambda},\mathcal{H},I)\,\mathrm{d}\boldsymbol{\lambda}$, called the marginal likelihood, or the evidence of the hypothesis $\mathcal{H}$.

If our hypothesis $\mathcal{H}$ is that the data contains a GW signal described by a GR waveform model $h(\boldsymbol{\lambda})$ and stationary Gaussian noise described by the power spectral density $S_{n}(f)$, then, as described in appendix A of [27], the likelihood function can be defined as:

[TABLE]

where $\langle\cdot|\cdot\rangle$ is the following noise-weighted inner product

[TABLE]

Above, $\tilde{B}(f)$ denotes the Fourier transform of $B(t)$ and a ∗ denotes complex conjugation. The limits of integration ${f_{\mathrm{low}}}$ and ${f_{\mathrm{cut}}}$ are dictated by the bandwidth of the detector sensitivity, the bandwidth of the signal, as well as the cutoff frequencies in our calculations (described below). Owing to the large dimensionality of the parameter set $\boldsymbol{\lambda}$, the posterior distribution $P(\boldsymbol{\lambda}|d,\mathcal{H},I)$ in equation (1) is computed by stochastically sampling the parameter space using techniques such as Markov-Chain Monte Carlo [28] or nested sampling [29]. For this paper, we use the LALInferenceNest [30, 26] code that provides an implementation of the nested sampling algorithm in the LSC Algorithms Library [31] for computing the posterior distributions.

Using the Bayesian framework described above, we start by computing the joint posterior probability distribution on the initial masses and dimensionless spins $P(m_{1},m_{2},\boldsymbol{a}_{1},\boldsymbol{a}_{2}\,|\,d)$222Here onwards, we drop $\mathcal{H},I$ from the posteriors to simplifying the notation. Also, unless otherwise noted, the masses we consider are the ones in the detector frame, including the cosmological redshift [32]. of the binary black hole system (with $m_{1}\geq m_{2}$ by convention). This is done by marginalizing the posterior over the remaining parameters describing the signal. Assuming quasicircular inspirals, fitting formulas calibrated to NR simulations then give us predictions of the mass $M_{f}$ and dimensionless spin magnitude $a_{f}$ of the final (remnant) black hole as functions of the initial masses and spins,

[TABLE]

which allow us to obtain the posterior probability distribution $P(M_{f},a_{f}\,|\,d)$ of the final mass and spin.

First, we estimate the posterior $P_{\textsc{imr}}(M_{f},a_{f}\,|\,d)$ using the full observed signal. We choose to demarcate the inspiral and post-inspiral parts of the signal using the $m=2$ mode GW frequency $f_{\textsc{isco}}$ of the innermost stable circular orbit (ISCO) of the remnant Kerr black hole [33], with mass and spin given by the median values of $P_{\textsc{imr}}(M_{f},a_{f}\,|\,d)$. We estimate the same parameters from the data containing only the inspiral (low-frequency) part of the observed signal; that is, by setting $f_{\mathrm{cut}}=f_{\textsc{isco}}$ in equations (2) and (3). This allows us to compute the posterior distribution $P_{\textsc{i}}(M_{f},a_{f}\,|\,d)$ of the mass and spin of the remnant purely from the inspiral part of the signal. Similarly, using only the merger–ringdown (high-frequency) part of the observed signal [by setting $f_{\mathrm{low}}=f_{\textsc{isco}}$ in equations (2) and (3)], we can get yet another estimate of the posterior $P_{\textsc{mr}}(M_{f},a_{f}\,|\,d)$ of the mass and spin of the remnant. If the observed signal is well described by GR, the two independent estimates $P_{\textsc{i}}(M_{f},a_{f}\,|\,d)$ and $P_{\textsc{mr}}(M_{f},a_{f}\,|\,d)$ have to be consistent with each other as well as with the estimate $P_{\textsc{imr}}(M_{f},a_{f}\,|\,d)$ using the full data (see, e.g., the left plot in figure 1).

To quantify the consistency of the observed signal with a binary black hole system predicted by GR, we define two parameters that describe the fractional difference between the two estimates of the remnant’s mass and spin

[TABLE]

where

[TABLE]

and

[TABLE]

In the absence of departures from GR, we expect the posterior $P(\epsilon,\sigma\,|\,d)$ of $\{\epsilon,\sigma\}$ to be consistent with zero (within the expected statistical fluctuations due to noise).333We used the full IMR posteriors for $M_{f}$ and $a_{f}$ to normalize the fractional parameters in Paper I. We have switched to using the average of the inspiral and merger–ringdown posteriors here so that we only are using posteriors from the two independent measurements, and not mixing in ones from the case where the parameters describing the inspiral and merger–ringdown parts of the waveform are not allowed to vary independently. A describes the calculation of $P(\epsilon,\sigma\,|\,d)$ from the posteriors $P_{\textsc{i}}(M_{f},a_{f}\,|\,d)$ and $P_{\textsc{mr}}(M_{f},a_{f}\,|\,d)$. An example of the posterior distribution $P(\epsilon,\sigma\,|\,d)$ from a simulated GR signal is shown in figure 1 (right plot).

Finally, if we assume that $\{\epsilon,\sigma\}$ take the same values in multiple events, the posterior from one event could be treated as the prior for computing the posterior for the second event (see section 4.3 for a discussion on the validity of this assumption). Such a combined posterior can more sensitively identify certain deviations from GR. In practice, the posteriors $P_{j}(\epsilon,\sigma\,|\,d_{j})$ from $N$ such observations can be combined to produce an overall posterior in the following way:

[TABLE]

where $P(\epsilon,\sigma)$ is the overall prior distribution on $\epsilon$ and $\sigma$ and $P_{j}(\epsilon,\sigma)$ is the prior distribution used to compute the posterior $P_{j}(\epsilon,\sigma\,|\,d_{j})$.

2.1 Choice of priors

Given our poor understanding of the astrophysical distribution of source parameters (such as the component masses and spins), we use uninformative priors for the parameter estimation [26, 4]. Sources are assumed to be uniformly distributed in volume and isotropically oriented. We use uniform priors in component masses $m_{1,2}\in[1,300]{M}_{\odot}$, where $m_{1}\geq m_{2}$, and a uniform prior in the spin magnitudes $|\boldsymbol{a}_{1,2}|\in[0,0.98]$. We use isotropic priors on the spin orientation when precessing templates are used in the parameter estimation. When nonprecessing templates are used, the priors are set in such a way that the dimensionless spin projection $\boldsymbol{a}_{1,2}\cdot\hat{\mathbf{L}}$ onto the orbital angular momentum $\mathbf{L}$ has a uniform distribution in $[-0.98,0.98]$. These choices induce a non-uniform prior in $\{\epsilon,\sigma\}$, which is shown by the dashed contours in Figure 1 (right panel).

3 Simulations of GR signals in Gaussian noise

First, we estimate the expected constraints that we can place on parameters describing deviations from GR, when the actual signals are well described by GR. We inject simulated GW signals modeling inspiral, merger and ringdown of binary black holes (based on GR, modeled by the SEOBNRv2_ROM_DoubleSpin approximant [34, 35]) into colored Gaussian noise with the design power spectrum of the Advanced LIGO detectors in the high-power, zero-detuning configuration [36]. Although we perform the analysis assuming the $3$-detector Advanced LIGO–Virgo network, for simplicity, we assume that Advanced Virgo has the same noise power spectrum as the Advanced LIGO detectors. All detectors are assumed to have a low frequency cutoff $f_{\mathrm{low}}=10~{}\mathrm{Hz}$. We simulate a population of binary black holes that is uniformly distributed in comoving four-volume in the redshift interval $z=[0,2]$, with isotropic orientations (see Section II C of [37] for more details). Binaries have (source-frame) component masses uniformly distributed in the range $m^{\mathrm{s}}_{1,2}=[10,80]M_{\odot}$ and component spin magnitudes uniformly distributed in the range $[0,0.98]$. The GW signals are redshifted so that the effective masses describing the signal in the detector frame become $m_{1,2}=m^{\mathrm{s}}_{1,2}(1+z)$ [32]. Black hole spins are always assumed to be aligned/antialigned with the orbital angular momentum, as is required by the waveform model used: This is not an overly restrictive assumption, because the spin components along the orbital angular momentum have the dominant effect on the binary’s dynamics and GW emission [38]. Figure 2 shows the distribution of the network SNR of the injected population of binaries. For our analysis, we consider signals with network SNR greater than $8$ for both the inspiral and merger–ringdown parts, using the Kerr ISCO frequency to split the signal.444Here, for simplicity, we use the Kerr ISCO frequency computed from the injected values of masses and spins to split the signal. In an actual observation, we use the median estimate of the Kerr ISCO frequency estimated from the observed IMR signal. In addition, we only consider binaries with total redshifted mass less than $150M_{\odot}$, so that the observed signals contain imprints of the inspiral, merger and ringdown. Hence, starting from an initial population of $5000$ injections, the final analysis is performed on around $150$ events which survive the above thresholds on SNR and total redshifted mass.

On the injections that pass the SNR and mass thresholds, we perform the analysis described in section 2, using SEOBNRv2_ROM_DoubleSpin as our template for parameter estimation and the final mass and spin fits from Healy et al. [39]. Parameter estimation using LALInferenceNest is performed using priors that are uniform in component masses in the interval $m_{1,2}\in[1,300]M_{\odot}$ and component (aligned/antialigned) spins in the interval $[0,0.98]$. This corresponds to a non-uniform prior in the deviation parameters $\{\epsilon,\sigma\}$ that is centered around the GR value of $\{0,0\}$. See figure 1 for an example result. The left panel shows the estimated posterior distributions of the final mass and spin from the inspiral and post-inspiral parts of the data, along with that estimated from the full data. The right panel shows the posterior distribution on the parameters describing deviations from GR, along with the prior used for computing this posterior.

Constraints on the deviations expected from single observations of binary black holes by second-generation ground-based detectors are expected to be modest, due to the large statistical errors from the relatively small SNRs expected for most systems (see, e.g., figure 2). However, the constraints can be improved by combining the likelihoods of $\{\epsilon,\sigma\}$ from multiple events.555We expect the combined errors to go down roughly as $N^{-1/2}$, where $N$ is the number of events while the individual errors to go down as $\rho^{-1}$, where $\rho$ is the SNR of the event. Thus, combined constraints from $N$ events with roughly same SNR should be comparable to the constraints from a single event with an SNR of $\rho\,N^{1/2}$. Figure 3 shows the combined posteriors $P(\epsilon,\sigma)$ as a function of the number of simulated events. The constraints on the deviation parameters $\{\epsilon,\sigma\}$ become narrower when multiple events are combined. These constraints are within the reach of Advanced LIGO–Virgo observations in the next few years [21, 5].

From the posterior distribution $P(\epsilon,\sigma)$, one can obtain the credible level corresponding to the GR value of the deviation parameters, i.e., $\{\epsilon,\sigma\}=\{0,0\}$. The credible level for a certain set of parameters $\{\epsilon_{0},\sigma_{0}\}$ is the total probability enclosed within an isoprobability contour passing through $\{\epsilon_{0},\sigma_{0}\}$:

[TABLE]

See, e.g., equation (42) in [26] for a general definition. The area enclosed within the isoprobability contour is the credible area (or credible interval in $1$ dimension). A lower credible level for the GR parameters indicates a greater agreement with GR. However, this GR credible level is not expected to be zero even when the data is correctly described by GR, because the peak of posterior on the deviation parameters $\{\epsilon,\sigma\}$ is randomly shifted away from $\{0,0\}$ in the presence of noise.

If the parameter values $\{\epsilon,\sigma\}$ were sampled from the prior distribution used to compute their posterior, we would expect the $p$ credible interval of the posterior to include the true value approximately $p$ of the time, for the case of an appropriately chosen prior [40, 41]. However, for the GR injections that we performed, we do not expect the true value of $\{\epsilon=0,\sigma=0\}$ to be found at credible level $p$ for a fraction $p$ of events. This is because of our prior distribution for $\{\epsilon,\sigma\}$ does not match the intrinsic distribution (which would be a delta function): we are allowing for the possibility of mismatches between inspiral and merger–ringdown even though this is not the case for the injected GR signals. To demonstrate this for our set of GR injections, in figure 4 we show the fraction of events for which the (true) GR point is found within the $p$ credible region. From figure 4, we see the true GR value is often near the peak of the distribution, meaning that it is found within the $p$ credible region for more than a fraction $p$ of events.666For the marginalised $\epsilon$ distribution, we see that the true GR value is approximately found within the $p$ credible interval a fraction $p$ of the time. In this case, the scatter due to noise is sufficient to balance the mismatch in the prior used. In order to make a frequentist statement about the significance of the measurement of non-zero $\epsilon$ or $\sigma$ (how common such a deviation is if GR were correct), a similar study would need to be done for the region of mass and spin space of interest in order to calibrate expectations.

4 Simulations of modified GR signals in Gaussian noise

4.1 Generation of modified GR waveforms

We also test the sensitivity of our pipeline towards certain kinds of deviations from GR, using signals that differ slightly from the predictions of GR. Specifically, we generate kludge modified GR waveforms that are qualitatively similar to binary black hole waveforms in GR, but differ in their energy and angular momentum loss as a function of frequency. These waveforms are similar enough to GR binary black hole waveforms that they would likely be detected with a standard detection pipeline, as we discuss below. However, we demonstrate that the consistency test described in this paper is able to identify such deviations from GR by combining the posteriors from multiple observations.

We start from the nonspinning effective one-body (EOB) waveform model given in [42], which is available as Matlab code at IHES [43], and modify the GW flux to yield our kludge modified GR waveform. Specifically, we multiply all the modes of the waveform that first enter the flux at second PN order [i.e. $O(v^{4}/c^{4})$, where $v$ is the binary’s orbital velocity], viz the $(\ell,m)=(3,\pm 2)$, $(4,\pm 4)$, and $(4,\pm 2)$ modes, by a constant factor $\sqrt{\alpha_{2}}$, so that their contributions to the flux are multiplied by $\alpha_{2}$. This also modifies the flux at $3$PN and all higher PN orders. We start with the modes that first enter at $2$PN rather than the ones that first enter at $1$PN since the latter modes vanish for equal-mass nonspinning binaries, and we still want a nonvanishing modification in this limit.

We retain the termination condition for the inspiral–plunge phase of the waveform used the original code. This uses the maximum of the orbital frequency (calculated from the EOB Hamiltonian) to mark the termination of the calculation of the inspiral–plunge waveform (by evolving the EOB Hamilton equations) and the start of the matching to the QNMs used to model the merger and ringdown. This termination condition coincides with crossing the light ring in the extreme mass-ratio limit. We also keep the next-to-quasicircular parameters set to the values determined from NR simulations that are already given in the code, for simplicity. Additionally, we keep all the other parameters set to their default values, except for a parameter that keeps one from probing unphysical regions of the EOB potential while the EOB equations are solved numerically. This was necessary to change in order to compute waveforms for higher mass ratios, even without other modifications to the code.

With these choices, we find that the waveforms look similar to a GR waveform up to the values of $\alpha_{2}=20$ we consider here. This is illustrated in figure 5 for an equal-mass binary; the unequal-mass waveforms we consider in the next subsection also look similar to GR waveforms. Additionally, we have checked that this modification to the waveform does not induce any excess eccentricity: The eccentricity of the modified waveforms is small, and the same size as the eccentricities of the unmodified waveforms ($\lesssim 10^{-5}$ for the initial frequencies we are considering).111We estimate the eccentricity of the waveforms as in section IV C of [44] [Eq. (4.13) and following], though with two additional higher-order PN terms added to the fitting function, since the eccentricity is small. The results did not change upon adding a third additional PN term. Since this estimate only uses the leading-order GW effects, it does not need to be modified to account for the modification to the flux, which only starts at $2$PN.

Since the final mass and spin in the original EOB waveform are set by a fit to NR results, we replace this determination by demanding self-consistency of the radiated energy and angular momentum. That is, we choose the final mass and spin by minimizing the difference between the final mass and spin used to obtain the merger–ringdown part of the model and the final mass and spin obtained from energy and angular momentum balance using the initial energy and angular momentum and the radiated energy and angular momentum (calculated from the waveform). Here we are able to obtain fractional disagreements of less than $\sim 10^{-4}$ between the final mass and spin calculated from the radiated energy and angular momentum and those used to compute the merger-ringdown part of the signal. These fractional disagreements are considerably better than the agreement of these quantities for the original IHES EOB waveform for some mass ratios, and much smaller than the differences in radiated energy and angular momentum between the different waveforms we consider. We compute the energy and angular momentum loss using multipoles up to $\ell=7$, since we are only able to fit for the ringdown up to $\ell=7$, the highest $\ell$ for which the tabulated QNM results used in the code are available [14, 45]. We take the standard GR expressions for the radiated energy and angular momentum to still be valid for this modified gravity waveform; these expressions indeed remain valid for a significant number of modified theories [46]. We find that the final mass and spin falls outside the region allowed in GR in many cases, though there is a region of overlap for $\alpha_{2}=20$ and mass ratio $q:=m_{1}/m_{2}\gtrsim 4$; see figure 6. We also plot the fractional difference in the final mass and spin between GR and the kludge waveforms with $\alpha_{2}=20$ in this figure, and see that the fractional difference in the final spin is significantly larger than the fractional difference in the final mass.

In addition to enhanced energy and angular momentum fluxes, these kludge modified GR waveforms have an enhanced linear momentum flux. This leads to larger GW recoils than one obtains for nonspinning systems in GR, where the maximum is only $\sim 180~{}\mathrm{km\,s^{-1}}$ (see, e.g., [47]). For the $\alpha_{2}=20$ deviation we consider here, the maximum recoil is $\sim 440~{}\mathrm{km\,s^{-1}}$.222The original IHES EOB waveform model (without deviations) already has a maximum recoil of $\sim 350~{}\mathrm{km\,s^{-1}}$. We do not include any potential effects (e.g., the Doppler shifts discussed in [48]) from these enhanced recoils, which are still small (at most $\sim 10^{-3}$ times the speed of light).333Enhanced recoils were not an issue for the larger deviation considered in Paper I because the recoil vanishes, by symmetry, for the equal-mass nonspinning system considered there.

We have not changed the QNM spectrum of the final black hole in these kludge modified gravity waveforms for simplicity. However, one expects the QNM spectrum to be modified in any alternative theory that predicts a different GW energy loss for a binary black hole, even those for which the unperturbed Kerr metric is a solution of the field equations, since the behavior of perturbations of Kerr will generically be different in an alternative theory (see, e.g., [49]). Extending this analysis to the case where one modifies the QNM spectrum would be a useful additional check of the consistency test presented here, though the requisite catalogs of QNMs are not yet available, even for the case of a Kerr–Newman black hole, which could be used as a stand-in for a modified gravity black hole: See, however, Dias et al. [50] for some recent results that could be used to obtain such a QNM catalogue.

We only inject the $(2,\pm 2)$ modes of the waveform here, since these are the only modes included in the waveform model we use for the test, as well as for the injections in section 3. The modifications to the higher modes appear in the $(2,\pm 2)$ modes by their effect on the waveform’s phasing, since the higher modes are used to compute the modified energy flux. We have run the test on two GR injections with higher modes in section 5.4, and found that the inclusion of higher modes does not bias the test in those cases.

Finally, almost all of our modified GR signals will be confidently detected by the standard matched-filter based searches for binary black holes, which include the chi-squared discriminatory test [51]. To demonstrate this, we have run the PyCBC [52] matched filter based detection pipeline on these injections (using a single GR template from the SEOBNRv2_ROM_DoubleSpin family, corresponding to the parameters of the simulated modGR signal). We see that the chi-squared weighted SNR recovered by the detection pipeline is close to the optimal SNR of the signal (see table 1 and figure 7 for some examples). The SNR maximized over the full template bank is expected to be even larger than this, practically ensuring confident detections.

4.2 Results from the simulations of modified GR signals

Here we demonstrate that the IMR consistency test is able to identify (at least certain types of) deviations from GR by performing our analysis on a population of simulated signals using the modified GR waveforms described above. In Paper I, we have demonstrated the ability of this test to identify a relatively large modification to the binary’s energy flux ($\alpha_{2}=400$). Such a deviation from GR was easily detectable with high confidence from a single observation of moderate SNR. Here we consider the ability of the test to discern much smaller deviations from GR by combining results from multiple observations. Specifically, we consider the same population of binary parameters considered to simulate the GR case in section 3 except that the waveforms are generated with a modification of the GR energy flux as described in section 4.1. Also, for simplicity, we consider binary black holes with zero spins, since the EOB waveform family that we employ to produce modified GR waveforms is a nonspinning model [42, 43]. However, we still perform parameter estimation using the same SEOBNRv2_ROM_DoubleSpin aligned-spin model employed in section 3.

Modified GR waveforms for binaries with different mass ratios (as well as distances, sky locations, and other extrinsic parameters) were constructed using the prescription presented in section 4.1, with deviation parameter $\alpha_{2}=20$. An example of a modified GR waveform (along with the corresponding GR waveform from the same binary) is shown in figure 5. Although such a modified GR waveforms appears similar to a GR waveform, we show below that the IMR consistency test is able to identify deviations from GR by combining multiple events. Figure 8 shows the joint posteriors on the parameters that describe deviations from GR. It can be seen that, although individual posteriors are unable to identify deviations from GR, the combined posteriors show a clear departure from the GR predictions, even for the relatively small deviation that we consider.444Taking $\alpha_{2}=20$ corresponds to multiplying the $2$PN term in the frequency domain phase expression by a factor of $\sim 0.4$ for equal masses, with the factor decreasing monotonically to $\sim-4.5$ for a mass ratio of $7$ (the largest mass ratio in the detected population). We also see that the combined posterior on $\sigma:=\Delta a_{f}/\bar{a}_{f}$ differs significantly more from [math] than does the combined posterior on $\epsilon:=\Delta M_{f}/\bar{M}_{f}$. This might be expected, due to the larger fractional differences from GR in the final spin seen in figure 6.

4.3 Caveats

The way that we combine the posteriors from multiple events implicitly assumes that $\{\epsilon,\sigma\}$ take the same values in multiple events. This is certainly true if the true theory is GR, when they are always zero, but is not necessarily true for all modified theories of gravity. If the $\{\epsilon,\sigma\}$ values induced by the true theory are strongly dependent on the intrinsic parameters of the black holes (such as their masses and spins), then the combined posteriors from multiple events may not converge to a non-zero value. Hence, this test (like any other test of GR) is not sensitive to all kinds of possible deviations from GR. However, if the deviation is only weakly dependent on the intrinsic parameters, it is a reasonable assumption that the values of $\{\epsilon,\sigma\}$ observed by LIGO-Virgo will be clustered in a small range, due to the rather narrow range of masses ($\sim 10$–$80M_{\odot}$) of black hole binaries employed in this test (and hence the reasonably narrow range of length scales probed by these black holes).555One might also worry about the spin-dependence of possible deviations. However, if the size of the deviation is controlled by the maximum curvature outside the horizons of the black holes, the spin only changes the length scale given by the curvature at the horizon of a Kerr black hole by a factor of $2\sqrt{2}\simeq 2.8$; see, e.g., Eq. (5.47) of [54] for the expression for the Kretschmann scalar in Kerr, whose inverse fourth root gives the length scale we are considering. This neglects possible propagation effects, though even in that case LIGO is not expected to detect binaries with orders of magnitude different distances. This is the case for simulations presented in section 4. Here the modification from GR has a dependence on the intrinsic parameters (see, e.g., figure 6); still, $\{\epsilon,\sigma\}$ estimates from individual posteriors tend to cluster around a reasonably small range of values, and the combined posteriors converge to a non-zero value in figure 8. On the other hand, if the deviation is a strong function of the intrinsic parameters, then combining the posteriors of the entire population may not be the most sensitive method of detecting the deviation. One possibility to address this would be to combine the posteriors of subsets of detections with similar inferred parameters. Such combined posteriors from different regions of the parameter space (e.g., low-mass or high-mass) will provide us some indication of the nature of the deviation as a function of some scales (e.g., curvature). We leave such studies as future work.

Note that, in Figure 8, the fractional difference in the final spin has the opposite sign to what might naïvely be expected, given that the modGR waveforms correspond to smaller final spins than their GR counterparts. Since most of the angular momentum loss comes during the merger and ringdown, it seems reasonable that the inspiral values for the final mass and spin would be close to their GR values for the system’s mass ratio, while those inferred from the ringdown would correspond to the true values used to construct the waveform. However, such an interpretation is possible only when the parameters are estimated from only the early inspiral and late ringdown parts of the waveform. Such a test would require both the early inspiral and late ringdown to be observed with high SNR. While such tests might be possible in the future, currently we are limited by the moderate SNRs of the observed signals, and hence are forced to include the late inspiral and merger parts in our analysis. For the case of the modGR waveforms that we use, the estimated parameters from the inspiral part are biased (typically to more massive, asymmetric binaries) which correspond to a smaller value of the final spin, as compared to the same estimated from the post-inspiral part, producing the observed sign in the fractional difference in the final spin.

5 Robustness of the consistency test

We also test the robustness of the IMR consistency test to the choice of various parameters used in the analysis, for example, the cutoff frequency used to demarcate the division between the low-frequency (inspiral) and high-frequency (merger–ringdown) parts of the waveform, the choice of the particular approximant waveform, and the fitting formula for the mass and spin of the final black hole, etc. For this we use a simulated GR signal with waveform approximant SEOBNRv2_ROM_DoubleSpin (except in section 5.3 where we use approximant IMRPhenomPv2 [55, 56, 57] and in section 5.4 where we use NR waveforms). We also study the effect of neglecting higher modes and spin precession on the IMR consistency test. These are performed by injecting NR waveforms from the publicly available Simulating eXtreme Spacetimes (SXS) catalog [58]. All these studies are performed assuming a binary black hole signal with parameters close to that of the first LIGO event GW150914 — masses $m_{1}=36{M}_{\odot}$, $m_{2}=29{M}_{\odot}$ and aligned spins, $a_{1,z}=-0.32$, $a_{2,z}=0.58$, with optimum sky position and orientation, producing an optimal SNR of $25$ in the Advanced LIGO Hanford–Livingston and Advanced Virgo network. These studies demonstrate the robustness of the IMR consistency test only for the case of signals with parameters similar to that of GW150914, that is, with moderate mass ratios. If we observe signals with large mass ratios and large, misaligned spins in the future, the robustness of the IMR consistency test will need to be reinvestigated. Also, these robustness studies have been restricted to the case of single events: when a large number of events are combined to produce precise constraints on the deviations from GR, we will have to worry about even small systematic errors affecting our analysis. Such studies have to be performed in the near future, in anticipation of the large number of binary black hole signals that Advanced LIGO is expected to observe.

5.1 Cutoff frequency between the inspiral and merger–ringdown

There is no well defined transition frequency (or time) between the inspiral and post-inspiral (merger–ringdown) parts of the waveform [59]. We have chosen the ISCO frequency of the final Kerr black hole with mass and spin inferred from the full IMR waveform at this cutoff frequency. Increasing this cutoff frequency will increase (decrease) the SNR of the inspiral (post-inspiral) part and hence will improve (worsen) the parameter estimation from the inspiral (post-inspiral) analysis, thus affecting our statistical errors. There thus is likely an optimal choice for the cutoff frequency as far as the test’s sensitivity to a given deviation from GR is concerned. However, if our analysis is free from major systematic errors, the two independent estimates always have to be consistent with each other when it is applied to a binary black hole coalescence in GR. Here we illustrate that our analysis produces consistent results for reasonable choices of the cutoff frequency.

The ISCO frequency of a Kerr black hole with mass $M_{f}=62M_{\odot}$ and dimensionless spin $a_{f}=0.68$ (remnant of the merger of two black holes with masses $m_{1}=36{M}_{\odot}$, $m_{2}=29{M}_{\odot}$ and aligned spins, $a_{1,z}=-0.32$, $a_{2,z}=0.58$) is $146.4~{}\mathrm{Hz}$. We repeat the analysis with 4 different choices of the cutoff frequency in the interval $50$–$150~{}\mathrm{Hz}$. The results are shown in figure 9 (top panels). It can be seen that, while the spread of the posteriors (that is, the statistical errors) depends on the choice of the cutoff frequency, the inspiral and post-inspiral estimates are always consistent with each other. Also, the posteriors on the parameters describing deviations from GR are always consistent with zero, indicating the robustness of the test against the specific choice of the cutoff frequency.

5.2 Waveform approximant

The Bayesian inference described in section 2 is performed by employing semi-analytical gravitational waveforms (computed by a combination of analytical and numerical relativity) as the GR model. The inherent assumption is that these are faithful representations of the actual signals produced by nature. However, depending on the particular methods used to construct these waveforms, there can be minor differences between different gravitational waveform families (approximants). Here, we demonstrate the robustness of the IMR consistency test employing two IMR waveform families, namely, the non-precessing spin reduced-order EOB model SEOBNRv2_ROM_DoubleSpin [34, 35] and the precessing single effective-spin phenomenological model IMRPhenomPv2 [55, 56, 57]. The systematic errors of these models in the vicinity of GW150914 and their effects on parameter estimation were studied in [60].

Figure 9 (second row) shows the results from simulated GR signals from a binary black hole with parameters described in the previous section. Apart from performing the analysis where the same waveform family is employed in estimating the parameters from the inspiral and post-inspiral parts, the figure also demonstrates the robustness of the consistency test when we employ one waveform family for estimating the parameters from the inspiral part and another to estimate the same from the post-inspiral part. In each case, the waveform family used to create a portion of the waveform, inspiral or merger-ringdown, is the same used to recover the parameters from that portion.

5.3 Fit formulas for the mass and spin of the final black hole

The consistency test we have developed relies on NR fitting formulas to map the binary’s initial masses and spins into the mass and spin of the remnant black hole. There are a number of such fits available in the literature, of varying degrees of accuracy and generality. Here we consider a variety of recent fits for the final mass and spin.666While there are fits for the final mass and spin used in the two waveform models we employ, we do not just use those here, since they are not the most accurate and up-to-date ones: see Jiménez-Forteza et al. [61], which compares the final mass fit used in SEOBNRv2 and the final mass and spin fits used in IMRPhenomPv2 with more accurate fits. Note, however, that the additional calibrated terms in the merger-ringdown parts of these models may compensate for some of the slight inaccuracies of the final mass and spin fits used. We consider one set that is only applicable to nonspinning binaries (so we do not use the estimated spins in applying this) given in equations (29) of Pan et al. [62], and two sets applicable to binaries with aligned spins. The first set is from Healy, Lousto and Zlochower (HLZ) [39] and involves reasonably complicated implicit expressions given in Eq. (14) and Eq. (16) in that paper, with coefficients given in Table XI (we use the fourth-order fits). The second, simpler, though somewhat less accurate set is from Eq. (3.6) and Eq. (3.8) in Husa et al. [63] and only uses a single effective spin. The latter two fits were used in [3, 4] to infer the final mass and spin of GW150914 and the HLZ fits were used in the implementation of the current test in Paper I and [9]. For these fits, we use the projection of the spins along the orbital angular momentum when using a precessing waveform model, e.g., IMRPhenomPv2. The final mass and spin used to determine the ISCO frequency is not changed throughout the analysis, even when we are changing the fitting formulas. The calculation of the ISCO frequency uses the HLZ fits.

It is also possible to augment aligned-spin fits for the final spin to include the contribution from the in-plane spins, as was done for the HLZ fit in [64] and significantly increases its accuracy for precessing systems. As that reference shows for the HLZ fit, the basic aligned-spin final mass fit is already accurate for precessing systems when evaluated using the components of the spins along the orbital angular momentum. This augmentation was applied to the HLZ fit to infer the final spin of the binary black hole events from Advanced LIGO’s first observing run in [5, 6, 65]. We thus also consider this extension of the final spin expression for the HLZ and Husa et al. fits (HLZ aug and Husa et al. aug), as well as a new final spin fit for precessing systems from Hofmann, Barausse and Rezzolla (HBR) [66] (we use the one with the smallest number of coefficients, which the authors recommend), which we pair with the HLZ final mass fit. We also consider the augmented versions of both new aligned-spin fits that extend those in Husa et al. to include the spin difference (UIB aug) [61], as well as a small update to the HLZ fits from Healy and Lousto (HL aug) [47]. An average of these fits (HBR, UIB aug, and HL aug for the final spin, and just the second two for the final mass) was used to infer the final mass and spin of the binary coalescence that produced GW170104, as well as in the application of the IMR consistency test to that system [7].

In figure 9 (third row), we show that for the simple test case we consider, the two aligned-spin fits give almost identical results, as do the extensions using the in-plane spins, while the nonspinning fits still give a posterior that is consistent with zero. Here we use the precessing waveform model IMRPhenomPv2 [55, 56, 57] for both the injection and recovery. We use it for the recovery in order to have nonzero in-plane spins, and use it (in its aligned-spin limit) for the injection so we are only testing different fits, not different waveform models here. Since we apply the same fitting formulas to infer the final mass and spin from both the inspiral and post-inspiral portions of the waveform, systematic errors in the fitting formulas do not translate into systematic errors in the test. Nevertheless, having more accurate fitting formulas might help the test to identify modified GR signals.

5.4 Effects of spin precession and higher modes

Although GW emission in the leading order is quadrupolar, higher order (nonquadrupole) modes can contain nonnegligible power if the system has large asymmetries (e.g., highly unequal masses or large, misaligned spins). Due the unavailability of waveform models that are accurate and computationally fast to generate, currently the parameter estimation is mostly performed making use of waveform templates that only model the quadrupole ($\ell=2,m=\pm 2$) modes. This can potentially introduce systematic errors in the estimated parameters, which could be falsely taken for a deviation from GR. However, the systematic errors due to neglecting nonquadrupole modes are expected to be negligible for the case of binary black holes with moderate mass ratios and aligned/antialigned spins [60, 67]. In this paper’s population study we considered binaries with spins aligned/antialigned with the orbital angular momentum so that there is no precession. However, fast-to-evaluate waveform templates that take into account the dominant precession effects are available [55, 56, 57], which have been used for performing the IMR consistency test on the LIGO events GW150914 [9] and GW170104 [7].

We investigate the robustness of the IMR consistency test against the effects of precession and higher modes, making use of injected NR waveforms from the SXS waveform catalog [68, 58, 69]. We select two unequal mass ($m_{1}/m_{2}=1.228$) waveforms used in [60]: SXS:BBH:0307, which has aligned spins ($a_{1,z}=0.32$, $a_{2,z}=-0.5798$777The aligned spin components for the NR injections used in section 5.4 have signs opposite to the ones used in the rest of the section; this is not a typographical error.) binary black hole system, and SXS:BBH:0308, which is a precessing system with spin magnitudes $a_{1}=0.3406$, $a_{2}=0.6696$ and aligned components $a_{1,z}=0.3224$, $a_{2,z}=-0.5761$. We inject both systems with total masses of $65M_{\odot}$ (so their individual masses are $36$ and $29{M}_{\odot}$, the same as for the previous injections), and with inclination angles of [math]. For parameter estimation, we use SEOBNRv2_ROM_DoubleSpin as our template. The results are summarized in figure 9 (bottom panel): the test is robust against the presence of higher modes and spin precession in the comparable-mass regime.

5.5 Effect of splitting the signal in the frequency domain

In the IMR consistency test, we check the consistency of the posteriors of the final mass and spin estimated from the early part of the signal (inspiral) with the same estimated from the late part of the signal (merger–ringdown). The split between the inspiral and merger–ringdown parts is done in the Fourier domain. This choice is purely made for convenience, since, in the Fourier domain, this amounts to simply setting the lower and upper limits of the likelihood integral in equations (2) and (3).888This is different from splitting the signals in time domain and then taking the Fourier transform, which will introduce additional artifacts in the Fourier transform due to edge effects or windowing. There is a possibility that a fraction of the power from the inspiral can get deposited at high Fourier frequencies, or the power from the merger–ringdown parts to get deposited at low Fourier frequencies. This effect will be particularly pronounced for the ringdown signal, whose power is spread over a large range of Fourier frequencies. One may worry that this spectral leakage would cause the posteriors derived from the two regions of Fourier frequency to be automatically consistent (although, we have already demonstrated that when certain kinds of deviations from GR are present in the signal, the IMR consistency test is able to identify them; see section 4.2).

Below we show that the effect of spectral leakage is small for the case of cutoff frequencies that we choose. In order to demonstrate this, we do the following: We compute the Fourier transform of the NR waveform SXS:BBH:0253 (an aligned-spin system with $m_{1}/m_{2}=2$ and $a_{1,z}=a_{2,z}=0.5$) [68, 58] making use of the stationary phase approximation (SPA) [e.g., equation (3.7) of [70]]. The approximation here is that the power at a Fourier frequency $f$ comes entirely from the time $t_{f}$ (the saddle point) when the instantaneous frequency $F(t):=(\mathrm{d}\varphi(t)/\mathrm{d}t)/2\pi$ is equal to $f$. Here $\varphi(t)$ is the phase of the waveform. In this approximation, the magnitude of the Fourier transform of the NR waveform can be computed as

[TABLE]

where $A(t)$ is the time-domain amplitude of the NR waveform, and a dot denotes a time derivative.

If the Fourier transform computed using the SPA agrees well with numerically computed Fast Fourier Transform (FFT), this gives a strong indication that the spectral leakage is negligible. Figure 10 compares the Fourier transform computed using the SPA with the FFT. The corresponding ISCO frequency $f_{\mathrm{ISCO}}$ is also shown, which is used to demarcate the inspiral and merger–ringdown parts. The plot suggests that the power at Fourier frequencies less than $f_{\mathrm{ISCO}}$ can be fully explained to be coming from the early times of the waveform (i.e., with instantaneous frequency $F(t)<f_{\mathrm{ISCO}}$). The excellent agreement between the SPA and FFT at frequencies less than $f_{\mathrm{ISCO}}$ suggests that there is no appreciable spectral leakage between the two bands — almost all the power in the $f<f_{\mathrm{ISCO}}$ ($f>f_{\mathrm{ISCO}}$) band should come from the early (late) times. However, if we choose a transition frequency much higher than $f_{\mathrm{ISCO}}$ to demarcate the inspiral and merger–ringdown (say, the dominant QNM frequency, which is also shown in the figure) we expect significant spectral leakage between the two bands. Thus, if we need to perform this same analysis purely on the ringdown part (without including the merger), the analysis must be performed in the time domain. There is ongoing work in this direction.

6 Conclusions and future work

In this paper we have provided a detailed description of a test of GR based on the consistency between the inspiral, merger and ringdown parts of observed GW signals from coalescing binary black holes. First presented in [10], based on the idea proposed by Hughes and Menou [24], this was one of the handful of tests used to establish the consistency of the GW150914 event with a binary black hole system described by GR [9]. Here, we demonstrate how the likelihoods from multiple binary black hole events can be combined to produce tighter constraints on parameters describing a deviation from GR. In order to illustrate how this test might be able to detect certain types of deviations from GR, we generate internally consistent kludge waveforms which disagree with the GR prediction of the radiated energy and angular momentum, by increasing the flux radiated into GWs. We then simulate a population of binaries modeled by these modified GR waveforms. By combining the results from multiple events, we demonstrate how this test might be able to detect small deviations from GR. For the case of single events with parameters similar to LIGO’s first binary black hole observation, we also demonstrate the robustness of the consistency test against specific choices of various analysis parameters, such as the choice of the transition frequency used to demarcate the inspiral and merger–ringdown parts, the specific waveform approximant, fitting formulas for the mass and spin of the remnant black hole, etc.

While we have investigated the reliability of our test to ensure its robustness, there remain further extensions of the test and these analyses to pursue. In particular, we have assumed that black holes in binaries have nonprecessing spins in the populations we simulate in this paper. We have also neglected the effect of nonquadrupole modes in the parameter estimation (although the fitting formulas for the remnant mass and spin include their contributions to the energy and angular momentum loss). While the effects of precession and nonquadrupole modes are expected to be negligible for the parameter estimation of comparable-mass binaries (see, e.g., [60, 67]) their effects might be appreciable for binaries with large mass ratios and high spins. Waveform templates describing the spin precession through an effective spin parameter, notably the IMRPhenomPv2 template [55, 56, 57], are already available and have been used in parameter estimation. IMRPhenomPv2 has been employed in the application of this test to GW150914 in [9], to GW170104 in [7], and in some of the robustness tests in section 5. More recently, waveform templates describing the full double-spin precession effects have also become available [71, 72]. These templates can be employed in parameter estimation when precession effects are expected to be apparent. Fitting formulas for the remnant mass and spin valid for the case of precessing initial spins are already available; some of them have been used in this paper itself. Waveform templates for binaries of spinning black holes that include the effect of nonquadrupole modes have started to become available [73], and they may help to avoid possible systematic errors due to neglecting nonquadrupole modes.

The robustness studies performed in this paper are restricted to the case of single events with modest SNR. When a large number ($\sim 100$) of events are combined to produce precise constraints on deviations from GR, even small systematic errors could become significant. Apart from the aspects discussed above, possible sources of errors include the finite accuracy of the GR waveform templates and the calibration of the GW detectors. We already have preliminary indications that combined posteriors from a large number of events could be dominated by such errors. Therefore, a careful characterization of various systematic errors is required before we accumulate a large number of detections to analyze. We leave this as future work.

Acknowledgments.

We are grateful to M. Vallisneri and the anonymous referees for their insightful comments on the manuscript. We thank J. Veitch, A. Nagar and P. Schmidt for assistance with the LALInference, IHES EOB and NR injection codes, respectively. We also thank the SXS Collaboration for creating a public archive of NR waveforms. We have benefited from useful discussions with K. G. Arun, A. Buonanno, N. Christensen, M. K. Haris, B. R. Iyer, S. Kumar, A. K. Mehta, C. Messenger, A. Mukherjee, B. S. Sathyaprakash, M. Vallisneri, C. Van Den Broeck, S. Vitale and several members of the LIGO Scientific and Virgo Collaboration’s working group for strong gravity. Ar. G., N. K. J.-M., and P. A. acknowledge support from the AIRBUS Group Corporate Foundation through a chair in “Mathematics of Complex Systems” at ICTS. P. A.’s research was, in addition, supported by a Ramanujan Fellowship from the Science and Engineering Research Board (SERB), India, the SERB FastTrack fellowship SR/FTP/PS-191/2012, and by the Max Planck Society and the Department of Science and Technology, India through a Max Planck Partner Group at ICTS. W. D. P. was partly supported by a Leverhulme Trust research project grant. C. P. L. B. was supported by the Science and Technology Facilities Council. Computations were performed at the ICTS clusters Mowgli, Dogmatix, and Alice. This is LIGO document LIGO-P1700006.

Appendix A Calculation of the posterior of the parameters describing deviations from GR

Here we describe the calculation of the posterior $P(\epsilon,\sigma\,|\,d)$ of the parameters $(\epsilon,\sigma)$ describing the fractional difference between the two independent estimates of the mass and spin of the remnant. We first make a change of variables in the joint posterior $P_{\textsc{i}}(M^{\textsc{i}}_{f},a^{\textsc{i}}_{f}\,|\,d)P_{\textsc{mr}}(M^{\textsc{mr}}_{f},a^{\textsc{mr}}_{f}\,|\,d)$ from $\{M^{\textsc{i}}_{f},a^{\textsc{i}}_{f},M^{\textsc{mr}}_{f},a^{\textsc{mr}}_{f}\}$ to $\{\epsilon,\sigma,\bar{M}_{f},\bar{a}_{f}\}$ and then marginalize over $\bar{M}_{f}$ and $\bar{a}_{f}$. Starting from the definitions given in equations (5) and (7), we find that

[TABLE]

so that the Jacobian of the transformation from $\{\epsilon,\sigma,\bar{M}_{f},\bar{a}_{f}\}$ to $\{M^{\textsc{i}}_{f},a^{\textsc{i}}_{f},M^{\textsc{mr}}_{f},a^{\textsc{mr}}_{f}\}$ is $\bar{M}_{f}\bar{a}_{f}$. Thus, the final expression is

[TABLE]