Testing general relativity using golden blackhole binaries
Abstract
The coalescences of stellarmass blackhole binaries through their inspiral, merger, and ringdown are among the most promising sources for groundbased gravitationalwave (GW) detectors. If a GW signal is observed with sufficient signaltonoise ratio, the masses and spins of the black holes can be estimated from just the inspiral part of the signal. Using these estimates of the initial parameters of the binary, the mass and spin of the final black hole can be uniquely predicted making use of generalrelativistic numerical simulations. In addition, the mass and spin of the final black hole can be independently estimated from the merger–ringdown part of the signal. If the binary black hole dynamics is correctly described by general relativity (GR), these independent estimates have to be consistent with each other. We present a Bayesian implementation of such a test of general relativity, which allows us to combine the constraints from multiple observations. Using kludge modified GR waveforms, we demonstrate that this test can detect sufficiently large deviations from GR, and outline the expected constraints from upcoming GW observations using the secondgeneration of groundbased GW detectors.
I Introduction
The coalescence of blackhole binaries, driven by the emission of gravitational radiation, is perhaps the most luminous phenomenon occurring in the Universe after Big Bang. During the final stages of the coalescence, up to of the massenergy of the binary is radiated as gravitational waves (GWs) over the last few orbits of the inspiral and merger (see, e.g., Centrella et al. (2010) for a review). This will allow the secondgeneration groundbased GW observatories Aasi et al. (2015); Acernese et al. (2015); Somiya (2012); Iyer et al. (2011); Abbott et al. (2016) to detect such phenomena up to distances of several gigaparsecs Ajith et al. (2008), making binary black hole coalescences some of the most promising sources of GWs for these observatories. Expected detection rates for secondgeneration detectors vary from a handful to several thousands per year, as predicted by population synthesis models Dominik et al. (2015); Rodriguez et al. (2015). Thirdgeneration detectors Punturo et al. (2010); Miao et al. (2014); Dwyer et al. (2015) are expected to extend the range even further.
GW observations of binary black holes will enable us to test general relativity (GR) in a regime that is currently inaccessible by astronomical observations and laboratory tests. Apart from putting bounds on parameters of specific alternative theories, proposed tests include constraining parametrized deviations from postNewtonian gravity, tests of the nohair theorem by observing multiple quasinormal modes or by constraining deviations from the expected multipolar structure of black holes, etc. (see, e.g., Berti et al. (2015); Yunes and Siemens (2013) for reviews). Here we present a test of GR based on GW observations of “golden” blackhole binaries Hughes and Menou (2005); Nakano et al. (2015) – binaries with total mass –, so that the signals observed by groundbased GW observatories cover the inspiral, merger and ringdown (IMR) phases of the coalescence. During the inspiral, the two black holes spiralin under gravitational radiation reaction, and eventually merge to form a common horizon. In the ringdown stage, the newly formed horizon settles into a Kerr black hole with the emission dominated by a spectrum of quasinormal modes. According to the nohair theorem, the final black hole is fully characterized by its mass and spin angular momentum.
The idea of the proposed test is that, if a GW signal is observed with sufficient signaltonoise ratio (SNR), the masses and spins of the black holes can be estimated just from the inspiral part of the signal. Given the estimates of the initial parameters of the binary, the mass and spin of the final black hole can be uniquely predicted making use of fits to numericalrelativity (NR) simulations. In the same way, the mass and spin of the final black hole can be independently estimated from the mergerringdown portion of the signal.^{1}^{1}1The original test proposed in Hughes and Menou (2005) in the context of LISA makes use of only the inspiral and ringdown parts. In the case of secondgeneration groundbased detectors, the ringdown SNR is unlikely to be large for most events. Luckily, recent advances in NR have allowed us to model the merger accurately. Hence, our implementation makes use of the merger part as well. However, it is possible to restrict our test solely to the inspiral and ringdown parts by appropriate choice of the cutoff frequencies defined in Eq. (2). If the binary black hole dynamics is correctly described by GR, these independent estimates have to be consistent with each other. The consistency of the parameters estimated from the highly relativistic postinspiral regime with those inferred from the weakly relativistic inspiral regime is a nontrivial test of the ability of GR in modeling this complex phenomenon.
Ii Formulation of the test
The set of parameters of the binary, such as the masses () and spin angular momenta () of the black holes, are imprinted on the gravitational waveform. Given data containing an observed GW signal, and assuming the GR model , the posterior distribution of these parameters can be estimated making use of Bayes’ theorem
(1) 
where is the prior distribution of , is a normalization constant (called the evidence) and is the likelihood of observing the data given the signal model and the set of parameters ,
(2) 
Above, is the Fourier transform of the data, is the frequencydomain signal waveform corresponding to the set of parameters , and is the power spectral density of the detector noise, while and are the lower and upper cutoff frequencies used in the calculation. The sampling of the likelihood function over the (typically large dimensional) parameter space often makes use of stochastic sampling methods such as Markovchain MonteCarlo or nested sampling Veitch et al. (2015).
First, we estimate the joint posterior probability (marginalized over all other parameters of the binary) from the complete observed IMR signal.^{2}^{2}2From here onwards, we drop the explicit reference to the data and the GR model in the posteriors, for simplicity. This allows us to infer the posterior on the mass and dimensionless spin of the final black hole, using fitting formulas (e.g., Healy et al. (2014)) calibrated to NR simulations
(3) 
We use these estimates of and to split the signal into an inspiral part and a merger–ringdown part. In this paper, we define the inspiral [merger–ringdown] part as Fourier frequencies less [greater] than that of the innermost stable circular orbit (ISCO) of a Kerr black hole with mass and spin equal to that given by the median value of .^{3}^{3}3While we split the signal in the Fourier domain, we have checked that almost all the power below [above] our split frequency indeed comes from the early [late] portions of the waveform; the effect of the spectral leakage is negligible. However, this choice is not unique; alternative ways of splitting the signal are possible, and reasonable alternatives do not have a significant effect on the test.
We can now independently estimate the posterior from the inspiral part of the signal and compute the corresponding posterior of the mass and spin of the final black hole using the fitting formula Eq. (3). We independently estimate the posterior from the merger–ringdown part of the signal. In the absence of any deviations from GR (or significant systematic errors), we expect the two posteriors and to overlap (see, e.g., the top left panel in Fig. 1).
To constrain possible departures from GR, we define two parameters that describe departures from the GR prediction of the mass and spin of the final black hole
(4) 
whose posterior distribution can be computed as
(5)  
In the absence of departures from GR, we expect to be consistent with zero. We define two quantities and that describe the fractional differences in the two predictions of the mass and spin of the final black hole. The posteriors on these can be computed as
(6) 
Here, the posterior denotes our best estimate of the mass and spin of the final black hole assuming GR, which is estimated from the full IMR waveform. An example of the posterior distribution from a simulated GR signal is shown in the bottom left panel of Fig. 1. Finally, the posteriors from multiple observations of binary black holes can be combined to construct a single posterior that can better constrain deviations from GR (see, e.g., Fig. 2).
Iii Implementation
To compute the posterior distributions, we employ the LALInference Veitch et al. (2015) stochastic samplers available in the LIGO Algorithm Library (LAL) url . In particular we use the LALInferenceNest code Veitch and Vecchio (2010), which implements a nested sampling algorithm Skilling (2004) in the context of GW data analysis. As the GR signal model we employ the gravitational waveform family SEOBNRv2_ROM_DoubleSpin Pürrer (2016) which describes the inspiral, merger and ringdown waveform of blackhole binaries with nonprecessing spins. This is a reducedorder model version Pürrer (2016) of the effectiveonebody (EOB) waveform family Taracchini et al. (2014) calibrated to NR simulations. We use the fitting formulas proposed in Healy et al. (2014) to compute the mass and spin of the final black hole from the initial masses and (nonprecessing) spins.
From the (simulated) data containing a GW signal, we compute the posterior distributions of and in three different ways:

is computed from the full data: we set and in Eq. (2), where is the lowfrequency cutoff of the detector and is the Nyquist frequency of the data. From the median value of the posterior , we compute the frequency of the Kerr ISCO (). This is used as the characteristic frequency to delineate the inspiral and merger–ringdown parts of the signal in our current analysis.

is computed from the inspiral part of the data: we set and in Eq. (2).

is computed from the merger–ringdown part of the data: we set and in Eq. (2).
All posteriors are computed by assuming a prior distribution that is uniform in and . The posterior is computed from Eq. (5) using SciPy’s correlate2d function and is computed by numerically integrating Eq. (6).
Iv GR simulations
We have performed simulations where we inject simulated GW signals modelling inspiral, merger and ringdown of binary black holes (based on GR, as modelled by SEOBNRv2_ROM_DoubleSpin) into colored Gaussian noise with the design power spectrum of the Advanced LIGO detectors in the highpower, zerodetuning configuration aLI , with a low frequency cutoff Hz. Binaries had component masses (detector frame) uniformly distributed in the range and nonprecessing spins in the range . Sources were distributed uniformly in the sky with isotropic orientations in such a way that the observed signals will have a network SNR of . The estimated posterior from a single simulated event is shown in the bottom left panel of Fig. 1. We also combine posteriors from multiple events; Figure 2 shows the combined posteriors as a function of the number of simulated events. The constraints on the deviation parameters become narrower when multiple events are combined. The width of the credible region could be as low as a few percent when observations are combined. This is within the reach of one year of Advanced LIGO observation, according to several population synthesis models Dominik et al. (2015); Rodriguez et al. (2015).
V Modified GR simulations
We also test our analysis pipeline using simulated GW signals that show departures from GR. To obtain waveforms whose energy and angular momentum loss differs from that predicted by GR, we have chosen to make kludge waveforms based on a simple modification of EOB waveforms. Specifically, we take the IHES EOB waveform model described in Damour et al. (2013), which is given as publicly available code at IHE , and modify the GW flux starting at second postNewtonian (PN) order by multiplying the six modes that first enter at PN [viz., the , , and modes] by a constant factor .^{4}^{4}4The corresponding change in the PN phasing coefficients will depend on the mass ratio. For equal masses, the 2PN term in the frequency domain phase expression will be modified by a factor of . Such a PN modification to the flux is unconstrained by measurements of the GW energy loss from the double pulsar J07373039 Kramer et al. (2006); Yunes and Hughes (2010). We also multiply those modes of the waveform by a consistent factor . However, only the dominant modes are used for simulating the observation. As in the original code, we use the maximum of the orbital frequency (calculated from the EOB Hamiltonian) to mark the termination of the inspiral (and the start of the matching to the quasinormal modes to give the merger and ringdown). The eccentricity of our modified waveforms remains as small () as for the unmodified waveforms.
Since the final mass and spin in the original EOB waveform are set by a fit to NR results, for the modified waveform we replace this determination by demanding selfconsistency of the radiated energy and angular momentum. That is, we choose the final mass and spin by minimizing the difference between the values we set for the final black hole and those obtained by energy and angularmomentum balance using the initial data and the radiated quantities calculated from the waveform (through ). This treatment assumes that the standard GR expressions for the radiated energy and angular momentum remain valid for this modified gravity waveform, which is indeed the case for a large range of modified theories Stein and Yunes (2011). We have not changed the quasinormal mode spectrum of the final black hole, for simplicity.
The right panels of Fig. 1 show the estimated posteriors on the mass and spin of the final black hole from one modified GR simulation (equalmass, nonspinning binary), for which the final mass and spin are and , compared to and in the analogous GR case.^{5}^{5}5The GR waveform used in the left panels of Fig. 1 was computed using the unmodified IHES EOB code, to allow a direct comparison with the modified GR result, though the differences between SEOBNRv2 and IHES EOB are very small for this equalmass, nonspinning case. We also show the posterior on the parameters describing deviations from the GR predictions. It can be seen that the GR value (marked by a “+” sign) is well outside the credible region of the estimated posterior. In this example, GR can be ruled out with confidence. We have verified that this signal, having an optimal SNR of 18.7, produces a chisquare weighted SNR when filtering with the bestfit GR waveform and would thus likely be detected by a standard detection pipeline Usman et al. (2015).
Vi Conclusions
The test that we propose assumes the validity of GR and tests the null hypothesis by computing the posterior distribution for the parameters that quantify a deviation from the result in GR, where both parameters are identically zero. Multiple observations could be combined to produce better constraints on the deviation. We have seen that this test is able to detect deviations from GR that are not constrained by radio observations of the orbital decay of the double pulsar – the tightest constraint available. The test is not based on a specific theory and, consequently, could work in any theory in which massive compact binaries inspiral, merge, and then ringdown. Conversely, if the data were inconsistent with the null hypothesis, then they would not be able to give any direct indication of which modified theory is responsible for the deviation from GR. We expect this test to complement other GWbased tests of GR, including those looking for specific modifications to GR and those looking for generic parametrized deviations, providing confidence in any statements of whether a given signal (or population of signals) is consistent with GR.
Although we have used the ISCO frequency of the final Kerr black hole to delineate between inspiral and merger–ringdown in this paper, alternative ways of splitting the signal are possible. We have verified that the main results are robust against (reasonable) choices of cutoff frequencies. We have neglected the effect of spin precession and subdominant modes in this paper. However, they can be readily included in this method by incorporating these effects in our GR model and also (in the case of precession) in the fitting formulas for the final mass and spin. Systematic errors due to waveform inaccuracies could be mitigated or quantified by using waveform models that are better calibrated to NR simulations as they become available. Methods for mitigating the systematic errors due to detector calibration errors have been independently developed which involve marginalizing the posterior distributions of the masses and spins over additional parameters that model calibration errors Farr et al. . Studies pertaining to these aspects are to be reported in a forthcoming paper Ghosh et al. .
The test introduced in this paper has already had its first application: This was one of the tests used to establish the consistency of LIGO’s first gravitational wave detection with a binary black hole signal as predicted by GR B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration) (2016a, b).
Acknowledgements.
We thank J. Veitch and A. Nagar for assistance with the LALInference and IHES EOB codes, respectively. We also thank K. G. Arun, A. Buonanno, N. Christensen, B. R. Iyer, C. Messenger, A. Mukherjee, B. S. Sathyaprakash and C. Van Den Broeck for useful discussions. 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/PS191/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. Y. C.’s research is supported by NSF Grant PHY1404569, and D.N.’s by NSF grant PHY1404105. C. P. L. B. was supported by the Science and Technology Facilities Council. Computations were performed at the ICTS clusters Mowgli, Dogmatix, and Alice.References
 Centrella et al. (2010) J. Centrella, J. G. Baker, B. J. Kelly, and J. R. van Meter, Rev. Mod. Phys. 82, 3069 (2010).
 Aasi et al. (2015) J. Aasi et al., Class. Quantum Grav. 32, 074001 (2015).
 Acernese et al. (2015) F. Acernese et al., Class. Quantum Grav. 32, 024001 (2015).
 Somiya (2012) K. Somiya, Class. Quantum Grav. 29, 124007 (2012).
 Iyer et al. (2011) B. Iyer et al., LIGOIndia, proposal of the consortium for Indian initiative in gravitationalwave observations (IndIGO) (2011), URL https://dcc.ligo.org/LIGOM1100296/public.
 Abbott et al. (2016) B. P. Abbott et al., Living Rev. Relativity 19, 1 (2016).
 Ajith et al. (2008) P. Ajith et al., Phys. Rev. D 77, 104017 (2008).
 Dominik et al. (2015) M. Dominik et al., Astrophys. J. 806, 263 (2015).
 Rodriguez et al. (2015) C. L. Rodriguez et al., Phys. Rev. Lett. 115, 051101 (2015).
 Punturo et al. (2010) M. Punturo et al., Class. Quantum Grav. 27, 194002 (2010).
 Miao et al. (2014) H. Miao et al., Class. Quantum Grav. 31, 165010 (2014).
 Dwyer et al. (2015) S. Dwyer et al., Phys. Rev. D 91, 082001 (2015).
 Berti et al. (2015) E. Berti et al., Class. Quantum Grav. 32, 243001 (2015).
 Yunes and Siemens (2013) N. Yunes and X. Siemens, Living Rev. Rel. 16, 9 (2013).
 Hughes and Menou (2005) S. A. Hughes and K. Menou, Astrophys. J. 623, 689 (2005).
 Nakano et al. (2015) H. Nakano, T. Tanaka, and T. Nakamura, Phys. Rev. D 92, 064003 (2015).
 Veitch et al. (2015) J. Veitch et al., Phys. Rev. D 91, 042003 (2015).
 Healy et al. (2014) J. Healy, C. O. Lousto, and Y. Zlochower, Phys. Rev. D 90, 104004 (2014).
 (19) http://www.lscgroup.phys.uwm.edu/lal.
 Veitch and Vecchio (2010) J. Veitch and A. Vecchio, Phys. Rev. D 81 (2010).
 Skilling (2004) J. Skilling, in American Institute of Physics Conference Series (2004), vol. 735, pp. 395–405.
 Pürrer (2016) M. Pürrer, Phys. Rev. D 93, 064041 (2016).
 Taracchini et al. (2014) A. Taracchini et al., Phys. Rev. D 89, 061502(R) (2014).
 (24) Advanced LIGO anticipated sensitivity curves, URL https://dcc.ligo.org/LIGOT0900288/public.
 Damour et al. (2013) T. Damour, A. Nagar, and S. Bernuzzi, Phys. Rev. D 87, 084035 (2013).
 (26) http://eob.ihes.fr; we used the “1202” version of the code.
 Kramer et al. (2006) M. Kramer et al., Science 314, 97 (2006).
 Yunes and Hughes (2010) N. Yunes and S. A. Hughes, Phys. Rev. D 82, 082002 (2010).
 Stein and Yunes (2011) L. C. Stein and N. Yunes, Phys. Rev. D 83, 064038 (2011).
 Usman et al. (2015) S. A. Usman et al. (2015), 1508.02357.
 (31) W. M. Farr, B. Farr, and T. Littenberg, URL https://dcc.ligo.org/LIGOT1400682/public.
 (32) A. Ghosh et al., in preparation.
 B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration) (2016a) B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 061102 (2016a).
 B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration) (2016b) B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 116, 221101 (2016b).