{'fraction_non_alphanumeric': 0.03258785942492013, 'fraction_numerical': 0.012328509678631836, 'mean_word_length': 4.633178064789329, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Testing general relativity using gravitational wave signals from the inspiral, merger and ringdown of binary black holes

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.


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 ∼ 36M and 29M ) at a distance of ∼ 400 Mpc from the Earth [4,5]. This observation was followed by a second one on December 26, 2015 [6]. The second signal, termed GW151226, was produced by the coalescence of two black holes less massive (∼ 14M and 8M ) than those in GW150914 [5]. 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 ( 25M ) [7]. 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 [8,5] -a regime with binary orbital velocities ∼ 0.5 c, which is inaccessible by other astronomical observations and laboratory tests.
A handful of tests were 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 on GW150914 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 [9]. 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 [10,11,12,13]. 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 [14,15]. 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 [16]. In addition, the residual after subtracting the best fit GR signal was found to be consistent with detector noise. Within the statistical uncertainties, these investigations provided no evidence for deviations from GR [8].
The Advanced LIGO detectors, which are currently operating at less than their design sensitivities, are expected to approach design sensitivity over the next few years [17]. Several observing runs will be conducted over this period. Based on the rate of GW events in the first observing run, 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 [18,5]. The Advanced Virgo [19] detector in Italy and KAGRA [20] in Japan are expected to join the international network of advanced GW detectors in the next few months to years, and the LIGO-India [21] detector is planned to join this network by 2024. 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 sister projects, we present here a detailed follow up of the analysis method of the inspiral-merger-ringdown (IMR) consistency test for binary black holes presented in our previous paper [9], which we will henceforth refer to as Paper I.
The IMR consistency test was developed based on the original idea by Hughes and Menou [22] proposed in the context of the space-based LISA observatory. † 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 (lowfrequency) 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 (highfrequency) 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 [22] 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). 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 lowfrequency (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.

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 [24]). 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 λ := {λ i } of a model hypothesis H, given data d and any other information I, is: where P(λ|H, I) is the prior probability of the parameter set λ given H and I, while L(d|λ, H, I) is called the likelihood function, which is the probability of observing the data d, given λ, H and I. What appears in the denominator is a normalization constant P(d|H, I) := P(λ|H, I) L(d|λ, H, I) dλ, called the marginal likelihood, or the evidence of the hypothesis H. If our hypothesis H is that the data contains a GW signal described by a GR waveform model h(λ) and stationary Gaussian noise described by the power spectral density S n ( f ), then, as described in appendix A of [25], the likelihood function can be defined as: where ·|· is the following noise-weighted inner product Above,B( f ) denotes the Fourier transform of B(t) and a * denotes complex conjugation. The limits of integration f low and f 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 λ, the posterior distribution P(λ|d, H, I) in equation (1) is computed by stochastically sampling the parameter space using techniques such as Markov-Chain Monte Carlo [26] or nested sampling [27]. For this paper, we use the LALInferenceNest [28,24] code that provides an implementation of the nested sampling algorithm in the LSC Algorithms Library [29] 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 , a 1 , a 2 | d) † of the binary black hole system (with m 1 ≥ 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, 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 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 isco of the innermost stable circular orbit (ISCO) of the remnant Kerr black hole [31], with mass and spin given by the median values of P 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 cut = f isco in equations (2) and (3). This allows us to compute the posterior distribution P 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 low = f isco in equations (2) and (3)], we can get yet another estimate of the posterior P 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 i (M f , a f | d) and P mr (M f , a f | d) have to be consistent with each other as well as with the estimate P 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 where In the absence of departures from GR, we expect the posterior P( , σ | d) of { , σ} to be consistent with zero (within the expected statistical fluctuations due to noise). ‡ Appendix A describes the calculation of P( , σ | d) from the posteriors P i (M f , a f | d) and P mr (M f , a f | d). An example of the posterior distribution P( , σ | d) from a simulated GR signal is shown in figure 1 (right plot). † Here onwards, we drop 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 [30]. ‡ We 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. Finally, the likelihoods from multiple observations of binary black holes can be combined (multiplicatively) to construct a combined posterior on { , σ} that can better constrain deviations from GR. Here, we use the posterior from previous observations as the prior for the current observation. In practice, the posteriors P j ( , σ | d j ) from N such observations can be combined to produce an overall posterior in the following way: where P( , σ) is the overall prior distribution on and σ and P j ( , σ) is the prior distribution used to compute the posterior P j ( , σ | d j ).

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 [32,33]) into colored Gaussian noise with the design power spectrum of the Advanced LIGO detectors in the high-power, zero-detuning configuration [34]. 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 low = 10 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 [35]  the detector frame become m 1,2 = m s 1,2 (1 + z) [30]. 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 [36]. 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 †. In addition, we only consider binaries with total redshifted mass less than 150M , so that the observed signals contain imprints of the inspiral, merger and ringdown.
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. [37]. Parameter estimation using LALInferenceNest is performed using priors that are uniform in component masses in the interval m 1,2 ∈ [1, 300]M and component (aligned/antialigned) spins in the interval [0, 0.98]. This corresponds to a non-uniform prior in the deviation parameters { , σ} 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 { , σ} from multiple events. Figure 3 shows the combined posteriors P( , σ) as a function of the number of simulated events. The constraints on the deviation parameters { , σ} become narrower when multiple events are combined. These constraints are within the reach of Advanced LIGO-Virgo observations in the next few years [18,5].
From the posterior distribution P( , σ), one can obtain the credible level corresponding to the GR value of the deviation parameters, i.e., { , σ} = {0, 0}. The credible level for a certain set of parameters { 0 , σ 0 } is the total probability enclosed within an isoprobability contour passing through { 0 , σ 0 }: See, e.g., equation (42) in [24] 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 { , σ} is randomly shifted away from {0, 0} in the presence of noise.
If the parameter values { , σ} 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 [38,39]. However, for the GR injections that we performed, we do not expect the true value of { = 0, σ = 0} to be found at credible level p for a fraction p of events. This is because of our prior distribution for { , σ} 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. † In order to make a frequentist statement about the significance of the measurement of non-zero or σ (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.

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 [40], which is available as Matlab code at IHES [41], 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 ( , m) = (3, ±2), (4, ±4), and (4, ±2) modes, by a constant factor √ α 2 , so that their contributions to the flux are multiplied by α 2 . This also modifies the flux at 3PN and all higher PN orders. We start with the modes that first enter at 2PN rather than the ones that first enter at 1PN 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    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 α 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 ( 10 −5 for the initial frequencies we are considering). † 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 ∼ 10 −4 , which is considerably better than the agreement of 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 = 7, since we are only able to fit for the ringdown up to = 7, the highest for which the tabulated QNM results used in the code are available [13,43]. 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 [44]. 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 α 2 = 20 and mass ratio q := m 1 /m 2 4; see figure 6. We also plot the fractional difference in the final mass and spin between GR and the kludge waveforms with α 2 = 20 in this figure, and see that the † We estimate the eccentricity of the waveforms as in section IV C of [42] [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 2PN. 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 ∼ 180 km s −1 (see, e.g., [45]). For the α 2 = 20 deviation we consider here, the maximum recoil is ∼ 440 km s −1 . † We do not include any potential effects (e.g., the Doppler shifts discussed in [46]) from these enhanced recoils, which are still small (at most ∼ 10 −3 times the speed of light). ‡ 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., [47]). 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. [48] for some recent results that could be used to obtain such a QNM catalogue.

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 (α 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 [40,41].
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 α 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 7 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. † We also see that the combined posterior on σ := ∆a f /ā f differs significantly more from 0 than does the combined posterior on := ∆M f /M f . This might be expected, due to the larger fractional differences from GR in the final spin seen in figure 6. ‡ † The original IHES EOB waveform model (without deviations) already has a maximum recoil of ∼ 350 km s −1 .
‡ Enhanced recoil 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. † Taking α 2 = 20 corresponds to multiplying the 2PN term in the frequency domain phase expression by a factor of ∼ 0.4 for equal masses, with the factor decreasing monotonically to ∼ −4.5 for a mass ratio of 7 (the largest mass ratio in the detected population). ‡ This 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 reasonably 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, since the analysis uses the merger-ringdown (and not just the ringdown) in order to have sufficient SNR, we do not expect such a direct correspondence to the values used to construct the waveform. For high-SNR events, the final mass and spin can be inferred directly using just the ringdown, where the only complication is the Doppler shift from any recoil. In this case, the initial masses and spins could also be inferred using only frequencies well below the merger. This would then allow the test to be interpreted as a missing energy experiment, where it could be concluded that additional energy (than predicted by GR) had been lost if the final black hole had a smaller mass than was  Table 1. The table shows the parameters of some of the injected modified GR signals (D 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 [50]. We only inject the (2, ±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, ±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 [49]. To demonstrate this, we have run the PyCBC [50] 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 8 for some examples). The SNR maximized over the full template bank is expected to be even larger than this, practically ensuring confident detections.

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 expected from the inspiral. 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 [51,52,53] 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 [54]. 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 = 36M , m 2 = 29M 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.

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 [55]. 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 and dimensionless spin a f = 0.68 (remnant of the merger of two black holes with masses m 1 = 36M , m 2 = 29M and aligned spins, a 1,z = −0.32, a 2,z = 0.58) is 146.4 Hz. We repeat the analysis with 4 different choices of the cutoff frequency in the interval 50-150 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.

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 [32,33] and the precessing single effective-spin phenomenological model IMRPhenomPv2 [51,52,53]. The systematic errors of these models in the vicinity of GW150914 and their effects on parameter estimation were studied in [56]. 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.

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. † 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. [58], and two sets applicable to binaries with aligned spins. The first set is from Healy, Lousto and Zlochower (HLZ) [37] 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. [59] 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 [8]. 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 [60] 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,61]. 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) [62] (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) [57], as well as a small update to the HLZ fits from Healy and Lousto (HL aug) [45].
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 [51,52,53] 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 distinguish modified GR signals from GR signals.

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., large mass ratios 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 ( = 2, m = ±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 [56,63]. 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 [51,52,53], which have been used for performing the IMR consistency test on the first LIGO event, GW150914 [8].
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 [64,54]. We select two unequal mass (m 1 /m 2 = 1.228) waveforms used in [56]: SXS:BBH:0307, which has aligned spins (a 1,z = 0.32, a 2,z = −0.5798 †) 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 (so their individual masses are 36 and 29M , the same as for the previous injections), and with inclination angles of 0. 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.

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). ‡ 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 then; 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) [64,54] making use of the stationary phase approximation (SPA) [e.g., equation (3.7) of [65]]. 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) := (dϕ(t)/dt)/2π 40   . Robustness studies against the choice of cutoff frequency (top row), approximants (second row), fitting formulas (third row) and effects of higher modes and precession (bottom row). Left panels show the 68% credible region of the posterior distributions P( , σ). The GR value is marked by a + sign. Right panels show the 68% credible region of the posterior distributions P i (M f , a f ), P mr (M f , a f ) and P imr (M f , a f ) of the mass and spin of the final black hole estimated from the inspiral, merger-ringdown parts and the full IMR signal for each simulated case, respectively. In the third panel, the contours for the five fits that include the contributions from in-plane spins to the final spin (HLZ aug [37], Husa et al. aug [59], HBR [62], UIB aug [57] and HL aug [45]; "aug" denotes augmentation following [60]) lie almost on top of each other, thus making the individual contours indistinguishable. The same is true for the contours for the purely aligned-spin fits (HLZ [37] and Husa et al. [59]). The nonspinning Pan et al. [58] fit has a distinct contour. The corresponding ISCO frequency f 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 ISCO can be fully explained to be coming from the early times of the waveform. Also shown is the dominant QNM frequency, for reference. At high frequencies, the agreement between SPA and FFT becomes poor, indicating spectral leakage. The noise in the SPA at lower frequencies is due to the numerical noise in computing the frequency derivative.
is equal to f . Here ϕ(t) is the phase of the waveform. In this approximation, the magnitude of the Fourier transform of the NR waveform can be computed as 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 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 ISCO can be fully explained to be coming from the early times of the waveform (i.e., with instantaneous frequency F(t) < f ISCO ). The excellent agreement between the SPA and FFT at frequencies less than f ISCO suggests that there is no appreciable spectral leakage between the two bands -almost all the power in the f < f ISCO ( f > f ISCO ) band should come from the early (late) times. However, if we choose a transition frequency much higher than f 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.

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 [9], based on the idea proposed by Hughes and Menou [22], 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 [8]. 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., [56,63]) 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 [51,52,53], are already available and have been used in parameter estimation. IMRPhenomPv2 has been employed in the application of this test to GW150914 in [8] 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 [66,67]. 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. However, waveform templates for binaries of spinning black holes that include the effect of nonquadrupole modes are currently under development. We anticipate them to be available in the near future, 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 (∼ 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.