Tests of General Relativity with the Binary Black Hole Signals from the LIGO-Virgo Catalog GWTC-1

TL;DR

This paper tests general relativity using gravitational wave data from binary black hole mergers, finding no deviations and setting new bounds on graviton mass and alternative polarization constraints.

Contribution

It presents four novel tests of general relativity with LIGO-Virgo data, including residual analysis, frequency consistency, waveform deviations, and gravitational wave dispersion constraints.

Findings

No inconsistencies with general relativity found.

Improved bounds on graviton mass to 4.7 x 10^{-23} eV/c^2.

First constraints on alternative polarizations from GWTC-1 events.

Abstract

The detection of gravitational waves by Advanced LIGO and Advanced Virgo provides an opportunity to test general relativity in a regime that is inaccessible to traditional astronomical observations and laboratory tests. We present four tests of the consistency of the data with binary black hole gravitational waveforms predicted by general relativity. One test subtracts the best-fit waveform from the data and checks the consistency of the residual with detector noise. The second test checks the consistency of the low- and high-frequency parts of the observed signals. The third test checks that phenomenological deviations introduced in the waveform model (including in the post-Newtonian coefficients) are consistent with zero. The fourth test constrains modifications to the propagation of gravitational waves due to a modified dispersion relation, including that from a massive graviton. We present results both for individual events and also results obtained by combining together particularly strong events from the first and second observing runs of Advanced LIGO and Advanced Virgo, as collected in the catalog GWTC-1. We do not find any inconsistency of the data with the predictions of general relativity and improve our previously presented combined constraints by factors of 1.1 to 2.5. In particular, we bound the mass of the graviton to be $m_g \leq 4.7 \times 10^{-23} \text{eV}/c^2$ ($90\%$ credible level), an improvement of a factor of 1.6 over our previously presented results. Additionally, we check that the four gravitational-wave events published for the first time in GWTC-1 do not lead to stronger constraints on alternative polarizations than those published previously.

Tables

Table 1. Table 1: The GW events considered in this paper, separated by observing run. The first block of columns gives the names of the events and lists some of their relevant properties obtained using GR waveforms (luminosity distance D L subscript 𝐷 L D_{\text{L}} , source frame total mass M tot subscript 𝑀 tot M_{\text{tot}} and final mass M f subscript 𝑀 f M_{\text{f}} , and dimensionless final spin a f subscript 𝑎 f a_{\text{f}} ). The next block of columns gives the significance, measured by the false-alarm-rate (FAR), with which each event was detected by each of the three searches employed, as well as the matched filter signal-to-noise ratio from the stochastic sampling analyses with GR waveforms. A dash indicates that an event was not identified by a search. The parameters and SNR values give the medians and 90% credible intervals. All the events except for GW151226 and GW170729 are consistent with a binary of nonspinning black holes (when analyzed assuming GR). See  Abbott  et al. ( 2019 ) for more details about all the events. The last block of columns indicates which GR tests are performed on a given event: RT = residuals test (Sec.  V.1 ); IMR = inspiral-merger-ringdown consistency test (Sec.  V.2 ); PI & PPI = parameterized tests of GW generation for inspiral and post-inspiral phases (Sec.  VI ); MDR = modified GW dispersion relation (Sec.  VII ). The events with bold names are used to obtain the combined results for each test.

Event	Properties					FAR				SNR		GR tests performed
	$D_{\text{L}}$	$M_{\text{tot}}$	$M_{\text{f}}$	$a_{\text{f}}$		PyCBC	GstLAL	cWB				RT	IMR	PI	PPI	MDR
	[Mpc]	[$M_{\odot }$]	[$M_{\odot }$]			[yr-1]	[yr-1]	[yr-1]
GW150914555 The FARs for these events differ from those in Abbott et al. (2016) because the data were re-analyzed with the new pipeline statistics used in O2 (see Abbott et al. (2019) for more details).,,	${440}_{-170}^{+150}$	${66.1}_{-3.3}^{+3.8}$	${63.1}_{-3.0}^{+3.4}$	${0.69}_{-0.04}^{+0.05}$						${25.3}_{-0.2}^{+0.1}$		✓	✓	✓	✓	✓
GW1510125,,	${1080}_{-490}^{+550}$	${37.2}_{-3.9}^{+10.6}$	${35.6}_{-3.8}^{+10.8}$	${0.67}_{-0.11}^{+0.13}$		$0.17$	$7.9\times {10}^{-3}$	–		${9.2}_{-0.4}^{+0.3}$		✓	–	–	✓	✓
GW1512265,666At least one black hole has dimensionless spin $>0.28$ ($99\%$ credible level).	${450}_{-190}^{+180}$	${21.5}_{-1.5}^{+6.2}$	${20.5}_{-1.5}^{+6.4}$	${0.74}_{-0.05}^{+0.07}$				$0.02$		${12.4}_{-0.3}^{+0.2}$		✓	–	✓	–	✓
GW170104,,,,	${990}_{-430}^{+440}$	${51.0}_{-4.1}^{+5.3}$	${48.9}_{-4.0}^{+5.1}$	${0.66}_{-0.11}^{+0.08}$				$2.9\times {10}^{-4}$		${14.0}_{-0.3}^{+0.2}$		✓	✓	✓	✓	✓
GW170608,,,,	${320}_{-110}^{+120}$	${18.6}_{-0.7}^{+3.2}$	${17.8}_{-0.7}^{+3.4}$	${0.69}_{-0.04}^{+0.04}$				$1.4\times {10}^{-4}$		${15.6}_{-0.3}^{+0.2}$		✓	–	✓	✓	✓
GW170729777This event has a higher significance in the unmodeled search than in the modeled searches. Additionally, at least one black hole has dimensionless spin $>0.27$ ($99\%$ credible level).,,	${2840}_{-1360}^{+1400}$	${84.4}_{-11.1}^{+15.8}$	${79.5}_{-10.2}^{+14.7}$	${0.81}_{-0.13}^{+0.07}$		$1.4$	$0.18$	$0.02$		${10.8}_{-0.5}^{+0.4}$		✓	✓	–	✓	✓
GW170809,,,,	${1030}_{-390}^{+320}$	${59.0}_{-4.1}^{+5.4}$	${56.3}_{-3.8}^{+5.2}$	${0.70}_{-0.09}^{+0.08}$		$1.4\times {10}^{-4}$		–		${12.7}_{-0.3}^{+0.2}$		✓	✓	–	✓	✓
GW170814,,,,	${600}_{-220}^{+150}$	${55.9}_{-2.6}^{+3.4}$	${53.2}_{-2.4}^{+3.2}$	${0.72}_{-0.05}^{+0.07}$						${17.8}_{-0.3}^{+0.3}$		✓	✓	✓	✓	✓
GW170818,,,,	${1060}_{-380}^{+420}$	${62.2}_{-4.1}^{+5.2}$	${59.4}_{-3.8}^{+4.9}$	${0.67}_{-0.08}^{+0.07}$		–	$4.2\times {10}^{-5}$	–		${11.9}_{-0.4}^{+0.3}$		✓	✓	–	✓	✓
GW170823,,,,	${1940}_{-900}^{+970}$	${68.7}_{-8.1}^{+10.8}$	${65.4}_{-7.4}^{+10.1}$	${0.72}_{-0.12}^{+0.09}$				$2.1\times {10}^{-3}$		${12.0}_{-0.3}^{+0.2}$		✓	✓	–	✓	✓

Table 2. Table 2: Results of the residuals analysis. For each event, this table presents the 90%-credible upper limit on the reconstructed network SNR after subtraction of the best-fit GR waveform ( SNR 90 subscript SNR 90 \mathrm{SNR}_{90} ), a corresponding lower limit on the fitting factor ( FF 90 subscript FF 90 \mathrm{FF}_{90} in the text), and the SNR 90 subscript SNR 90 \mathrm{SNR}_{90} p 𝑝 p -value. SNR 90 subscript SNR 90 \mathrm{SNR}_{90} is a measure of the maximum possible coherent signal power not captured by the best-fit GR template, while the p 𝑝 p -value is an estimate of the probability that instrumental noise produced such SNR 90 subscript SNR 90 \mathrm{SNR}_{90} or higher. We also indicate which interferometers (IFOs) were used in the analysis of a given event, either the two Advanced LIGO detectors (HL) or the two Advanced LIGO detectors plus Advanced Virgo (HLV). See Sec.  V.1 in the main text for details.

Event	IFOs	Residual ${\mathrm{SNR}}_{90}$	Fitting factor	$p$-value
GW150914	HL	6.1	$\ge 0.97$	0.46
GW151012	HL	7.3	$\ge 0.79$	0.11
GW151226	HL	5.6	$\ge 0.91$	0.81
GW170104	HL	5.1	$\ge 0.94$	0.99
GW170608	HL	7.9	$\ge 0.89$	0.05
GW170729	HLV	6.5	$\ge 0.85$	0.74
GW170809	HLV	6.5	$\ge 0.88$	0.78
GW170814	HLV	8.9	$\ge 0.88$	0.16
GW170818	HLV	9.2	$\ge 0.78$	0.19
GW170823	HL	5.5	$\ge 0.90$	0.86

Table 3. Table 3: Results from the inspiral-merger-ringdown consistency test for selected binary black hole events. f c subscript 𝑓 c f_{\text{c}} denotes the cutoff frequency used to demarcate the division between the inspiral and post-inspiral regimes; ρ IMR subscript 𝜌 IMR \rho_{\mathrm{IMR}} , ρ insp subscript 𝜌 insp \rho_{\mathrm{insp}} , and ρ post − insp subscript 𝜌 post insp \rho_{\mathrm{post-insp}} are the median values of the SNR in the full signal, the inspiral part, and the post-inspiral part, respectively; and the GR quantile denotes the fraction of the posterior enclosed by the isoprobability contour that passes through the GR value, with smaller values indicating better consistency with GR. (Note, however, that the posterior distribution will be broader for smaller SNRs, and hence the GR quantile will typically be smaller in such cases. This effect is further complicated by the randomness of the noise.)

Event	$f_{\text{c}}$ [Hz]	${\rho }_{\mathrm{IMR}}$	${\rho }_{\mathrm{insp}}$	${\rho }_{\mathrm{post}-\mathrm{insp}}$	GR quantile [%]
GW150914	132	25.3	19.4	16.1	55.5
GW170104	143	13.7	10.9	8.5	24.4
GW170729	91	10.7	8.6	6.9	10.4
GW170809	136	12.7	10.6	7.1	14.7
GW170814	161	16.8	15.3	7.2	7.8
GW170818	128	12.0	9.3	7.2	25.5
GW170823	102	11.9	7.9	8.5	80.4

Table 4. Table 4: 90% credible level upper bounds on the graviton mass m g subscript 𝑚 𝑔 m_{g} and the absolute value of the modified dispersion relation parameter A α subscript 𝐴 𝛼 A_{\alpha} , as well as the GR quantiles Q GR subscript 𝑄 GR Q_{\text{GR}} . The < < and > > labels denote the bounds for A α < 0 subscript 𝐴 𝛼 0 A_{\alpha}<0 and > 0 absent 0 >0 , respectively, and we have defined the dimensionless quantity A ¯ α ≔ A α / eV 2 − α ≔ subscript ¯ 𝐴 𝛼 subscript 𝐴 𝛼 superscript eV 2 𝛼 \bar{A}_{\alpha}\coloneqq A_{\alpha}/\text{eV}^{2-\alpha} . Events with names in boldface are used to obtain the combined results.

	$m_g$		$\left|{\overline{A}}_0\right|$				$\left|{\overline{A}}_{0.5}\right|$				$\left|{\overline{A}}_1\right|$				$\left|{\overline{A}}_{1.5}\right|$				$\left|{\overline{A}}_{2.5}\right|$				$\left|{\overline{A}}_3\right|$				$\left|{\overline{A}}_{3.5}\right|$				$\left|{\overline{A}}_4\right|$
Event	[${10}^{-23}$			$>$	$Q_{\text{GR}}$			$>$	$Q_{\text{GR}}$			$>$	$Q_{\text{GR}}$			$>$	$Q_{\text{GR}}$			$>$	$Q_{\text{GR}}$			$>$	$Q_{\text{GR}}$			$>$	$Q_{\text{GR}}$			$>$	$Q_{\text{GR}}$
	eV$/c^2$]		[${10}^{-44}$]		[%]		[${10}^{-38}$]		[%]		[${10}^{-32}$]		[%]		[${10}^{-25}$]		[%]		[${10}^{-13}$]		[%]		[${10}^{-8}$]		[%]		[${10}^{-2}$]		[%]		[${10}^4$]		[%]
GW150914	$9.9$		$1.4$	$1.1$	$71$		$5.9$	$4.9$	$57$		$5.5$	$3.5$	$74$		$3.2$	$2.1$	$74$		$2.4$	$2.1$	$50$		$17$	$19$	$37$		$11$	$20$	$42$		$7.7$	$9.8$	$52$
GW151012	$17$		$3.8$	$3.5$	$40$		$3.6$	$11$	$35$		$6.6$	$9.5$	$41$		$1.9$	$2.5$	$46$		$2.8$	$1.5$	$56$		$21$	$9.7$	$56$		$18$	$15$	$61$		$14$	$6.7$	$58$
GW151226	$29$		$7.1$	$9.3$	$21$		$8.1$	$21$	$10$		$13$	$13$	$26$		$3.4$	$3.3$	$28$		$3.4$	$2.1$	$68$		$22$	$9.7$	$61$		$14$	$6.1$	$58$		$19$	$5.4$	$72$
GW170104	$9.2$		$2.5$	$0.99$	$62$		$4.2$	$2.4$	$72$		$7.0$	$2.9$	$76$		$1.9$	$0.80$	$83$		$1.1$	$2.8$	$23$		$10$	$13$	$35$		$6.4$	$10$	$37$		$6.5$	$8.7$	$43$
GW170608	$30$		$13$	$9.2$	$49$		$22$	$8.8$	$68$		$15$	$28$	$49$		$3.1$	$3.8$	$64$		$3.3$	$2.1$	$49$		$9.8$	$8.3$	$46$		$87$	$8.0$	$46$		$30$	$3.7$	$38$
GW170729	$7.4$		$0.29$	$0.64$	$17$		$0.93$	$1.1$	$26$		$2.1$	$4.5$	$16$		$0.79$	$1.5$	$18$		$4.2$	$1.4$	$94$		$36$	$8.6$	$96$		$26$	$8.9$	$94$		$43$	$7.1$	$95$
GW170809	$9.3$		$1.4$	$1.2$	$64$		$2.5$	$2.5$	$49$		$11$	$5.8$	$49$		$1.4$	$2.3$	$37$		$3.8$	$1.3$	$79$		$22$	$6.9$	$81$		$18$	$6.2$	$78$		$13$	$6.5$	$82$
GW170814	$8.5$		$4.0$	$1.1$	$94$		$5.2$	$1.8$	$92$		$15$	$5.2$	$93$		$3.0$	$1.3$	$96$		$1.3$	$3.9$	$5.7$		$7.4$	$25$	$6.2$		$5.4$	$29$	$7.5$		$4.2$	$13$	$9.6$
GW170818	$7.2$		$1.4$	$0.66$	$74$		$2.5$	$1.9$	$80$		$4.9$	$4.3$	$73$		$1.7$	$0.67$	$79$		$1.4$	$3.4$	$28$		$17$	$15$	$41$		$19$	$7.9$	$73$		$12$	$9.0$	$49$
GW170823	$6.3$		$1.2$	$0.49$	$61$		$1.1$	$1.6$	$51$		$2.7$	$2.3$	$49$		$0.99$	$1.2$	$46$		$2.5$	$1.3$	$53$		$11$	$16$	$41$		$8.8$	$11$	$37$		$14$	$11$	$46$
Combined	$4.7$		$0.80$	$0.34$	$79$		$1.2$	$0.70$	$73$		$2.5$	$1.2$	$70$		$0.70$	$0.37$	$86$		$0.50$	$0.80$	$28$		$2.9$	$3.7$	$25$		$2.0$	$3.7$	$35$		$1.4$	$2.3$	$34$

Table 5. Table 5: Median value and symmetric 90% credible level bounds of the waveform parameters δ ​ p ^ i 𝛿 subscript ^ 𝑝 𝑖 \delta\hat{p}_{i} , as well as the GR quantiles Q GR subscript 𝑄 GR Q_{\text{GR}} . For the inspiral parameters, we show results in pairs of rows for when the data from individual events are analyzed using IMRPhenomPv2 (P) and SEOBNRv4 (S), while for the post-inspiral phenomenological parameters, results are obtained only for IMRPhenomPv2 .

GW150914

GW151226

GW170104

GW170608

GW170814

Parameter

Model

${\tilde{X}}_{-}^{+}$

$Q_{\text{GR}}$ [%]

${\tilde{X}}_{-}^{+}$

$Q_{\text{GR}}$ [%]

${\tilde{X}}_{-}^{+}$

$Q_{\text{GR}}$ [%]

${\tilde{X}}_{-}^{+}$

$Q_{\text{GR}}$ [%]

${\tilde{X}}_{-}^{+}$

$Q_{\text{GR}}$ [%]

$\delta {\hat{\varphi }}_{-2}$

$\left[{10}^{-2}\right]$

P

$-{0.46}_{-0.59}^{+0.68}$

$88$

${0.14}_{-0.35}^{+0.52}$

$27$

${2.0}_{-1.9}^{+3.6}$

$3$

$-{0.06}_{-0.35}^{+0.61}$

$60$

$-{0.76}_{-0.81}^{+0.74}$

$96$

S

$-{0.63}_{-0.61}^{+0.62}$

$95$

${0.10}_{-0.31}^{+0.40}$

$30$

${1.4}_{-1.4}^{+2.6}$

$5.9$

${0.10}_{-0.34}^{+0.27}$

$26$

$-{0.99}_{-0.79}^{+0.89}$

$97$

$\delta {\hat{\varphi }}_0$

$\left[{10}^{-1}\right]$

P

$-{1.0}_{-1.1}^{+1.0}$

$95$

$-{0.1}_{-1.7}^{+1.7}$

$53$

${2.5}_{-2.3}^{+4.3}$

$3.1$

$-{0.5}_{-1.6}^{+1.4}$

$73$

$-{0.8}_{-1.2}^{+1.2}$

$88$

S

$-{1.1}_{-1.1}^{+1.0}$

$97$

${0.3}_{-2.8}^{+2.4}$

$41$

${1.8}_{-2.0}^{+3.5}$

$7.3$

$-{0.9}_{-2.2}^{+1.8}$

$77$

$-{1.1}_{-2.0}^{+1.5}$

$89$

$\delta {\hat{\varphi }}_1$

$\left[{10}^0\right]$

P

$-{0.39}_{-0.50}^{+0.40}$

$95$

${0.03}_{-0.44}^{+0.43}$

$45$

${1.0}_{-0.9}^{+1.7}$

$4.1$

$-{0.14}_{-0.33}^{+0.33}$

$76$

$-{0.31}_{-0.50}^{+0.47}$

$86$

S

$-{0.51}_{-0.41}^{+0.41}$

$98$

${0.02}_{-0.56}^{+0.55}$

$48$

${0.7}_{-0.8}^{+1.2}$

$6.6$

$-{0.16}_{-0.40}^{+0.38}$

$77$

$-{0.27}_{-0.61}^{+0.57}$

$79$

$\delta {\hat{\varphi }}_2$

$\left[{10}^0\right]$

P

$-{0.35}_{-0.31}^{+0.32}$

$97$

$-{0.01}_{-0.24}^{+0.29}$

$52$

${0.7}_{-0.6}^{+1.1}$

$3.5$

$-{0.07}_{-0.20}^{+0.21}$

$72$

$-{0.17}_{-0.36}^{+0.34}$

$80$

S

$-{0.34}_{-0.30}^{+0.28}$

$97$

${0.02}_{-0.32}^{+0.33}$

$47$

${0.47}_{-0.56}^{+0.83}$

$8.6$

$-{0.09}_{-0.25}^{+0.24}$

$75$

$-{0.10}_{-0.38}^{+0.39}$

$67$

$\delta {\hat{\varphi }}_3$

$\left[{10}^{-1}\right]$

P

${2.2}_{-1.9}^{+2.0}$

$3.4$

$-{0.1}_{-2.0}^{+1.5}$

$54$

$-{4.8}_{-6.5}^{+4.2}$

$97$

${0.5}_{-1.2}^{+1.2}$

$26$

${0.7}_{-2.2}^{+2.1}$

$30$

S

${2.2}_{-1.8}^{+2.0}$

$1.8$

$-{0.2}_{-2.0}^{+2.0}$

$55$

$-{3.2}_{-5.2}^{+3.8}$

$91$

${0.6}_{-1.4}^{+1.6}$

$26$

${0.1}_{-2.4}^{+2.4}$

$48$

$\delta {\hat{\varphi }}_4$

$\left[{10}^0\right]$

P

$-{1.9}_{-1.6}^{+1.7}$

$97$

${0.1}_{-1.3}^{+1.6}$

$47$

${3.7}_{-3.6}^{+4.6}$

$4.2$

$-{0.3}_{-1.1}^{+1.1}$

$67$

$-{0.5}_{-1.6}^{+1.7}$

$67$

S

$-{1.7}_{-1.5}^{+1.3}$

$98$

${0.2}_{-1.6}^{+1.8}$

$41$

${2.2}_{-2.9}^{+3.9}$

$10$

$-{0.4}_{-1.4}^{+1.3}$

$70$

${0.1}_{-1.7}^{+1.9}$

$45$

$\delta {\hat{\varphi }}_5^{\left(l\right)}$

$\left[{10}^0\right]$

P

${0.70}_{-0.58}^{+0.56}$

$2.2$

$-{0.03}_{-0.67}^{+0.49}$

$54$

$-{1.4}_{-1.6}^{+1.3}$

$97$

${0.09}_{-0.42}^{+0.42}$

$36$

${0.10}_{-0.62}^{+0.69}$

$40$

S

${0.61}_{-0.45}^{+0.49}$

$1.5$

${0.00}_{-0.77}^{+0.71}$

$50$

$-{1.0}_{-1.5}^{+1.1}$

$92$

${0.25}_{-0.54}^{+0.62}$

$23$

$-{0.16}_{-0.72}^{+0.66}$

$64$

$\delta {\hat{\varphi }}_6$

$\left[{10}^0\right]$

P

$-{1.2}_{-1.1}^{+1.0}$

$97$

${0.2}_{-1.1}^{+1.3}$

$40$

${2.6}_{-2.5}^{+2.9}$

$4.7$

$-{0.11}_{-0.84}^{+0.84}$

$59$

$-{0.1}_{-1.1}^{+1.2}$

$56$

S

$-{1.07}_{-0.92}^{+0.80}$

$99$

${0.1}_{-1.3}^{+1.3}$

$45$

${1.6}_{-2.1}^{+2.7}$

$11$

$-{0.3}_{-1.0}^{+1.0}$

$68$

${0.4}_{-1.1}^{+1.3}$

$29$

$\delta {\hat{\varphi }}_6^{\left(l\right)}$

$\left[{10}^1\right]$

P

${0.78}_{-0.70}^{+0.63}$

$3.5$

$-{0.05}_{-0.66}^{+0.56}$

$55$

$-{1.6}_{-1.9}^{+1.5}$

$96$

${0.09}_{-0.47}^{+0.44}$

$36$

${0.20}_{-0.78}^{+0.73}$

$32$

S

${0.72}_{-0.55}^{+0.57}$

$1.8$

$-{0.13}_{-0.76}^{+0.70}$

$61$

$-{1.0}_{-1.6}^{+1.1}$

$92$

${0.15}_{-0.52}^{+0.57}$

$31$

$-{0.07}_{-0.78}^{+0.69}$

$56$

$\delta {\hat{\varphi }}_7$

$\left[{10}^0\right]$

P

${3.2}_{-2.6}^{+2.7}$

$2.2$

$-{0.4}_{-3.5}^{+2.8}$

$58$

$-{7.0}_{-7.4}^{+6.9}$

$95$

${0.3}_{-2.2}^{+2.5}$

$42$

${0.2}_{-3.5}^{+3.2}$

$47$

S

${2.8}_{-2.0}^{+2.1}$

$1.3$

$-{0.1}_{-3.1}^{+1.7}$

$54$

$-{4.8}_{-7.2}^{+5.4}$

$93$

${1.1}_{-2.7}^{+1.8}$

$23$

$-{1.2}_{-3.2}^{+3.1}$

$73$

$\delta {\hat{\beta }}_2$

$\left[{10}^0\right]$

P

${0.03}_{-0.28}^{+0.35}$

$43$

$-$

$-$

$-{0.23}_{-0.35}^{+0.38}$

$85$

${0.13}_{-0.30}^{+0.40}$

$25$

$-{0.04}_{-0.23}^{+0.28}$

$61$

$\delta {\hat{\beta }}_3$

$\left[{10}^0\right]$

P

${0.00}_{-0.30}^{+0.36}$

$50$

$-$

$-$

$-{0.13}_{-0.41}^{+0.44}$

$71$

${0.12}_{-0.26}^{+0.28}$

$21$

${0.07}_{-0.29}^{+0.35}$

$35$

$\delta {\hat{\alpha }}_2$

$\left[{10}^0\right]$

P

$-{0.02}_{-0.42}^{+0.34}$

$54$

$-$

$-$

$-{0.13}_{-0.58}^{+0.67}$

$64$

$3_{-55}^{+62}$

$40$

$-{0.64}_{-0.59}^{+0.57}$

$97$

$\delta {\hat{\alpha }}_3$

$\left[{10}^0\right]$

P

$-{0.1}_{-1.4}^{+1.7}$

$55$

$-$

$-$

$-{0.3}_{-3.8}^{+3.8}$

$55$

$-1_{-89}^{+95}$

$52$

$-{1.9}_{-1.8}^{+2.1}$

$93$

$\delta {\hat{\alpha }}_4$

$\left[{10}^0\right]$

P

$-{0.01}_{-0.51}^{+0.54}$

$52$

$-$

$-$

${0.5}_{-1.3}^{+2.8}$

$28$

$-1_{-41}^{+45}$

$51$

$-{0.6}_{-1.5}^{+0.9}$

$85$

GW151012

GW170729

GW170809

GW170818

GW170823

Parameter

Model

${\tilde{X}}_{-}^{+}$

$Q_{\text{GR}}$ [%]

${\tilde{X}}_{-}^{+}$

$Q_{\text{GR}}$ [%]

${\tilde{X}}_{-}^{+}$

$Q_{\text{GR}}$ [%]

${\tilde{X}}_{-}^{+}$

$Q_{\text{GR}}$ [%]

${\tilde{X}}_{-}^{+}$

$Q_{\text{GR}}$ [%]

$\delta {\hat{\beta }}_2$

$\left[{10}^0\right]$

P

${0.3}_{-0.6}^{+3.8}$

$24$

${0.01}_{-0.51}^{+0.56}$

$49$

$-{0.13}_{-0.33}^{+0.37}$

$74$

$-{0.32}_{-0.29}^{+0.44}$

$89$

${0.00}_{-0.38}^{+0.49}$

$50$

$\delta {\hat{\beta }}_3$

$\left[{10}^0\right]$

P

${3.3}_{-3.4}^{+3.6}$

$5.6$

${0.31}_{-0.63}^{+0.81}$

$21$

$-{0.20}_{-0.36}^{+0.42}$

$80$

$-{0.41}_{-0.36}^{+0.58}$

$90$

$-{0.02}_{-0.57}^{+0.74}$

$52$

$\delta {\hat{\alpha }}_2$

$\left[{10}^0\right]$

P

$0_{-3}^{+14}$

$37$

${0.78}_{-0.73}^{+0.85}$

$4$

${0.42}_{-0.64}^{+0.69}$

$14$

${0.01}_{-0.65}^{+0.71}$

$49$

$-{0.04}_{-0.66}^{+0.70}$

$55$

$\delta {\hat{\alpha }}_3$

$\left[{10}^0\right]$

P

$5_{-33}^{+81}$

$29$

${4.4}_{-3.8}^{+4.6}$

$2.1$

${2.2}_{-3.0}^{+4.0}$

$12$

${0.6}_{-2.9}^{+3.5}$

$38$

${0.1}_{-2.3}^{+3.5}$

$49$

$\delta {\hat{\alpha }}_4$

$\left[{10}^0\right]$

P

$6_{-30}^{+24}$

$33$

${1.3}_{-1.1}^{+1.2}$

$2.2$

${0.9}_{-1.0}^{+1.1}$

$7.6$

${0.4}_{-1.1}^{+1.0}$

$26$

$-{0.09}_{-0.92}^{+0.97}$

$56$

Table 6. Table 6: 90% credible level upper bounds on the absolute value of the modified dispersion relation parameter A α subscript 𝐴 𝛼 A_{\alpha} , as well as the GR quantiles Q GR subscript 𝑄 GR Q_{\text{GR}} for IMRPhenomPv2 (P) and SEOBNRv4 (S) runs. The < < and > > labels denote the bounds for A α < 0 subscript 𝐴 𝛼 0 A_{\alpha}<0 and > 0 absent 0 >0 , respectively, with the given scalings and A ¯ α ≔ A α / eV 2 − α ≔ subscript ¯ 𝐴 𝛼 subscript 𝐴 𝛼 superscript eV 2 𝛼 \bar{A}_{\alpha}\coloneqq A_{\alpha}/\text{eV}^{2-\alpha} .

			GW170729							GW170814
Quantity					$>$		$Q_{\text{GR}}$ [%]					$>$		$Q_{\text{GR}}$ [%]
			P	S	P	S	P	S		P	S	P	S	P	S
${\overline{A}}_0$	[${10}^{-44}$]		$0.29$	$0.67$	$0.64$	$0.67$	$17$	$39$		$4.0$	$2.7$	$1.1$	$1.4$	$94$	$87$
${\overline{A}}_{0.5}$	[${10}^{-38}$]		$0.93$	$0.86$	$1.1$	$1.8$	$26$	$16$		$5.2$	$7.6$	$1.8$	$1.8$	$92$	$92$
${\overline{A}}_1$	[${10}^{-32}$]		$2.1$	$4.6$	$4.5$	$4.2$	$16$	$26$		$15$	$12$	$5.2$	$4.0$	$93$	$95$
${\overline{A}}_{1.5}$	[${10}^{-25}$]		$0.79$	$0.74$	$1.5$	$1.2$	$18$	$17$		$3.0$	$2.8$	$1.3$	$1.0$	$96$	$94$
${\overline{A}}_{2.5}$	[${10}^{-13}$]		$4.2$	$4.7$	$1.4$	$1.6$	$94$	$92$		$1.3$	$1.2$	$3.9$	$4.7$	$5.7$	$5.1$
${\overline{A}}_3$	[${10}^{-8}$]		$36$	$30$	$8.6$	$16$	$96$	$92$		$7.4$	$8.3$	$25$	$36$	$6.2$	$6.4$
${\overline{A}}_{3.5}$	[${10}^{-2}$]		$26$	$32$	$8.9$	$12$	$94$	$94$		$5.4$	$5.9$	$29$	$17$	$7.5$	$8.5$
${\overline{A}}_4$	[${10}^4$]		$43$	$47$	$7.1$	$9.1$	$95$	$95$		$4.2$	$4.7$	$13$	$20$	$9.6$	$11$

Equations

$$\tilde{h}(f)=A(f)e^{i\Phi (f)}.$$

$$\{\delta {\hat{\phi }}_{-2},\delta {\hat{\phi }}_0,\delta {\hat{\phi }}_1,\delta {\hat{\phi }}_2,\delta {\hat{\phi }}_3,\delta {\hat{\phi }}_4,\delta {\hat{\phi }}_{5l},\delta {\hat{\phi }}_6,\delta {\hat{\phi }}_{6l},\delta {\hat{\phi }}_7\}.$$

$$E^2=p^2{c}^2+A_{\alpha }{p}^{\alpha }{c}^{\alpha }.$$

$$\delta {\Phi }_{\alpha }(f)=\mathrm{sign}(A_{\alpha })\{\begin{array}{ll}\frac{\pi D_{\text{L}}}{\alpha -1}{\lambda }_{A,\text{eff}}^{\alpha -2}{\left(\frac{f}{c}\right)}^{\alpha -1},& \alpha \ne 1\\ \frac{\pi D_{\text{L}}}{{\lambda }_{A,\text{eff}}}\ln \left(\frac{\pi G{\mathcal{M}}^{\text{det}}f}{c^3}\right),& \alpha =1\end{array}.$$

$${\lambda }_{A,\text{eff}}:={\left[\frac{(1+z{)}^{1-\alpha }{D}_{\text{L}}}{D_{\alpha }}\right]}^{1/(\alpha -2)}{\lambda }_A.$$

$$D_{\alpha }=\frac{(1+z{)}^{1-\alpha }}{H_0}{\int }_0^z\frac{(1+\bar{z}{)}^{\alpha -2}}{\sqrt{{\Omega }_{\text{m}}(1+\bar{z}{)}^3+{\Omega }_{\Lambda }}}d\bar{z},$$

Full text

I Introduction

Einstein’s theory of gravity, general relativity (GR), has withstood a large number of experimental tests Will (2014). With the advent of gravitational-wave (GW) astronomy and the observations by the Advanced LIGO Aasi et al. (2015) and Advanced Virgo Acernese et al. (2015) detectors, a range of new tests of GR have become possible. These include both weak field tests of the propagation of GWs, as well as tests of the strong field regime of compact binary sources. See Abbott et al. (2016, 2016, 2017, 2017, 2019) for previous applications of such tests to GW data.

We report results from tests of GR on all the confident binary black hole GW events in the catalog GWTC-1 LIGO Scientific Collaboration and Virgo Collaboration (2018), i.e., those from the first and second observing runs of the advanced generation of detectors. Besides all of the events previously announced (GW150914, GW151012, GW151226, GW170104, GW170608, and GW170814) Abbott et al. (2016, 2016, 2016, 2016, 2017, 2017, 2017), this includes the four new GW events reported in Abbott et al. (2019) (GW170729, GW170809, GW170818, and GW170823). We do not investigate any of the marginal triggers in GWTC-1, which have a false-alarm rate (FAR) greater than one per year. Table 1 displays a complete list of the events we consider. Tests of GR on the binary neutron star event GW170817 are described in Abbott et al. (2019).

The search results in Abbott et al. (2019) originate from two modeled searches and one weakly modeled search Abbott et al. (2016, 2016, 2016, 2019). The modeled searches use templates based on GR to find candidate events and to assess their significance. However, detection by such searches does not in itself imply full compatibility of the signal with GR Vallisneri and Yunes (2013); Vitale and Del Pozzo (2014). The weakly modeled search relies on coherence of signals between multiple detectors, as expected for an astrophysical source. While it assumes that the morphology of the signal resembles a chirp (whose frequency increases with time), as expected for a compact binary coalescence, it does not assume that the detailed waveform shape agrees with GR. A transient signal strongly deviating from GR would likely be found by the weakly modeled search, even if missed by the modeled searches. So far, however, all significant [$\mathrm{FAR}<(1\,\mathrm{yr})^{-1}$] transient signals found by the weakly modeled search were also found by at least one of the modeled searches Abbott et al. (2019).

At present, there are no complete theories of gravity other than GR that are mathematically and physically viable and provide well defined alternative predictions for the waveforms arising from the coalescence of two black holes (if, indeed, these theories even admit black holes).111There are very preliminary simulations of scalar waveforms from binary black holes in the effective field theory (EFT) framework in alternative theories Okounkova et al. (2017); Witek et al. (2019), and the leading corrections to the gravitational waveforms in head-on collisions Okounkova et al. (2019), but these simulations require much more development before their results can be used in gravitational wave data analysis. There are also concerns about the mathematical viability of the theories considered when they are not treated in the EFT framework.

Thus, we cannot test GR by direct comparison with other specific theories. Instead, we can

(i) check the consistency of the GR predictions with the data and

(ii) introduce ad hoc modifications in GR waveforms to determine the degree to which the values of the deviation parameters agree with GR.

These methods are agnostic to any particular choice of alternative theory. For the most part, our results should therefore be interpreted as observational constraints on possible GW phenomenologies, independent of the overall suitability or well-posedness of any specific alternative to GR. These limits are useful in providing a quantitative indication of the degree to which the data is described by GR; they may also be interpreted more specifically in the context of any given alternative to produce constraints, if applicable.

In particular, with regard to the consistency of the GR predictions (i), we

(a) look for residual power after subtracting the best-fitting GR waveform from the data, and

(b) evaluate the consistency of the high and low frequency components of the observed signal.

With regard to deviations from GR (ii), we separately introduce parametrizations for

(a) the emitted waveform, and

(b) its propagation.

The former could be viewed as representing possible GR modifications in the strong-field region close to the binary, while the latter would correspond to weak-field modifications away from the source. Although we consider these independently, modifications to GW propagation would most likely be accompanied by modifications to GW generation in any given extension of GR. We have also checked that none of the events discussed here provide stronger constraints on models with purely vector and purely scalar GW polarizations than those previously published in Abbott et al. (2017, 2019). Our analyses do not reveal any inconsistency of the data with the predictions of GR. These results supersede all our previous testing GR results on the binary black hole signals found in O1 and O2 Abbott et al. (2016, 2016, 2017, 2017). In particular, the previously published residuals and propagation test results were affected by a slight normalization issue.

Limits on deviations from GR for individual events are dominated by statistical errors due to detector noise. These errors can be reduced by appropriately combining results from multiple events. Sources of systematic errors, on the other hand, include uncertainties in the detector calibration and power spectral density (PSD) estimation and errors in the modeling of waveforms in GR. Detector calibration uncertainties are modeled as corrections to the measured detector response function and are marginalized over. Studies on the effect of PSD uncertainties are currently ongoing. A full characterization of the systematic errors due to the GR waveform models that we employ is beyond the scope of this study; some investigations can be found in Bohé et al. (2017); Khan et al. (2016); Blackman et al. (2017); Khan et al. (2019); Abbott et al. (2017).

This paper is organized as follows. Section II provides an overview of the data sets employed here, while Sec. III details which GW events are used to produce the individual and combined results presented in this paper. In Sec. IV we explain the gravitational waveforms and data analysis formalisms which our tests of GR are based on, before we present the results in the following sections. Section V contains two signal consistency tests: the residuals test in V.1 and the inspiral-merger-ringdown consistency test in V.2. Results from parameterized tests are given in Sec. VI for GW generation, and in Sec. VII for GW propagation. We briefly discuss the study of GW polarizations in Sec. VIII. Finally, we conclude in Sec. IX. We give results for individual events and some checks on waveform systematics in the Appendix.

The results of each test and associated data products can be found in Ref. LIGO Scientific Collaboration and Virgo Collaboration (2019). The GW strain data for all the events considered are available at the Gravitational Wave Open Science Center LIGO Scientific Collaboration, Virgo Collaboration (2018).

II Data, calibration and cleaning

The first observing run of Advanced LIGO (O1) lasted from September 12th, 2015 to January 19th, 2016. The second observing run (O2) lasted from November 30th, 2016 to August 25th, 2017, with the Advanced Virgo observatory joining on August 1st, 2017. This paper includes all GW events originating from the coalescence of two black holes found in these two data sets and published in Abbott et al. (2016, 2019).

The GW detector’s response to changes in the differential arm length (the interferometer’s degree of freedom most sensitive to GWs) must be calibrated using independent, accurate, absolute references. The LIGO detectors use photon recoil (radiation pressure) from auxiliary laser systems to induce mirror motions that change the arm cavity lengths, allowing a direct measure of the detector response  Karki et al. (2016); Cahillane et al. (2017); Viets et al. (2018). Calibration of Virgo relies on measurements of Michelson interference fringes as the main optics swing freely, using the primary laser wavelength as a fiducial length. Subsequent measurements propagate the calibration to arrive at the final detector response Estevez et al. (2018); Acernese et al. (2018). These complex-valued, frequency-dependent measurements of the LIGO and Virgo detectors’ response yield the uncertainty in their respective estimated amplitude and phase of the GW strain output. The amplitude and phase correction factors are modeled as cubic splines and marginalized over in the estimation of astrophysical source parameters  Farr et al. (2015); Veitch et al. (2015); Abbott et al. (2019, 2019). Additionally, the uncertainty in the time stamping of Virgo data (much larger than the LIGO timing uncertainty, which is included in the phase correction factor) is also accounted for in the analysis.

Post-processing techniques to subtract noise contributions and frequency lines from the data around gravitational-wave events were developed in O2 and introduced in Abbott et al. (2017, 2017, 2017), for the astrophysical parameter estimation of GW170608, GW170814, and GW170817. This noise subtraction was achieved using optimal Wiener filters to calculate coupling transfer functions from auxiliary sensors Driggers et al. (2019). A new, optimized parallelizable method in the frequency domain Davis et al. (2019) allows large scale noise subtraction on LIGO data. All of the O2 analyses presented in this manuscript use the noise-subtracted data set with the latest calibration available. The O1 data set is the same used in previous publications, as the effect of noise subtraction is expected to be negligible. Reanalysis of the O1 events is motivated by improvements in the parameter estimation pipeline, an improved frequency-dependent calibration, and the availability of new waveform models.

III Events and significance

We present results for all confident detections of binary black hole events in GTWC-1 LIGO Scientific Collaboration and Virgo Collaboration (2018), i.e., all such events detected during O1 and O2 with a FAR lower than one per year, as published in Abbott et al. (2019). The central columns of Table 1 list the FARs of each event as evaluated by the three search pipelines used in Abbott et al. (2019). Two of these pipelines (PyCBC and GstLAL) rely on waveform templates computed from binary black hole coalescences in GR. Making use of a measure of significance that assumes the validity of GR could potentially lead to biases in the selection of events to be tested, systematically disfavoring signals in which a GR violation would be most evident. Therefore, it is important to consider the possibilities that

(1) there were GW signals with such large deviations from GR that they were missed entirely by the modeled searches, and

(2) there were events that were picked up by the modeled searches but classified as marginal (and thus excluded from our analysis) because of their significant deviations from GR.

These worries can largely be dispelled by considering the third GW search pipeline, the coherent WaveBurst (cWB) weakly modeled search presented in Abbott et al. (2019). This cWB search Klimenko et al. (2008, 2016); Abbott et al. (2016) was tuned to detect chirping signals—like those that would be expected from compact binary coalescences—but was not tuned to any specific GR predictions.222Chirping signals from compact binary coalescences are a feature of many theories of gravity. All that is required is that the orbital frequency increases as the binary radiates energy and angular momentum in GWs. cWB is most sensitive to short signals from high-mass binary black holes. It is still able to detect signals from lower-mass binaries (e.g., GW151226), though with reduced significance compared to the modeled searches. Thus, a signal from a low-mass binary, or a marginal event, with a significant departure from the GR predictions (hence not detected by the GR modeled searches) would not necessarily be detected by the cWB search with a $\mathrm{FAR}<(1\,\mathrm{yr})^{-1}$. However, if there is a population of such signals, they will not all be weak and/or from low-mass binaries. Thus, one would expect some of the signals in the population to be detected by cWB, even if they evade detection by the modeled searches.

All signals detected by the cWB search with $\mathrm{FAR}<(1\,\mathrm{yr})^{-1}$ were also found by at least one modeled search with $\mathrm{FAR}<(1\,\mathrm{yr})^{-1}$. Given the above considerations, this is evidence that our analysis does not exclude chirping GW signals that were missed in the modeled searches because of drastic departures from GR. Similarly, this is also evidence against the possibility of marginal events representing a population of GR-deviating signals, as none of them show high significance [$\mathrm{FAR}<(1\,\mathrm{yr})^{-1}$] in the cWB search only. Thus, we believe that we have not biased our analysis by considering only the ten events with $\mathrm{FAR}<(1\,\mathrm{yr})^{-1}$, as published in Abbott et al. (2019).

We consider each of the GW events individually, carrying out different analyses on a case-by-case basis. Some of the tests presented here, such as the inspiral-merger-ringdown (IMR) consistency test in Sec. V.2 and the parameterized tests in Sec. VI, distinguish between the inspiral and the post-inspiral regimes of the signal. The separation between these two regimes is performed in the frequency-domain, choosing a particular cutoff frequency determined by the parameters of the event. Larger-mass systems merge at lower frequencies, presenting a short inspiral signal in band; lower mass systems have longer observable inspiral signals, but the detector’s sensitivity decreases at higher frequencies and hence the post-inspiral signal becomes less informative. Therefore, depending on the total mass of the system, a particular signal might not provide enough information within the sensitive frequency band of the GW detectors for all analyses.

As a proxy for the amount of information that can be extracted from each part of the signal, we calculate the signal-to-noise ratio (SNR) of the inspiral and the post-inspiral parts of the signals separately. We only apply inspiral (post-inspiral) tests if the inspiral (post-inspiral) SNR is greater than 6. Each test uses a different inspiral-cutoff frequency, and hence they assign different SNRs to the two regimes (details provided in the relevant section for each test). In Table 1 we indicate which analyses have been performed on which event, based on this frequency and the corresponding SNR.333 While we perform these tests on all events with SNR $>6$ in the appropriate regime, in a few cases the results appear uninformative and the posterior distribution extends across the entire prior considered. Since the results are prior dependent, upper limits should not be set from these individual analyses. See Sec. A.3 of the Appendix for details.

In addition to the individual analysis of each event, we derive combined constraints on departures from GR using multiple signals simultaneously. Constraints from individual events are largely dominated by statistical uncertainties due to detector noise. Combining events together can reduce such statistical errors on parameters that take consistent values across all events. However, it is impossible to make joint probabilistic statements from multiple events without prior assumptions about the nature of each observation and how it relates to others in the set. This means that, although there are well-defined statistical procedures for producing joint results, there is no unique way of doing so.

In light of this, we adopt what we take to be the most straightforward strategy, although future studies may follow different criteria. First, in combining events we assume that deviations from GR are manifested equally across events, independent of source properties. This is justified for studies of modified GW propagation, since those effects should not depend on the source.444Propagation effects do depend critically on source distance. However, this dependence is factored out explicitly, in a way that allows for combining events as we do here (see Sec. VII). For other analyses, it is quite a strong assumption to take all deviations from GR to be independent of source properties. Such combined tests should not be expected to necessarily reveal generic source-dependent deviations, although they might if the departures from GR are large enough (see, e.g., Ghosh et al. (2018)). Future work may circumvent this issue by combining marginalized likelihood ratios (Bayes factors), instead of posterior probability distributions Agathos et al. (2014). More general ways of combining results are discussed and implemented in Refs. Zimmerman et al. (2019); Isi et al. (2019).

Second, we choose to produce combined constraints only from events that were found in both modeled searches (PyCBC Nitz et al. (2018); Dal Canton et al. (2014); Usman et al. (2016) and GstLAL Sachdev et al. (2019); Messick et al. (2017)) with a FAR of at most one per one-thousand years. This ensures that there is a very small probability of inclusion of a non-astrophysical event. The events used for the combined results are indicated with bold names in Table 1. The events thus excluded from the combined analysis have low SNR and would therefore contribute only marginally to tightened constraints. Excluding marginal events from our analyses amounts to assigning a null a priori probability to the possibility that those data contain any information about the tests in question. This is, in a sense, the most conservative choice.

In summary, we enforce two significance thresholds: $\mathrm{FAR}<(1\,\mathrm{yr})^{-1}$, for single-event analyses, and $\mathrm{FAR}<(1000\,\mathrm{yr})^{-1}$, for combined results. This two-tiered setup allows us to produce conservative joint results by including only the most significant events, while also providing information about a broader (less significant) set of triggers. This is intended to enable the interested reader to combine individual results with less stringent criteria and under different statistical assumptions, according to their specific needs and tolerance for false positives. In the future, we may adapt our thresholds depending on the rate of detections.

IV Parameter inference

The starting point for all the analyses presented here are waveform models that describe the GWs emitted by coalescing black-hole binaries. The GW signature depends on the intrinsic parameters describing the binary as well as the extrinsic parameters specifying the location and orientation of the source with respect to the detector network. The intrinsic parameters for circularized black-hole binaries in GR are the two masses $m_{i}$ of the black holes and the two spin vectors $\vec{S}_{i}$ defining the rotation of each black hole, where $i\in\{1,2\}$ labels the two black holes. We assume that the binary has negligible orbital eccentricity, as is expected to be the case when the binary enters the band of ground-based detectors Peters and Mathews (1963); Peters (1964) (except in some more extreme formation scenarios,888These scenarios could occur often enough, compared to the expected rate of detections, that the inclusion of eccentricity in waveform models is a necessity for tests of GR in future observing runs; see, e.g., Hinderer and Babak (2017); Cao and Han (2017); Hinder et al. (2018); Huerta et al. (2018); Klein et al. (2018); Moore et al. (2018); Moore and Yunes (2019); Tiwari et al. (2019) for recent work on developing such waveform models. e.g., Samsing (2018); Rodriguez et al. (2018); Zevin et al. (2019); Rodriguez et al. (2018)). The extrinsic parameters comprise four parameters that specify the space-time location of the binary black-hole, namely the sky location (right ascension and declination), the luminosity distance, and the time of coalescence. In addition, there are three extrinsic parameters that determine the orientation of the binary with respect to Earth, namely the inclination angle of the orbit with respect to the observer, the polarization angle, and the orbital phase at coalescence.

We employ two waveform families that model binary black holes in GR: the effective-one-body based SEOBNRv4 Bohé et al. (2017) waveform family that assumes non-precessing spins for the black holes (we use the frequency domain reduced order model SEOBNRv4_ROM for reasons of computational efficiency), and the phenomenological waveform family IMRPhenomPv2 Husa et al. (2016); Khan et al. (2016); Hannam et al. (2014) that models the effects of precessing spins using two effective parameters by twisting up the underlying aligned-spin model. We use IMRPhenomPv2 to obtain all the main results given in this paper, and use SEOBNRv4 to check the robustness of these results, whenever possible. When we use IMRPhenomPv2, we impose a prior $m_{1}/m_{2}\leq 18$ on the mass ratio, as the waveform family is not calibrated against numerical relativity simulations for $m_{1}/m_{2}>18$. We do not impose a similar prior when using SEOBNRv4, since it includes information about the extreme mass ratio limit. Neither of these waveform models includes the full spin dynamics (which requires 6 spin parameters). Fully precessing waveform models have been recently developed Taracchini et al. (2014); Babak et al. (2017); Blackman et al. (2017); Khan et al. (2019); Varma et al. (2019) and will be used in future applications of these tests.

The waveform models used in this paper do not include the effects of subdominant (non-quadrupole) modes, which are expected to be small for comparable-mass binaries Kamaretsos et al. (2012); London et al. (2014). The first generation of binary black hole waveform models including spin and higher order modes has recently been developed Blackman et al. (2017); London et al. (2018); Cotesta et al. (2018); Varma et al. (2019, 2019). Preliminary results in Abbott et al. (2019), using NR simulations supplemented by NR surrogate waveforms, indicate that the higher mode content of the GW signals detected by Advanced LIGO and Virgo is weak enough that models without the effect of subdominant modes do not introduce substantial biases in the intrinsic parameters of the binary. For unequal-mass binaries, the effect of the non-quadrupole modes is more pronounced Varma and Ajith (2017), particularly when the binary’s orientation is close to edge-on. In these cases, the presence of non-quadrupole modes can show up as a deviation from GR when using waveforms that only include the quadrupole modes, as was shown in Pang et al. (2018). Applications of tests of GR with the new waveform models that include non-quadrupole modes will be carried out in the future.

We believe that our simplifying assumptions on the waveform models (zero eccentricity, simplified treatment of spins, and neglect of subdominant modes) are justified by astrophysical considerations and previous studies. Indeed, as we show in the remainder of the paper, the observed signals are consistent with the waveform models. Of course, had our analyses resulted in evidence for violations of GR, we would have had to revisit these simplifications very carefully.

The tests described in this paper are performed within the framework of Bayesian inference, by means of the LALInference code Veitch et al. (2015) in the LIGO Scientific Collaboration Algorithm Library Suite (LALSuite) LIGO Scientific Collaboration and Virgo Collaboration (2018). We estimate the PSD using the BayesWave code Cornish and Littenberg (2015); Littenberg and Cornish (2015), as described in Appendix B of Abbott et al. (2019). Except for the residuals test described in Sec. V.1, we use the waveform models described in this section to estimate from the data the posterior distributions of the parameters of the binary. These include not only the intrinsic and extrinsic parameters mentioned above, but also other parameters that describe possible departures from the GR predictions. Specifically, for the parameterized tests in Secs. VI and VII, we modify the phase $\Phi(f)$ of the frequency-domain waveform

[TABLE]

For the GR parameters, we use the same prior distributions as the main parameter estimation analysis described in Abbott et al. (2019), though for a number of the tests we need to extend the ranges of these priors to account for correlations with the non-GR parameters, or for the fact that only a portion of the signal is analyzed (as in Sec. V.2). We also use the same low-frequency cutoffs for the likelihood integral as in Abbott et al. (2019), i.e., $20$ Hz for all events except for GW170608, where $30$ Hz is used for LIGO Hanford, as discussed in Abbott et al. (2017), and GW170818, where $16$ Hz is used for all three detectors. For the model agnostic residual test described in Sec. V.1, we use the BayesWave code Cornish and Littenberg (2015) which describes the GW signals in terms of a number of Morlet-Gabor wavelets.

V Consistency tests

V.1 Residuals test

One way to evaluate the ability of GR to describe GW signals is to subtract the best-fit template from the data and make sure the residuals have the statistical properties expected of instrumental noise. This largely model-independent test is sensitive to a wide range of possible disagreements between the data and our waveform models, including those caused by deviations from GR and by modeling systematics. This analysis can look for GR violations without relying on specific parametrizations of the deviations, making it a versatile tool. Results from a similar study on our first detection were already presented in Abbott et al. (2016).

In order to establish whether the residuals agree with noise (Gaussian or otherwise), we proceed as follows. For each event in our set, we use LALInference and the IMRPhenomPv2 waveform family to obtain an estimate of the best-fit (i.e., maximum likelihood) binary black-hole waveform based on GR. This waveform incorporates factors that account for uncertainty in the detector calibration, as described in Sec. II. This best-fit waveform is then subtracted from the data to obtain residuals for a 1 second window centered on the trigger time reported in Abbott et al. (2019).999The analysis is sensitive only to residual power in that 1 s window due to technicalities related to how BayesWave handles its sine-Gaussian basis elements Cornish and Littenberg (2015); Littenberg and Cornish (2015). If the GR-based model provides a good description of the signal, we expect the resulting residuals at each detector to lack any significant coherent SNR beyond what is expected from noise fluctuations. We compute the coherent SNR using BayesWave, which models the multi-detector residuals as a superposition of incoherent Gaussian noise and an elliptically-polarized coherent signal. The residual network SNR reported by BayesWave is the SNR that would correspond to such a coherent signal.

In particular, for each event, BayesWave produces a distribution of possible residual signals consistent with the data, together with corresponding a posteriori probabilities. This is trivially translated into a probability distribution over the coherent residual SNR. We summarize each of these distributions by computing the corresponding 90%-credible upper limit ($\mathrm{SNR}_{90}$). This produces one number per event that represents an upper bound on the coherent power that could be present in its residuals.

We may translate the $\mathrm{SNR}_{90}$ into a measure of how well the best-fit templates describe the signals in our data. We do this through the fitting factor Apostolatos (1995), $\mathrm{FF}\coloneqq\mathrm{SNR}_{\mathrm{GR}}\,/(\mathrm{SNR}_{\mathrm{res}}^{2}+\mathrm{SNR}_{\mathrm{GR}}^{2})^{1/2}$, where $\mathrm{SNR}_{\mathrm{res}}$ is the coherent residual SNR and $\mathrm{SNR}_{\mathrm{GR}}$ is the network SNR of the best-fit template (see Table 1 for network SNRs). By setting $\mathrm{SNR}_{\mathrm{res}}=\mathrm{SNR}_{90}$, we produce a 90%-credible lower limit on the fitting factor ($\mathrm{FF}_{90}$). Because FF is itself a lower limit on the overlap between the true and best-fit templates, so is $\mathrm{FF}_{90}$. As in Abbott et al. (2016), we may then assert that the disagreement between the true waveform and our GR-based template is at most $(1-\mathrm{FF}_{90})\times 100\%$. This is interesting as a measure of the sensitivity of our test, but does not tell us about the consistency of the residuals with instrumental noise.

To assess whether the obtained residual $\mathrm{SNR}_{90}$ values are consistent with detector noise, we run an identical BayesWave analysis on 200 different sets of noise-only detector data near each event. This allows us to estimate the $p$-value for the null hypothesis that the residuals are consistent with noise. The $p$-value gives the probability of noise producing coherent power with $\mathrm{SNR}_{90}^{\mathrm{n}}$ greater than or equal to the residual value $\mathrm{SNR}_{90}$, i.e., $p\coloneqq P(\mathrm{SNR}_{90}^{\mathrm{n}}\geq\mathrm{SNR}_{90}\mid\mathrm{noise})$. In that sense, a smaller $p$-value indicates a smaller chance that the residual power arose from instrumental noise only. For each event, our estimate of $p$ is produced from the fraction of noise instantiations that yielded $\mathrm{SNR}_{90}^{\mathrm{n}}\geq\mathrm{SNR}_{90}$ (that is, from the empirical survival function).101010Computing $p$-values would not be necessary if the noise was perfectly Gaussian, in which case we could predict the noise-only distribution of $\mathrm{SNR}_{90}^{\mathrm{n}}$ from first principles.

Our results are summarized in Table 2. For each event, we present the values of the residual $\mathrm{SNR}_{90}$, the lower limit on the fitting factor ($\mathrm{FF}_{90}$), and the $\mathrm{SNR}_{90}$ $p$-value. The background distributions that resulted in those $p$-values are shown in Fig. 1. In Fig. 1 we represent these distributions through the empirical estimate of their survival functions, i.e., $p(\mathrm{SNR}_{90})=1-\mathrm{CDF}(\mathrm{SNR}_{90})$, with “CDF” the cumulative distribution function. Fig. 1 also displays the actual value of $\mathrm{SNR}_{90}$ measured from the residuals of each event (dotted vertical line). In each case, the height of the curve evaluated at the $\mathrm{SNR}_{90}$ measured for the corresponding detection yields the $p$-value reported in Table 2 (markers in Fig. 1).

The values of residual $\mathrm{SNR}_{90}$ vary widely among events because they depend on the specific state of the instruments at the time of detection: segments of data with elevated noise levels are more likely to result in spurious coherent residual power, even if the signal agrees with GR. In particular, the background distributions for events seen by three detectors are qualitatively different from those seen by only two. This is both due to

(i) the fact that BayesWave is configured to expect the SNR to increase with the number of detectors and

(ii) the fact that Virgo data present a higher rate of non-Gaussianities than LIGO.

We have confirmed both these factors play a role by studying the background $\mathrm{SNR}_{90}$ distributions for real data from each possible pair of detectors, as well as distributions over simulated Gaussian noise. Specifically, removing Virgo from the analysis results in a reduction in the coherent residual power for GW170729 ($\mathrm{SNR}_{90}^{\rm HL}=6.5$), GW170809 ($\mathrm{SNR}_{90}^{\rm HL}=6.3$), GW170814 ($\mathrm{SNR}_{90}^{\rm HL}=6.0$), and GW170818 ($\mathrm{SNR}_{90}^{\rm HL}=7.2$).

The event-by-event variation of $\mathrm{SNR}_{90}$ is also reflected in the values of $\mathrm{FF}_{90}$. GW150914 provides the strongest result with $\mathrm{FF}_{90}={0.97}$, which corresponds to an upper limit of 3% on the magnitude of potential deviations from our GR-based template,111111This value is better than the one quoted in Abbott et al. (2016) by 1 percentage point. The small difference is explained by several factors, including that paper’s use of the maximum a posteriori waveform (instead of maximum likelihood) and 95% (instead of 90%) credible intervals, as well as improvements in data calibration. in the specific sense defined in Abbott et al. (2016) and discussed above. On the other hand, GW170818 yields the weakest result with $\mathrm{FF}_{90}=0.78$ and a corresponding upper limit on waveform disagreement of $1-\mathrm{FF}_{90}=22\%$. The average $\mathrm{FF}_{90}$ over all events is 0.88.

The set of $p$-values shown in Table 2 is consistent with all coherent residual power being due to instrumental noise. Assuming that this is indeed the case, we expect the $p$-values to be uniformly distributed over $[0,1]$, which explains the variation in Table 2. With only ten events, however, it is difficult to obtain strong quantitative evidence of the uniformity of this distribution. Nevertheless, we follow Fisher’s method Fisher (1948) to compute a meta $p$-value for the null hypothesis that the individual $p$-values in Table 2 are uniformly distributed. We obtain a meta $p$-value of 0.4, implying that there is no evidence for coherent residual power that cannot be explained by noise alone. All in all, this means that there is no statistically significant evidence for deviations from GR.

V.2 Inspiral-merger-ringdown consistency test

The inspiral-merger-ringdown consistency test for binary black holes Ghosh et al. (2016, 2018) checks the consistency of the low-frequency part of the observed signal (roughly corresponding to the inspiral of the black holes) with the high-frequency part (to a good approximation, produced by the post-inspiral stages).121212Note that this is not exactly equal to testing the consistency between the early and late part of the waveform in time domain, because the low-frequency part of the signal could be “contaminated” by power from late times and vice versa. In practice, this effect is negligible with our choice of cutoff frequencies. See Ghosh et al. (2018) for a discussion. The cutoff frequency $f_{\text{c}}$ between the two regimes is chosen as the frequency of the innermost stable circular orbit of a Kerr black hole Bardeen et al. (1972), with mass and dimensionless spin computed by applying NR fits Healy and Lousto (2017); Hofmann et al. (2016); Jiménez-Forteza et al. (2017); Johnson-McDaniel et al. (2016) to the median values of the posterior distributions of the initial masses and spherical coordinate components of the spins. This determination of $f_{\text{c}}$ is performed separately for each event and based on parameter inference of the full signal (see Table 3 for values of $f_{\text{c}}$).131313 The frequency $f_{\text{c}}$ was determined using preliminary parameter inference results, so the values in Table 3 are slightly different than those that would be obtained using the posterior samples in GWTC-1 LIGO Scientific Collaboration and Virgo Collaboration (2018). However, the test is robust against small changes in the cutoff frequency Ghosh et al. (2018). The binary’s parameters are then estimated independently from the low (high) frequency parts of the data by restricting the noise-weighted integral in the likelihood calculation to frequencies below (above) this frequency cutoff $f_{\text{c}}$. For each of these independent estimates of the source parameters, we make use of fits to numerical-relativity simulations given in Healy and Lousto (2017); Hofmann et al. (2016); Jiménez-Forteza et al. (2017) to infer the mass $M_{\mathrm{f}}$ and dimensionless spin magnitude $a_{\mathrm{f}}=c|\vec{S}_{\mathrm{f}}|/(GM_{\mathrm{f}}^{2})$ of the remnant black hole.141414As in Abbott et al. (2017), we average the $M_{\mathrm{f}},a_{\mathrm{f}}$ posteriors obtained by different fits Healy and Lousto (2017); Hofmann et al. (2016); Jiménez-Forteza et al. (2017) after augmenting the fitting formulae for aligned-spin binaries by adding the contribution from in-plane spins Johnson-McDaniel et al. (2016). However, unlike in Abbott et al. (2017); Johnson-McDaniel et al. (2016), we do not evolve the spins before applying the fits, due to technical reasons. If the data are consistent with GR, these two independent estimates have to be consistent with each other Ghosh et al. (2016, 2018). Because this consistency test ultimately compares between the inspiral and the post-inspiral results, posteriors of both parts must be informative. In the case of low-mass binaries, the SNR in the part $f>f_{\text{c}}$ is insufficient to perform this test, so that we only analyze seven events as indicated in Tables 1 and 3.

In order to quantify the consistency of the two different estimates of the final black hole’s mass and spin we define two dimensionless quantities that quantify the fractional difference between them: $\Delta M_{\mathrm{f}}/\bar{M}_{\mathrm{f}}\coloneqq 2\,(M_{\mathrm{f}}^{\text{insp}}-M_{\mathrm{f}}^{\text{post-insp}})/(M_{\mathrm{f}}^{\text{insp}}+M_{\mathrm{f}}^{\text{post-insp}})$ and $\Delta a_{\mathrm{f}}/\bar{a}_{\mathrm{f}}\coloneqq 2\,(a_{\mathrm{f}}^{\text{insp}}-a_{\mathrm{f}}^{\text{post-insp}})/(a_{\mathrm{f}}^{\text{insp}}+a_{\mathrm{f}}^{\text{post-insp}})$, where the superscripts indicate the estimates of the mass and spin from the inspiral and post-inspiral parts of the signal.151515For black hole binaries with comparable masses and moderate spins, as we consider here, the remnant black hole is expected to have $a_{\mathrm{f}}\gtrsim 0.5$; see, e.g., Healy and Lousto (2017); Hofmann et al. (2016); Jiménez-Forteza et al. (2017) for fitting formulae derived from numerical simulations, or Table 1 for values of the remnant’s spins obtained from GW events. Hence, $\Delta a_{\mathrm{f}}/\bar{a}_{\mathrm{f}}$ is expected to yield finite values. The posteriors of these dimensionless parameters, estimated from different events, are shown in Fig. 2. For all events, the posteriors are consistent with the GR value ($\Delta M_{\mathrm{f}}/\bar{M}_{\mathrm{f}}=0,\Delta a_{\mathrm{f}}/\bar{a}_{\mathrm{f}}=0)$. The fraction of the posterior enclosed by the isoprobability contour that passes through the GR value (i.e., the GR quantile) for each event is shown in Table 3. Figure 2 also shows the posteriors obtained by combining all the events that pass the stronger significance threshold $\mathrm{FAR}<(1000\,\mathrm{yr})^{-1}$, as outlined in Sec. III (see the same section for a discussion of caveats).

The parameter estimation is performed employing uniform priors in component masses and spin magnitudes and isotropic priors in spin directions Abbott et al. (2019). This will introduce a non-flat prior in the deviation parameters $\Delta M_{\mathrm{f}}/\bar{M}_{\mathrm{f}}$ and $\Delta a_{\mathrm{f}}/\bar{a}_{\mathrm{f}}$, which is shown as a thin, dashed contour in Fig. 2. Posteriors are estimated employing the precessing spin phenomenological waveform family IMRPhenomPv2. To assess the systematic errors due to imperfect waveform modeling, we also estimate the posteriors using the effective-one-body based waveform family SEOBNRv4 that models binary black holes with non-precessing spins. There is no qualitative difference between the results derived using the two different waveforms families (see Sec. A.2 of the Appendix).

We see additional peaks in the posteriors estimated from GW170814 and GW170823. Detailed follow-up investigations did not show any evidence of the presence of a coherent signal in multiple detectors that differs from the GR prediction. The second peak in GW170814 is introduced by the posterior of $M_{\mathrm{f}}^{\text{post-insp}}$, while the extra peak in GW170823 is introduced by the posterior of $M_{\mathrm{f}}^{\text{insp}}$. Injection studies in real data around the time of these events, using simulated GR waveforms with parameters consistent with GW170814 and GW170823, suggest that such secondary peaks occur for $\sim 10\%$ of injections. Features in the posteriors of GW170814 and GW170823 are thus consistent with expected noise fluctuations.

VI Parameterized tests of gravitational wave generation

A deviation from GR could manifest itself as a modification of the dynamics of two orbiting compact objects, and in particular, the evolution of the orbital (and hence, GW) phase. In an analytical waveform model like IMRPhenomPv2, the details of the GW phase evolution are controlled by coefficients that are either analytically calculated or determined by fits to numerical-relativity (NR) simulations, always under the assumption that GR is the underlying theory. In this section we investigate deviations from the GR binary dynamics by introducing shifts in each of the individual GW phase coefficients of IMRPhenomPv2. Such shifts correspond to deviations in the waveforms from the predictions of GR. We then treat these shifts as additional unconstrained GR-violating parameters, which we measure in addition to the standard parameters describing the binary.

The early inspiral of compact binaries is well modeled by the post-Newtonian (PN) approximation Blanchet et al. (1995, 2004, 2005); Blanchet (2014) to GR, which is based on the expansion of the orbital quantities in terms of a small velocity parameter $v/c$. For a given set of intrinsic parameters, coefficients for the different orders in $v/c$ in the PN series are uniquely determined. A consistency test of GR using measurements of the inspiral PN phase coefficients was first proposed in Arun et al. (2006, 2006); Mishra et al. (2010), while a generalized parametrization was motivated in Yunes and Pretorius (2009). Bayesian implementations based on such parametrized methods were presented and tested in Li et al. (2012, 2012); Agathos et al. (2014); Cornish et al. (2011) and were also extended to the post-inspiral part of the gravitational-wave signal Sampson et al. (2014); Meidam et al. (2018). These ideas were applied to the first GW observation, GW150914 Abbott et al. (2016), yielding the first bounds on higher-order PN coefficients Abbott et al. (2016). Since then, the constraints have been revised with the binary black hole events that followed, GW151226 in O1 Abbott et al. (2016) and GW170104 in O2 Abbott et al. (2017). More recently, the first such constraints from a binary neutron star merger were placed with the detection of GW170817 Abbott et al. (2019). Bounds on parametrized violations of GR from GW detections have been mapped, to leading order, to constraints on specific alternative theories of gravity (see, e.g., Yunes et al. (2016)). In this paper, we present individual constraints on parametrized deviations from GR for each of the GW sources in O1 and O2 listed in Table 1, as well as the tightest combined constraints obtained to date by combining information from all the significant binary black holes events observed so far, as described in Sec. III.

The frequency-domain GW phase evolution $\Phi(f)$ in the early-inspiral stage of IMRPhenomPv2 is described by a PN expansion, augmented with higher-order phenomenological coefficients. The PN phase evolution is analytically expressed in closed form by employing the stationary phase approximation. The late-inspiral and post-inspiral (intermediate and merger-ringdown) stages are described by phenomenological analytical expressions. The transition frequency161616This frequency is different than the cutoff frequency used in the inspiral-merger-ringdown consistency test, as was briefly mentioned in Sec. III. from inspiral to intermediate regime is given by the condition $GMf/c^{3}=0.018$, with $M$ the total mass of the binary in the detector frame, since this is the lowest frequency above which this model was calibrated with NR data Khan et al. (2016). Let us use $p_{i}$ to collectively denote all of the inspiral and post-inspiral parameters $\varphi_{i}$, $\alpha_{i}$, $\beta_{i}$, that will be introduced below. Deviations from GR in all stages are expressed by means of relative shifts $\delta\hat{p}_{i}$ in the corresponding waveform coefficients: $p_{i}\rightarrow(1+\delta\hat{p}_{i})\ p_{i}$, which are used as additional free parameters in our extended waveform models.

We denote the testing parameters corresponding to PN phase coefficients by $\delta\hat{\varphi}_{i}$, where $i$ indicates the power of $v/c$ beyond leading (Newtonian or [math]PN) order in $\Phi(f)$. The frequency dependence of the corresponding phase term is $f^{(i-5)/3}$. In the parametrized model, $i$ varies from 0 to 7, including the terms with logarithmic dependence at $2.5$PN and $3$PN. The non-logarithmic term at $2.5$PN (i.e., $i=5$) cannot be constrained, because of its degeneracy with a constant reference phase (e.g., the phase at coalescence). These coefficients were introduced in their current form in Eq. (19) of Li et al. (2012). In addition, we also test for $i=-2$, representing an effective $-1$PN term, which is motivated below. The full set of inspiral parameters are thus

[TABLE]

Since the $-1$PN term and the $0.5$PN term are absent in the GR phasing, we parametrize $\delta\hat{\varphi}_{-2}$ and $\delta\hat{\varphi}_{1}$ as absolute deviations, with a pre-factor equal to the [math]PN coefficient.

The $-1$PN term of $\delta\hat{\varphi}_{2}$ can be interpreted as arising from the emission of dipolar radiation. For binary black holes, this could occur in, e.g., alternative theories of gravity where an additional scalar charge is sourced by terms related to curvature Yagi et al. (2016); Yagi and Stein (2016). At leading order, this introduces a deviation in the $-1$PN coefficient of the waveform Barausse et al. (2016); Arun (2012). This effectively introduces a term in the inspiral GW phase, varying with frequency as $f^{-7/3}$, while the gravitational flux is modified as $\mathcal{F_{\mathrm{GR}}}\rightarrow\mathcal{F_{\mathrm{GR}}}(1+Bc^{2}/v^{2})$. The first bound on $\delta\hat{\varphi}_{-2}$ was published in Abbott et al. (2019). The higher-order terms in the above expansion also lead to a modification in the higher-order PN coefficients. Unlike the case of GW170817 (which we study separately in Abbott et al. (2019)), where the higher-order terms in the expansion of the flux are negligible, the contribution of higher-order terms can be significant in the binary black-hole signals that we study here. This prohibits an exact interpretation of the $-1$PN term as the strength of dipolar radiation. Hence, this analysis only serves as a test of the presence of an effective $-1$PN term in the inspiral phasing, which is absent in GR.

To measure the above GR violations in the post-Newtonian inspiral, we employ two waveform models: (i) the analytical frequency-domain model IMRPhenomPv2 which also provided the natural parametrization for our tests and (ii) SEOBNRv4, which we use in the form of SEOBNRv4_ROM, a frequency-domain, reduced-order-model of the SEOBNRv4 model. The inspiral part of SEOBNRv4 is based on a numerical evolution of the aligned-spin effective-one-body dynamics of the binary and its post-inspiral model is phenomenological. The entire SEOBNRv4 model is calibrated against NR simulations. Despite its non-analytical nature, SEOBNRv4_ROM can also be used to test the parametrized modifications of the early inspiral defined above. Using the method presented in Abbott et al. (2019), we add deviations to the waveform phase corresponding to a given $\delta\hat{\varphi}_{i}$ at low frequencies and then taper the corrections to zero at a frequency consistent with the transition frequency between early-inspiral and intermediate phases used by IMRPhenomPv2. The same procedure cannot be applied to the later stages of the waveform, thus the analysis performed with SEOBNRv4 is restricted to the post-Newtonian inspiral, cf. Fig. 3.

The analytical descriptions of the intermediate and merger-ringdown stages in the IMRPhenomPv2 model allow for a straightforward way of parametrizing deviations from GR, denoted by $\{\delta\hat{\beta}_{2},\delta\hat{\beta}_{3}\}$ and $\{\delta\hat{\alpha}_{2},\delta\hat{\alpha}_{3},\delta\hat{\alpha}_{4}\}$ respectively, following Meidam et al. (2018). Here the parameters $\delta\hat{\beta}_{i}$ correspond to deviations from the NR-calibrated phenomenological coefficients $\beta_{i}$ of the intermediate stage, while the parameters $\delta\hat{\alpha}_{i}$ refer to modifications of the merger-ringdown coefficients $\alpha_{i}$ obtained from a combination of phenomenological fits and analytical black-hole perturbation theory calculations Khan et al. (2016).

Using LALInference, we calculate posterior distributions of the parameters characterizing the waveform (including those that describe the binary in GR). Our parametrization recovers GR at $\delta\hat{p}_{i}=0$, so consistency with GR is verified if the posteriors of $\delta\hat{p}_{i}$ have support at zero. We perform the analyses by varying one $\delta\hat{p}_{i}$ at a time; as shown in Ref. Sampson et al. (2013), this is fully robust to detecting deviations present in multiple PN-orders. In addition, allowing for a larger parameter space by varying multiple coefficients simultaneously would not improve our efficiency in identifying violations of GR, as it would yield less informative posteriors. A specific alternative theory of gravity would likely yield correlated deviations in many parameters, including modifications that we have not considered here. This would be the target of an exact comparison of an alternative theory with GR, which would only be possible if a complete, accurate description of the inspiral-merger-ringdown signal in that theory was available.

We use priors uniform on $\delta\hat{p}_{i}$ and symmetric around zero. Figure 3 shows the combined posteriors of $\delta\hat{p}_{i}$ (marginalized over all other parameters) estimated from the combination of all the events that cross the significance threshold of $\mathrm{FAR<(1000~{}yr)}^{-1}$ in both modeled searches; see Table 1. Events with SNR$<6$ in the inspiral regime (parameters $\delta\hat{\varphi}_{i}$) or in the post-inspiral regime ($\delta\hat{\beta}_{i}$ and $\delta\hat{\alpha}_{i}$ for the intermediate and merger-ringdown parameters respectively) are not included in the results, since the data from those instances failed to provide useful constraints (see Sec. III for more details). This SNR threshold, however, is not equally effective in ensuring informative results for all cases; see Sec. A.3 in the Appendix for a detailed discussion. In all cases considered, the posteriors are consistent with $\delta\hat{p}_{i}=0$ within statistical fluctuations. Bounds on the inspiral coefficients obtained with the two different waveform models are found to be in good agreement with each other. Finally, we note that the event-combining analyses on $\delta\hat{p}_{i}$ assume that these parametrized violations are constant across all events considered. This assumption should not be made when testing a specific theory that predicts violations that depend on the binary’s parameters. Posterior distributions of $\delta\hat{p}_{i}$ for the individual-event analysis, also showing full consistency with GR, are provided in Sec. A.3 of the Appendix.

Figure 4 shows the $90\%$ upper bounds on $|\delta\hat{\varphi}_{i}|$ for all the individual events which cross the SNR threshold (SNR $>$ 6) in the inspiral regime (the most massive of which is GW150914). The bounds from the combined posteriors are also shown; these include the events which exceed both the SNR threshold in the inspiral regime as well as the significance threshold, namely GW150914, GW151226, GW170104, GW170608, and GW170814. The bound from the likely lightest mass binary black hole event GW170608 at 1.5PN is currently the strongest constraint obtained on a positive PN coefficient from a single binary black hole event, as shown in Fig. 4. However, the constraint at this order is about five times worse than that obtained from the binary neutron star event GW170817 alone Abbott et al. (2019). The $-1$PN bound is two orders of magnitude better for GW170817 than the best bound obtained here (from GW170608). The corresponding best $-1$PN bound coming from the double pulsar PSR J0737$-$3039, is a few orders of magnitude tighter still, at $|\delta\hat{\varphi}_{-2}|\lesssim 10^{-7}$ Barausse et al. (2016); Kramer et al. (2006). At [math]PN we find that the bound from GW170608 beats the one from GW170817, but remains weaker than the one from the double pulsar by one order of magnitude Yunes and Hughes (2010); Kramer et al. (2006). For all other PN orders, GW170608 also provides the best bounds, which at high PN orders are of the same order of magnitude as the ones from GW170817. Our results can be compared statistically to those obtained by performing the same tests on simulated GR and non-GR waveforms given in Meidam et al. (2018). The results presented here are consistent with those of GR waveforms injected into realistic detector data. The combined bounds are the tightest obtained so far, improving on the bounds obtained in Abbott et al. (2016) by factors between 1.1 and 1.8.

VII Parameterized tests of gravitational wave propagation

We now place constraints on a phenomenological modification of the GW dispersion relation, i.e., on a possible frequency dependence of the speed of GWs. This modification, introduced in Mirshekari et al. (2012) and first applied to LIGO data in Abbott et al. (2017), is obtained by adding a power-law term in the momentum to the dispersion relation $E^{2}=p^{2}c^{2}$ of GWs in GR, giving

[TABLE]

Here, $c$ is the speed of light, $E$ and $p$ are the energy and momentum of the GWs, and $A_{\alpha}$ and $\alpha$ are phenomenological parameters. We consider $\alpha$ values from [math] to $4$ in steps of $0.5$. However, we exclude $\alpha=2$, where the speed of the GWs is modified in a frequency-independent manner, and therefore gives no observable dephasing.171717For a source with an electromagnetic counterpart, $A_{2}$ can be constrained by comparison with the arrival time of the photons, as was done with GW170817/GRB170817A Abbott et al. (2017). Thus, in all cases except for $\alpha=0$, we are considering a Lorentz-violating dispersion relation. The group velocity associated with this dispersion relation is $v_{g}/c=(dE/dp)/c=1+(\alpha-1)A_{\alpha}E^{\alpha-2}/2+O(A_{\alpha}^{2})$. The associated length scale is $\lambda_{A}\coloneqq hc|A_{\alpha}|^{1/(\alpha-2)}$, where $h$ is Planck’s constant. $\lambda_{A}$ gives the scale of modifications to the Newtonian potential (the Yukawa potential for $\alpha=0$) associated with this dispersion relation.

While Eq. (2) is a purely phenomenological model, it encompasses a variety of more fundamental predictions (at least to leading order) Mirshekari et al. (2012); Yunes et al. (2016). In particular, $A_{0}>0$ corresponds to a massive graviton, i.e., the same dispersion as for a massive particle in vacuo Will (1998), with a graviton mass given by $m_{g}=A_{0}^{1/2}/c^{2}$.181818Thus, the Yukawa screening length is $\lambda_{0}=h/(m_{g}c)$. Furthermore, $\alpha$ values of $2.5$, $3$, and $4$ correspond to the leading predictions of multi-fractal spacetime Calcagni (2010); doubly special relativity Amelino-Camelia (2002); and Hořava-Lifshitz Hořava (2009) and extra dimensional Sefiedgar et al. (2011) theories, respectively. The standard model extension also gives a leading contribution with $\alpha=4$ Kostelecký and Mewes (2016), only considering the non-birefringent terms; our analysis does not allow for birefringence.

In order to obtain a waveform model with which to constrain these propagation effects, we start by assuming that the waveform extracted in the binary’s local wave zone (i.e., near to the binary compared to the distance from the binary to Earth, but far from the binary compared to its own size) is well-described by a waveform in GR.191919This is likely to be a good assumption for $\alpha<2$, where we constrain $\lambda_{A}$ to be much larger than the size of the binary. For $\alpha>2$, where we constrain $\lambda_{A}$ to be much smaller than the size of the binary, one has to posit a screening mechanism in order to be able to assume that the waveform in the binary’s local wave zone is well-described by GR, as well as for this model to evade Solar System constraints. Since we are able to bound these propagation effects to be very small, we can work to linear order in $A_{\alpha}$ when computing the effects of this dispersion on the frequency-domain GW phasing,202020The dimensionless parameter controlling the size of the linear correction is $A_{\alpha}f^{\alpha-2}$, which is $\lesssim 10^{-18}$ at the 90% credible level for the events we consider and frequencies up to $1$ kHz. thus obtaining a correction Mirshekari et al. (2012) that is added to $\Phi(f)$ in Eq. (1):

[TABLE]

Here, $D_{\text{L}}$ is the binary’s luminosity distance, $\mathcal{M}^{\text{det}}$ is the binary’s detector-frame (i.e., redshifted) chirp mass, and $\lambda_{A,\text{eff}}$ is the effective wavelength parameter used in the sampling, defined as

[TABLE]

The parameter $z$ is the binary’s redshift, and $D_{\alpha}$ is a distance parameter given by

[TABLE]

where $H_{0}=67.90\text{ km s}^{-1}\text{ Mpc}^{-1}$ is the Hubble constant, and $\Omega_{\text{m}}=0.3065$ and $\Omega_{\Lambda}=0.6935$ are the matter and dark energy density parameters; these are the TT+lowP+lensing+ext values from Ade et al. (2016).212121We use these values for consistency with the results presented in Abbott et al. (2019). If we instead use the more recent results from Aghanim et al. (2018), specifically the TT,TE,EE+lowE+lensing+BAO values used for comparison in Abbott et al. (2019), then there are very minor changes to the results presented in this section. For instance, the upper bounds in Table 4 change by at most $\sim 0.1\%$.

The dephasing in Eq. (3) is obtained by treating the gravitational wave as a stream of particles (gravitons), which travel at the particle velocity $v_{p}/c=pc/E=1-A_{\alpha}E^{\alpha-2}/2+O(A_{\alpha}^{2})$. There are suggestions to use the particle velocity when considering doubly special relativity, though there are also suggestions to use the group velocity $v_{g}$ in that case (see, e.g., Amelino-Camelia et al. (2006) and references therein for both arguments). However, the group velocity is appropriate for, e.g., multi-fractal spacetime theories (see, e.g., Calcagni (2017)). To convert the bounds presented here to the case where the particles travel at the group velocity, scale the $A_{\alpha}$ bounds for $\alpha\neq 1$ by factors of $1/(1-\alpha)$. The group velocity calculation gives an unobservable constant phase shift for $\alpha=1$.

We consider the cases of positive and negative $A_{\alpha}$ separately, and obtain the results shown in Table 4 and Fig. 5 when applying this analysis to the GW events under consideration. While we sample with a flat prior in $\log\lambda_{A,\text{eff}}$, our bounds are given using priors flat in $A_{\alpha}$ for all results except for the mass of the graviton, where we use a prior flat in the graviton mass. We also show the results from combining together all the signals that satisfy our selection criterion. We are able to combine together the results from different signals with no ambiguity, since the known distance dependence is accounted for in the waveforms.

Figure 6 displays the full $A_{\alpha}$ posteriors obtained by combining all selected events (using IMRPhenomPv2 waveforms). To obtain the full $A_{\alpha}$ posteriors, we combine together the positive and negative $A_{\alpha}$ results for individual events by weighting by their Bayesian evidences; we then combine the posteriors from individual events. We give the analogous plots for the individual events in Sec. A.4 of the Appendix. The combined positive and negative $A_{\alpha}$ posteriors are also used to compute the GR quantiles given in Table 4, which give the probability to have $A_{\alpha}<0$, where $A_{\alpha}=0$ is the GR value. Thus, large or small values of the GR quantile indicate that the distribution is not peaked close to the GR value. For a GR signal, the GR quantile will be distributed uniformly in $[0,1]$ for different noise realizations. The GR quantiles we find are consistent with such a uniform distribution. In particular, the (two-tailed) meta $p$-value for all events and $\alpha$ values obtained using Fisher’s method Fisher (1948) (as in Sec. V.1) is $0.9995$.

We find that the combined bounds overall improve on those quoted in Abbott et al. (2017) by roughly the factor of $\sqrt{7/3}\simeq 1.5$ expected from including more events, with the bounds for some quantities improving by up to a factor of $2.5$, due to the inclusion of several more massive and distant systems in the sample.222222While the results in Abbott et al. (2017) were affected by a slight normalization issue, and also had insufficiently fine binning in the computation of the upper bounds, we find improvements of up to a factor of $3.4$ when comparing to the combined GW150914 + GW151226 + GW170104 bounds we compute here. These massive and distant systems, notably GW170823 (and GW170729, which is not included in the combined results), generally give the best individual bounds, particularly for small values of $\alpha$, where the dephasing is largest at lower frequencies. Closer and less massive systems such as GW151226 and GW170608 provide weaker bounds, overall. However, their bounds can be comparable to those of the more massive, distant events for larger values of $\alpha$. The lighter systems have more power at higher frequencies where the dephasing from the modified dispersion is larger for larger values of $\alpha$.

The new combined bound on the mass of the graviton of $m_{g}\leq 4.7\times 10^{-23}\text{ eV}/c^{2}$ is a factor of $1.6$ improvement on the one presented in Abbott et al. (2017). It is also a small improvement on the bound of $m_{g}\leq 6.76\times 10^{-23}\text{ eV}/c^{2}$ ($90\%$ confidence level) obtained from Solar System ephemerides in Bernus et al. (2019).232323The much stronger bound in Will (2018) is deduced from a post-fit analysis (i.e., using the residuals of a fit to Solar System ephemerides performed without including the effects of a massive graviton). It may therefore overestimate Solar System constraints, as is indeed seen to be the case in Bernus et al. (2019). However, these bounds are complementary, since the GW bound comes from the radiative sector, while the Solar System bound considers the static modification to the Newtonian potential. See, e.g., de Rham et al. (2017) for a review of bounds on the mass of the graviton.

We find that the posterior on $A_{\alpha}$ peaks away from [math] in some cases (illustrated in Sec. A.4 of the Appendix), and the GR quantile is in one of the tails of the distribution. This feature is expected for a few out of 10 events, simply from Gaussian noise fluctuations. We have performed simulations of 100 GR sources with source-frame component masses lying between 25 and 45 $M_{\odot}$, isotropically distributed spins with dimensionless magnitudes up to $0.99$, and at luminosity distances between 500 and 800 Mpc. These simulations used the waveform model IMRPhenomPv2 and considered the Advanced LIGO and Virgo network, using Gaussian noise with the detectors’ design sensitivity power spectral densities. We found that in about 20 – 30% of cases, the GR quantile lies in the tails of the distribution (i.e., $<10\%$ or $>90\%$), when the sources injected are analyzed using the same waveform model (IMRPhenomPv2).

In order to assess the impact of waveform systematics, we also analyze some events using the aligned-spin SEOBNRv4 model. We consider GW170729 and GW170814 in depth in this study because the GR quantiles of the IMRPhenomPv2 results lie in the tails of the distributions, and find that the $90\%$ upper bounds and GR quantiles presented in Table 4 differ by at most a factor of $2.3$ for GW170729 and $1.5$ for GW170814 when computed using the SEOBNRv4 model. These results are presented in Sec. A.4 of the Appendix.

There are also uncertainties in the determination of the 90% bounds due to the finite number of samples and the long tails of the distributions. As in Ref. Abbott et al. (2017), we quantify this uncertainty using Bayesian bootstrapping Rubin (1981). We use 1000 bootstrap realizations for each value of $\alpha$ and sign of $A_{\alpha}$, obtaining a distribution of 90% bounds on $A_{\alpha}$. We consider the 90% credible interval of this distribution and find that its width is $<30\%$ of the values for the 90% bounds on $A_{\alpha}$ given in Table 4 for all but $10$ of the $160$ cases we consider (counting positive and negative $A_{\alpha}$ cases separately). For GW170608, $A_{4}<0$, the width of the 90% credible interval from bootstrapping is $91\%$ of the value in Table 4. This ratio is $\leq 47\%$ for all the remaining cases. Thus, there are a few cases where the bootstrapping uncertainty in the bound on $A_{\alpha}$ is large, but for most cases, this is not a substantial uncertainty.

VIII Polarizations

Generic metric theories of gravity may allow up to six polarizations of gravitational waves Eardley et al. (1973): two tensor modes (helicity $\pm 2$), two vector modes (helicity $\pm 1$), and two scalar modes (helicity 0). Of these, only the two tensor modes ($+$ and $\times$) are permitted in GR. We may attempt to reconstruct the polarization content of a passing GW using a network of detectors Will (2014); Chatziioannou et al. (2012); Isi et al. (2017); Callister et al. (2017); Isi and Weinstein (2017). This is possible because instruments with different orientations will respond differently to signals from a given sky location depending on their polarization. In particular, the strain signal in detector $I$ can be written as $h_{I}(t)=\sum_{A}F^{A}{}_{I}h_{A}(t)$, with $F^{A}{}_{I}$ the detector’s response function and $h_{A}(t)$ the $A$-polarized part of the signal Will (2014); Błaut (2012).

In order to fully disentangle the polarization content of a transient signal, at least 5 detectors are needed to break all degeneracies Chatziioannou et al. (2012).242424Differential-arm detectors are only sensitive to the traceless scalar mode, meaning we can only hope to distinguish five, not six, polarizations. This limits the polarization measurements that are currently feasible. In spite of this, we may extract some polarization information from signals detected with both LIGO detectors and Virgo Isi and Weinstein (2017). This was done previously with GW170814 and GW170817 to provide evidence that GWs are tensor polarized, instead of fully vector or fully scalar Abbott et al. (2017, 2019). Besides GW170814, there are three binary black hole events that were detected with the full network (GW170729, GW170809, and GW170818). Of these events, only GW170818 has enough SNR and is sufficiently well localized to provide any relevant information (cf. Fig. 8 in Abbott et al. (2019)). The Bayes factors (marginalized likelihood ratios) obtained in this case are $12\pm 3$ for tensor vs vector and $407\pm 100$ for tensor vs scalar, where the error corresponds to the uncertainty due to discrete sampling in the evidence computations. These values are comparable to those from GW170814, for which the latest recalibrated and cleaned data (cf. Sec. II) yield Bayes factors of $30\pm 4$ and $220\pm 27$ for tensor vs vector and scalar respectively.252525These values are less stringent than those previously published in Abbott et al. (2017). This is solely due to the change in data, which impacted the sky locations inferred under the non-GR hypotheses. Values from these binary black holes are many orders of magnitude weaker than those obtained from GW170817, where we benefited from the precise sky-localization provided by an electromagnetic counterpart Abbott et al. (2019).

IX Conclusions and outlook

We have presented the results from various tests of GR performed using the binary black hole signals from the catalog GWTC-1 LIGO Scientific Collaboration and Virgo Collaboration (2018), i.e., those observed by Advanced LIGO and Advanced Virgo during the first two observing runs of the advanced detector era. These tests, which are among the first tests of GR in the highly relativistic, nonlinear regime of strong gravity, do not reveal any inconsistency of our data with the predictions of GR. We have presented full results on four tests of the consistency of the data with gravitational waveforms from binary black hole systems as predicted by GR. The first two of these tests check the self-consistency of our analysis. One checks that the residual remaining after subtracting the best-fit waveform is consistent with detector noise. The other checks that the final mass and spin inferred from the low- and high-frequency parts of the signal are consistent. The third and fourth tests introduce parameterized deviations in the waveform model and check that these deviations are consistent with their GR value of zero. In one test, these deviations are completely phenomenological modifications of the coefficients in a waveform model, including the post-Newtonian coefficients. In the other test, the deviations are those arising from the propagation of GWs with a modified dispersion relation, which includes the dispersion due to a massive graviton as a special case. In addition, we also check whether the observed polarizations are consistent with being purely tensor modes (as expected in GR) as opposed to purely scalar modes or vector modes.

We present results from all binary black hole events that are detected with a false alarm rate better than $(1\,\mathrm{yr})^{-1}$. This includes results from the re-analysis of some of the events which were published earlier Abbott et al. (2016, 2016, 2017, 2017), with better calibration and data quality. Assuming that the parameters that describe deviations from GR take values that are independent of source properties, we can combine results from multiple events. We choose to combine results only from highly significant events, detected with a false alarm rate better than $(1000\,\mathrm{yr})^{-1}$ in both modeled searches. Combining together these results has allowed us to significantly reduce the statistical errors on constraints on deviations from GR predictions, as compared to those from individual events. The combined constraints presented here improve our previously presented constraints by factors of {\color[rgb]{0,0,0}1.1} to {\color[rgb]{0,0,0}2.5}, with the largest improvements obtained for certain cases of the modified dispersion test. Notable constraints include that on the graviton’s mass $m_{g}\leq 4.7\times 10^{-23}\text{ eV}/c^{2}$ (an improvement of a factor of $1.6$ over previously presented constraints) and the first constraint on the $-1$PN coefficient obtained from binary black holes.

With the expected observations of additional binary black hole merger events in the upcoming LIGO/Virgo observing runs Abbott et al. (2018, 2019), the statistical errors of the combined results will soon decrease significantly. A number of potential sources of systematic errors (due to imperfect modeling of GR waveforms, calibration uncertainties, noise artifacts, etc.) need to be understood for future high-precision tests of strong gravity using GW observations. However, work to improve the analysis on all these fronts is well underway, for instance the inclusion of full spin-precession dynamics Taracchini et al. (2014); Babak et al. (2017); Blackman et al. (2017); Khan et al. (2019); Varma et al. (2019), non-quadrupolar modes Blackman et al. (2017); London et al. (2018); Cotesta et al. (2018); Varma et al. (2019, 2019), and eccentricity Hinderer and Babak (2017); Cao and Han (2017); Hinder et al. (2018); Huerta et al. (2018); Klein et al. (2018); Moore et al. (2018); Moore and Yunes (2019); Tiwari et al. (2019) in waveform models, as well as analyses that compare directly with numerical relativity waveforms Abbott et al. (2016); Lange et al. (2017). Additionally, a new, more flexible parameter estimation infrastructure is currently being developed Ashton et al. (2019), and this will allow for improvements in, e.g., the treatment of calibration uncertainties or PSD estimation to be incorporated more easily. We thus expect that tests of general relativity using the data from upcoming observing runs will be able to take full advantage of the increased sensitivity of the detectors.

Acknowledgements.

The authors gratefully acknowledge the support of the United States National Science Foundation (NSF) for the construction and operation of the LIGO Laboratory and Advanced LIGO as well as the Science and Technology Facilities Council (STFC) of the United Kingdom, the Max-Planck-Society (MPS), and the State of Niedersachsen/Germany for support of the construction of Advanced LIGO and construction and operation of the GEO600 detector. Additional support for Advanced LIGO was provided by the Australian Research Council. The authors gratefully acknowledge the Italian Istituto Nazionale di Fisica Nucleare (INFN), the French Centre National de la Recherche Scientifique (CNRS) and the Foundation for Fundamental Research on Matter supported by the Netherlands Organisation for Scientific Research, for the construction and operation of the Virgo detector and the creation and support of the EGO consortium. The authors also gratefully acknowledge research support from these agencies as well as by the Council of Scientific and Industrial Research of India, the Department of Science and Technology, India, the Science & Engineering Research Board (SERB), India, the Ministry of Human Resource Development, India, the Spanish Agencia Estatal de Investigación, the Vicepresidència i Conselleria d’Innovació, Recerca i Turisme and the Conselleria d’Educació i Universitat del Govern de les Illes Balears, the Conselleria d’Educació, Investigació, Cultura i Esport de la Generalitat Valenciana, the National Science Centre of Poland, the Swiss National Science Foundation (SNSF), the Russian Foundation for Basic Research, the Russian Science Foundation, the European Commission, the European Regional Development Funds (ERDF), the Royal Society, the Scottish Funding Council, the Scottish Universities Physics Alliance, the Hungarian Scientific Research Fund (OTKA), the Lyon Institute of Origins (LIO), the Paris Île-de-France Region, the National Research, Development and Innovation Office Hungary (NKFIH), the National Research Foundation of Korea, Industry Canada and the Province of Ontario through the Ministry of Economic Development and Innovation, the Natural Science and Engineering Research Council Canada, the Canadian Institute for Advanced Research, the Brazilian Ministry of Science, Technology, Innovations, and Communications, the International Center for Theoretical Physics South American Institute for Fundamental Research (ICTP-SAIFR), the Research Grants Council of Hong Kong, the National Natural Science Foundation of China (NSFC), the Leverhulme Trust, the Research Corporation, the Ministry of Science and Technology (MOST), Taiwan and the Kavli Foundation. The authors gratefully acknowledge the support of the NSF, STFC, MPS, INFN, CNRS and the State of Niedersachsen/Germany for provision of computational resources. The authors would like to thank Clifford Will and Nicolás Yunes for useful discussions.

Appendix A Individual results and systematics studies

In the main body of the paper, for most analyses, we present only the combined results from all events. Here we present the posteriors from various tests obtained from individual events. In addition, we offer a limited discussion on systematic errors in the analysis, due to the specific choice of a GR waveform approximant.

A.1 Residuals test

As mentioned in Sec. V.1, the residuals test is sensitive to all kinds of disagreement between the best-fit GR-based waveform and the data. This is true whether the disagreement is due to actual deviations from GR or more mundane reasons, like physics missing from our waveform models (e.g., higher-order modes). Had we found compelling evidence of coherent power in the residuals that could not be explained by instrumental noise, further investigations would be required to determine its origin. However, given our null result, we can simply state that we find no evidence for shortcomings in the best-fit waveform, neither from deviations from GR nor modeling systematics.

As the sensitivity of the detectors improves, the issue of systematics will become increasingly more important. To address this, future versions of this test will be carried out by subtracting a best-fit waveform produced with more accurate GR-based models, including numerical relativity.

A.2 Inspiral-merger-ringdown consistency test

In order to gauge the systematic errors in the IMR consistency test results due to imperfect waveform modeling, we have also estimated the posteriors of the deviation parameters $\Delta M_{\mathrm{f}}/\bar{M}_{\mathrm{f}}$ and $\Delta a_{\mathrm{f}}/\bar{a}_{\mathrm{f}}$ using the effective-one-body based waveform family SEOBNRv4 that models binary black holes with non-precessing spins. This analysis uses the same priors as used in the main analysis presented in Sec.V.2, except that spins are assumed to be aligned/antialigned with the orbital angular momentum of the binary. The resulting posteriors are presented in Fig. 7 and are broadly consistent with the posteriors using IMRPhenomPv2 presented in Fig. 2. The differences in the posteriors of some of the individual events are not surprising, due to the different assumptions on the spins. For all events, the GR value is recovered in the 90% credible region of the posteriors.

A.3 Parametrized tests of gravitational wave generation

Figures 8 and 9 report the parameterized tests of waveform deviations for the individual events, augmenting the results shown in Fig. 3. A statistical summary of the posterior PDFs, showing median and symmetric 90% credible level bounds for the measured parameters is given in Table 5. Sources with low SNR in the inspiral regime yield uninformative posterior distributions on $\delta\hat{\varphi}_{i}$. These sources are the ones further away and with higher mass, which merge at lower frequencies. For instance, although GW170823 has a total mass close to that of GW150914, being much further away (and redshifted to lower frequencies) makes it a low-SNR event, leaving very little information content in the inspiral regime. The same holds true for GW170729, which has a larger mass. Conversely, low-mass events like GW170608, having a significantly larger SNR in the inspiral regime and many more cycles in the frequency band provide very strong constraints in the $\delta\hat{\varphi}_{i}$ parameters (especially the low-order ones) while providing no useful constraints in the merger-ringdown parameters $\delta\hat{\alpha}_{i}$.

The choice of the $\mathrm{SNR}>6$ threshold explained in Sec. III ensures that most analyses are informative. However, this is not true in all cases, as not all parameters are as easily determined from the data (cf. the good constraints one obtains on the chirp mass with the much weaker constraints on the mass ratio). The two events for which the SNR threshold is insufficient are GW151012 and GW170608, where some post-inspiral parameters are largely unconstrained. The post-inspiral regime is itself divided into the intermediate and merger-ringdown regimes, and for both these events we find the intermediate regime parameters ($\delta\hat{\beta}_{i}$) to be informative; however, the merger–ringdown $\delta\hat{\alpha}_{3}$ for GW170608 and all $\delta\hat{\alpha}_{i}$ for GW151012 extend across the entire prior range considered in the analyses, far beyond the range constrained by other events (as can be seen in Fig. 9). Although we use the results from GW170608 $\delta\hat{\alpha}_{i}$ for combining posteriors in Fig. 3, the combined bounds remain unaffected by adding these results. For future tests, a more discerning threshold than a simple SNR cut, for example including information like the number of cycles of the signal in band, may be used to select which events will provide useful constraints.

Both here and in Sec. VI we report results on the parametrized deviations in the PN regime using two waveform models, IMRPhenomPv2 and SEOBNRv4. There is a subtle difference between the ways deviations from GR are introduced and parametrized in the two models. With IMRPhenomPv2, we directly constrain $\delta\hat{\varphi}_{i}$, which represent fractional deviations in the non-spinning portion of the $(i/2)$PN phase coefficients. The SEOBNRv4 analysis instead uses a parameterization that also applies the fractional deviations to spin contributions, as described in Abbott et al. (2019). The results are then mapped post-hoc from this native parameterization to posteriors on $\delta\hat{\varphi}_{i}$, shown in Figs. 3 and 8 (black solid lines).

In the SEOBNRv4 analysis at $3.5$PN, the native (spin-inclusive) posteriors contain tails that extend to the edge of the prior range. This is due to a zero-crossing of the $3.5$PN term in the $(\eta,a_{1},a_{2})$ parameter space, which makes the corresponding relative deviation ill-defined. After the post-hoc mapping to posteriors on $\delta\hat{\varphi}_{7}$, no tails appear and we find good agreement with the IMRPhenomPv2 analysis, as expected. By varying the prior range, we estimate a systematic uncertainty of at most a few percent on the quoted $90\%$ bounds due to the truncation of tails.

A.4 Parameterized tests of gravitational wave propagation

Posteriors on $A_{\alpha}$ for individual events are shown in Fig. 10, with data for positive and negative $A_{\alpha}$ combined into one violin plot. We provide results for all events with the IMRPhenomPv2 waveform model and also show results of the analysis with the SEOBNRv4 waveform model for GW170729 and GW170814. In Table 6 we compare the 90% bounds on $A_{\alpha}$ and GR quantiles obtained with IMRPhenomPv2 and SEOBNRv4 for GW170729 and GW170814. We focus on these two events because the GR quantiles obtained with IMRPhenomPv2 lie in the tails of the distributions, and we find that this remains true for most $\alpha$ values in the analysis with SEOBNRv4. For GW170729 and $\alpha\in\{0,0.5,1\}$, the GR quantiles obtained using the two waveforms differ by factors of $\sim 2$; the two waveforms give values that are in much closer agreement for the other cases.

Additionally, for the GW151012 event and certain $\alpha$ values, a technical issue with our computation of the likelihood meant that specific points with relatively large values of $A_{\alpha}$ had to be manually removed from the posterior distribution. In particular, for computational efficiency, the likelihood is calculated on as short a segment of data as practical, with duration set by the longest waveform to be sampled. Large values of $A_{\alpha}$ yield highly dispersed waveforms that are pushed beyond the confines of the segment we use, causing the waveform templates to wrap around the boundaries. This invalidates the assumptions underlying our likelihood computation and causes an artificial enhancement of the SNR as reported by the analysis. As expected, recomputing the SNRs for these points on a segment that properly fits the waveform results in smaller values that are consistent with noise. Therefore, we exclude from our analysis parameter values yielding waveforms that would not be contained by the data segment used, which is equivalent to using a stricter prior on $A_{\alpha}$. Failure to do this may result in the appearance of outliers with spuriously high likelihood for large values of $A_{\alpha}$, as we have seen in our own analysis.