Testing General Relativity with Gravitational Waves: An Overview

Abstract

The detections of gravitational-wave (GW) signals from compact binary coalescence by ground-based detectors have opened up the era of GW astronomy. These observations provide opportunities to test Einstein’s general theory of relativity at the strong-field regime. Here we give a brief overview of the various GW-based tests of General Relativity (GR) performed by the LIGO-Virgo collaboration on the detected GW events to date. After providing details for the tests performed in four categories, we discuss the prospects for each test in the context of future GW detectors. The four categories of tests include the consistency tests, parametrized tests for GW generation and propagation, tests for the merger remnant properties, and GW polarization tests.

1. Introduction

The binary evolution in General Relativity (GR) is described differently than in Newtonian gravity (NG). In GR, the binary orbit shrinks due to the emission of energy: angular and linear momenta through gravitational waves (GWs) [1,2,3,4]. Whereas in NG, there is no concept of radiation reaction and the orbital period is constant over time. Even though Albert Einstein predicted the existence of GWs more than a century ago, their detection remained a puzzle due to their weak interaction with matter. The indirect evidence of GWs came from the decades-long observations of orbital decay of a binary pulsar by Russell Alan Hulse and Joseph Taylor [5,6,7,8]. They found that the observed orbital decay of the binary system, known as PSR B1913+16(PSR J1915+1606, or PSR 1913+16), due to the emission of GWs, is consistent with the predictions of GR. That is, the rate of decay of the orbital period ($P_{orb}$) from observations ${\dot{P}}_{orb}\sim {10}^{-14}$–${10}^{-12}$ agreed to the GR predicted rate obtained from analytical calculations based on GR, leading the team to win the Physics Nobel Prize in 1993.

The direct detection of GWs had to wait until the LIGO detectors at Hanford, Washington, and Livingston, Louisiana, made their first detection of a binary merger on 14 September 2015 [9,10,11,12,13,14,15,16,17]. This discovery opened the era of GW astronomy, noting the first highly relativistic strong-field observation of GWs. Within the subsequent years of observation runs, the LIGO-Virgo collaborations announced the detection of more than fifty binary merger events [12,18,19,20,21,22,23,24,25,26,27,28,29].

Among the many significant contributions to fundamental physics and astrophysics, GW observations test GR at the relativistic, strong-field regime. A set of testing GR analyses conducted by the LIGO-Virgo Scientific Collaboration (LVC) on the GW150914 event established that GW150914 is consistent with a binary black hole (BBH) signal described in GR [9].

This set of tests include consistency tests, parameterized tests, tests to confirm the non-dispersive nature of the radiation, and tests on the remnant properties. Consistency tests check for the agreement of observed data with the signal predicted from GR. Parameterized tests introduce model-agnostic parametric GR deviations in the waveform and constrain those from the data to put statistical bounds on these parameters. The list of the GW events was extended further with more binary merger detections by LIGO and Virgo detectors’ first, second, and third observation runs. All the tests applied on GW150914 were also performed on these events with appropriate modifications to the above-mentioned tests, including additional tests. Here we will go through them in detail.

The GW-based tests of GR on the BBH coalescence events detected by LIGO and Virgo until 1 October 2019 are available in Reference [30]. A generic binary system evolves from its early inspiral weak-field regime to a highly relativistic merger and then the final ringdown stage. In the case of BBHs in GR, the object formed after they merge (i.e., the merger remnant) is another black hole (BH). On the other hand, for non-BH binaries [31,32],1 the merger remnant is not necessarily a BH but could be another compact star depending upon the properties of the binary.

In a model agnostic way, there were four broad classes of tests conducted in Reference [30]. These tests aim to look at different regimes of binary evolution or to the full inspiral-merger-ringdown signals. The first set consists of the residual analysis and the inspiral-merger-ringdown consistency test. Both of these tests check the consistency of GR predictions with the observed data (as in the case of Reference [9]). The second category of tests is the parametrized tests for GW generation and propagation. Here one sets statistical bounds on the parametrized deviations from GR, assuming GR is the correct theory of gravity, employing GR waveform models with parametric deviations present. On the third category of tests, one looks for any violation of GR by analyzing the merger remnant properties. The GW polarization tests look for extra polarization modes present in the data and comprise the fourth set of tests. This analysis provides statistical evidence for alternative theories of gravity that predict vector and scalar polarization modes along with the tensor modes. An overview of these tests is provided in Figure 1.

The probability of there being astrophysical origin of a candidate event plays an important role in determining whether that event is considered for the testing GR analyses or not. Usually, a higher threshold is assumed so that the events analyzed have higher chances of being of astrophysical origin. For instance, in Reference [30], events satisfying a false alarm rate (FAR) less than ${10}^{-3}$ per year are chosen to analyze. Once the set of events is chosen based on the detection significance, additional criteria are applied depending upon the strategies followed by each test. Bayesian formalism-based techniques are employed to get meaningful bounds from each test. The pipelines widely used for this purpose are, LALInference [48] available in the LIGO Scientific Collaboration’s algorithm library suite (LALSuite) [49], Bayeswave [50,51], parallel bilby (pBilby) [52,53,54], PYRING [55,56] and Bantam [57]. Reference [30] demonstrated the possibility of performing tests of GR on binary black hole (BBH) events, employing mainly two different waveform models, IMRPhenomPv2 (phenomenological waveform model for a precessing BBH system) [58,59,60] and SEOBNRv4_ROM (reduced-order effective one body (EOB, waveform model for a non-precessing binary system) [61].

There is a significant increase in the detection rate as the detectors improve their sensitivities through first, second, and third observing runs of LIGO-Virgo detectors. Interestingly, it is possible to infer information from multiple events by combining the data from each event. The combined bounds help to improve our understanding of binary population properties in general. As we combine results, the statistical uncertainty that arises due to instrumental noise lessens. Notice that this instrumental noise does not include the uncertainty contributions from the systematic errors of gravitational waveform modeling [62,63,64,65,66,67,68], calibration of the detectors, and power spectral density (PSD) estimation uncertainties [69,70,71,72,73]. Sometimes, systematic errors can dominate the statistical errors and lead to false identification of GR violations, which we do not discuss here. Previous studies in References [9,68,74] discuss two different statistical approaches to estimate combined information on GR test parameters from multiple events. The first one (also called restricted or simple combining) assumes equal GR deviations across all the events independent of the physical parameters characterizing the binary, and this technique is well described and demonstrated for GWTC-1 events [74]. This assumption is generally incorrect as there are cases when the waveform model can arbitrarily deviate from GR depending upon the binary source properties. The second method, the hierarchical combining strategy, tries to overcome the issue of universality assumption by relaxing it. In this case, instead of assuming uniform GR deviation for all events, a Gaussian distribution models the non-GR parameter. The statistical properties (mean ($\mu$) and standard deviation ($\sigma$)) of this distribution are obtained from the data itself, and the estimates are different for different models of gravity. We call the parameters $\mu$ and $\sigma$ hyperparameters. If GR is the correct theory, the Gaussian distribution should center around zero. The astrophysical population properties of sources play a crucial role in estimating these statistical quantities [30].

The organization of this article is as follows. Section 2 is dedicated to various tests of GR performed already on the GW events observed by the LIGO-Virgo detectors, including tests of consistency with GR (Section 2.1), parametrized tests (Section 2.2), tests based on the merger remnant properties (Section 2.3), and tests for GW polarizations (Section 2.4). We conclude with a summary section, Section 3.

2. Model-Agnostic Tests of General Relativity from Gravitational-Wave Observations

2.1. Tests of Consistency with General Relativity

Consistency tests do not need to assume any particular alternative theory to GR, nor do they test specific deviations. They address the simple question: can the observed data be fully explained by assuming GR? Or put differently, is there any statistically significant “trace” in the data that is unlikely to be explained as either part of an astrophysical signal (assuming GR) or the terrestrial instruments’ noise? So far two different kinds of consistency tests have been performed on the detected GW events [30]: the residual analysis and the Inspiral-merger-ringdown (IMR) consistency test.

2.1.1. Residual Test

The residual test checks for signatures left in the data after subtracting the best-fit GR template. If GR is the correct theory and we have subtracted the astrophysical signal completely, the residuals in each detector should be consistent with instrumental noises.

The best fit model of the astrophysical signal is obtained by a detailed parameter estimation analysis using a stochastic sampling of the signal’s parameter space. Typically, the best-fit parameters are taken to be those that maximize the likelihood of observing the recorded data assuming this signal is present in the data. This set of parameters is not necessarily the one that describes the most probable source configuration a posteriori, as the likelihood alone does not take prior assumptions into account. Nevertheless, the maximum-likelihood parameters are those that minimize the difference between the data $d_i$ and template $h_i$ by definition of the likelihood $\mathsf{\Lambda }$,

$$log\mathsf{\Lambda }=-\frac{1}{2}\sum _i{\parallel d_i-h_i\parallel }^2.$$

(1)

Here, the index i enumerates the different detectors; $h_i$ is the signal projected onto each detector, respectively, and the norm ${\parallel ·\parallel }^2=\langle ·,·\rangle$ is induced by the following inner product

$$〈a,b〉=4\Re \int \frac{\tilde{a}\left(f\right)\tilde{b}{\left(f\right)}^{*}}{S_n\left(f\right)}df.$$

(2)

The detector noise spectral density $S_n\left(f\right)$ acts as a weight in an integral over the Fourier transformed functions $\tilde{a}$ and $\tilde{b}$; ${}^{*}$ denotes complex conjugation.

By this construction, the residuals $d_i-h_i$ are small in the sense of Equation (1), but they could still contain a coherent signal that cannot be captured by the GR model. To look for such a potential signal, the method employed in References [9,30,68] is Bayeswave [50]: a transient search algorithm looking for coherent excess power in the (residual) detector data. This part of the analysis is model-independent. Bayeswave uses Morlet–Gabor wavelets to look for coherent, elliptically polarized signatures that rely on no further model assumption. In addition to generic signals, it employs models for stationary and non-stationary noise simultaneously.

To quantitatively explain the results, we require various definitions. First, the optimal network signal-to-noise-ratio (SNR) of a signal h is derived from its norm,

$$\mathrm{SNR}\left(h\right)=\sqrt{\sum _i\langle h_i,h_i\rangle }=\sqrt{\sum _i{\parallel h_i\parallel }^2},$$

(3)

using the inner product Equation (2). If we take h as the best-fit GR template, we obtain ${\mathrm{SNR}}_{\mathrm{GR}}$. The residual modelled by Bayeswave is not a single signal. The uncertainty in what constitutes a coherent residual signal and what is instrument noise leads Bayeswave to provide a discretized probability distribution in the parameter space of wavelets. However, each point in this distribution corresponds to a residual signal that has a well-defined SNR. Consequently, we can map the probability distribution of residual signals to a probability of their optimal network SNRs.

As is standard, we characterize the probability distribution by credible intervals that enclose a certain amount of probability (we use this quantity more frequently in this article). Specifically, we report the SNR of the residual at which the cumulative probability distribution is 90%. Put differently, we infer a 90% probability that the residual signal after subtracting the best-fit GR template has an optimal network SNR $\le {\mathrm{SNR}}_{90}$.

The left panel of Figure 2 shows ${\mathrm{SNR}}_{90}$ as a function of the best-fit ${\mathrm{SNR}}_{\mathrm{GR}}$ for the observed binary mergers from O1, O2, and O3a. The SNR of the GR signals ranges between 9.24 (GW151012) and 25.71 (GW190521). The 90% upper credible bound of the residual SNR ranges between 4.88 (GW190727_060333) and 9.24 (GW170818). No clear correlation is visible between the two quantities. If our GR models would consistently be unable to capture an ubiquitous deviation from GR, we might expect that stronger signals correlate with stronger residuals. The current data shows no indications of such correlations.

To systematically assess if the residual SNRs are consistent with the detector noise, we can define a p-value under the hypothesis that the residual is consistent with detector noise (this is the null hypothesis). This p-value is estimated by running identical Bayeswave analyses on a large number of noise-only data segments around, but not including, each event. The p-value then provides the probability of pure noise producing an ${\mathrm{SNR}}_{90}^n$ greater than or equal to the residual SNR found after subtracting the best-fit GR template, ${\mathrm{SNR}}_{90}$. That implies, p-value: = $P({\mathrm{SNR}}_{90}^n\ge {\mathrm{SNR}}_{90})$. A large p-value indicates that there is a high chance that the residual power originates from the instrumental noise. Small p-values, on the other hand, indicate that it is less likely for noise alone to yield such high values of residual SNR.

The p-values for all events considered here are included in Figure 2 as a color scale. They also span a large range between $0.07$ (GW90421_213856) and $0.97$ (GW190727_060333). This is to be expected in repeated uncorrelated experiments. In fact, assuming the residuals are pure noise, the p-values of the population should be uniformly distributed between $[0,1]$. As further discussed and quantified in [30], the p-values found for O1, O2, and O3a events are consistent with this expectation.

As a final interpretation of the residuals, one can ask: how well does the GR model fit the signal in the data? Obviously, if the model would be perfect, and we could unambiguously identify the coherent signal in the data, then the agreement between the model and data should be perfect, too, assuming GR is the correct theory. In reality, our GR models are very accurate, but not perfect, and we only have a probabilistic measure of the signal in the data. Therefore, we can only expect to obtain a lower bound on the fitting factor $\mathrm{FF}$ between the model h and the signal s by calculating

$$\begin{array}{cc}\hfill \mathrm{FF}& =\frac{\langle h,s\rangle }{\parallel h\parallel \parallel s\parallel }=\frac{\langle h,h+s_r\rangle }{\parallel h\parallel \parallel h+s_r\parallel }=\frac{\parallel h\parallel }{\sqrt{{\parallel h\parallel }^2+{\parallel s_r\parallel }^2}}\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \Rightarrow {\mathrm{FF}}_{90}& =\frac{{\mathrm{SNR}}_{\mathrm{GR}}}{\sqrt{{\mathrm{SNR}}_{\mathrm{GR}}^2+{\mathrm{SNR}}_{90}^2}}.\hfill \end{array}$$

(5)

Here we used that the coherent signal s is the sum of the GR model h and any residual $s_r$ that is perpendicular to h. The latter assumption is justified because h was chosen by maximizing the agreement between the data and the model. ${\mathrm{FF}}_{90}=1$ would indicate perfect agreement.

We plot $1-{\mathrm{FF}}_{90}$ on the right panel of Figure 2 for the observations we considered so far. They show more clearly a correlation with ${\mathrm{SNR}}_{\mathrm{GR}}$. Strong signals with large values of ${\mathrm{SNR}}_{\mathrm{GR}}$ tend to yield higher fitting factors. This is not because the models describe the actual signal better for louder events. It is because our confidence in how well the model agrees with the actual signal increases with increasing ${\mathrm{SNR}}_{\mathrm{GR}}$.

Another notable fact about the fitting factors is that they are larger than related quantities one often finds in the GW literature. For example, discrete template banks for GW searches are often constructed such that fitting factors of at least $0.97$ are guaranteed between any signal and the closest template in the bank [75]. Waveform models for BBHs are commonly tuned to $\mathrm{FF}\lesssim {10}^{-3}$ between the most accurate predictions and the full model [76,77]. Parameter estimation poses strict demands on waveform accuracy of the order of ${\mathrm{SNR}}_{\mathrm{GR}}^{-2}$ [78,79] (i.e., waveform differences of $<{10}^{-2}$). Is that level of accuracy necessary, given that for most observations, we cannot put stronger constraints on the fitting factor than ${\mathrm{FF}}_{90}\sim \mathcal{O}\left({10}^{-1}\right)$? The answer is, of course, yes! While we cannot be sure about the true signal for individual events, it is worth emphasizing that the Bayeswave analysis has great freedom and a large parameter space to identify virtually any residual signature as a potential coherent signal. Much more restricted measurements that only look for specific, lower-dimensional deviations are sensitive to significantly smaller signal differences, because only those signal differences that are consistent with both the assumed variation and the data are considered. Therefore, the residuals test is a very generic baseline test for anything that cannot be modeled with GR. It is, however, much less sensitive to specific deviations that can be tested more accurately in a dedicated test (see the rest of this paper). A more detailed discussion of the relation between various tests of GR can be found in [80].

2.1.2. Inspiral-Merger-Ringdown Consistency Test

In GR, the time evolution of BBH mergers is uniquely determined. Hence the final mass and spin of the remnant BH are uniquely determined from the initial mass and spin parameters of Kerr BHs. The inspiral-merger-ringdown consistency test (IMRCT) tests the consistency of inspiral and merger-ringdown parts of the signal by comparing two independent estimates of binary parameters. More specifically, the binary’s final mass and spin parameters are measured separately from both low- and high-frequency parts of the signals and then compared to the two measurements to check their agreement.

Given the final mass $\left(M_f\right)$ and spin $\left({\chi }_f\right)$, one can estimate the spin-dependent innermost stable circular frequency2 for a Kerr BH ($f_{\mathrm{cut}}$) [81,82,83]. The full BBH signal can be divided into two parts using this frequency, the low-frequency part and the high-frequency part3. By restricting the noise-weighted integral in the likelihood calculation to frequencies below $(f<f_{\mathrm{cut}})$ and above this frequency cutoff $(f>f_{\mathrm{cut}})$, the binary parameters are estimated using stochastic sampling algorithms based on Bayesian inference. The merger remnant properties are calculated by averaging NR-calibrated final state fits given the posterior median values [81,82,83] from the two independent mass-spin estimates above. This calculation assumes an aligned-spin binary system. If the data is consistent with GR, both estimates should agree [12,20,84,85].

The frequency, $f_{\mathrm{cut}}$, roughly divides the signal into inspiral and merger-ringdown (post-inspiral) regimes. To calculate $f_{\mathrm{cut}}$, the binary parameters inferred from the full signal are used. As the test relies on independent parameter inference from the low and high frequency parts of the full signal, one requires enough SNR in both these regimes of the signal. For the selected events, a detailed parameter estimation analysis is performed in Reference [30] focusing on the mass-spin parameters. If $M_f^{insp}$ and $M_f^{post-insp}$ denote the final mass estimates obtained from low and high frequency parts of the signal, we can define a dimension-less quantity that measures the fractional deviation from the final-mass estimate as,

$$\delta M_f=\frac{\Delta M_f}{\overline{M_f}}=2\frac{M_f^{insp}-M_f^{post-insp}}{M_f^{insp}+M_f^{post-insp}},$$

(6)

where subscript, ‘f’ denotes merger remnant parameters, and ‘insp’ and ‘post-insp’ correspond to the estimates coming from the low and high-frequency regimes, respectively. $\overline{M_f}$ and $\overline{{\chi }_f}$ are the symmetric combinations of $M_f$ and ${\chi }_f$ estimates from inspiral and post-inspiral regimes. A similar expression can be written for dimension-less spin parameter,

$$\delta {\chi }_f=\frac{\Delta {\chi }_f}{\overline{{\chi }_f}}=2\frac{{\chi }_f^{insp}-{\chi }_f^{post-insp}}{{\chi }_f^{insp}+{\chi }_f^{post-insp}}.$$

(7)

In principle, these fractional deviations vanish if the data is consistent with GR4. However, one must perform a detailed parameter estimation analysis and estimate the statistical confidence that the GR deviation vanishes. This is illustrated in Figure 3 and see Reference [86] for more details of the method and demonstration on simulated binary signals. See References [87,88,89] for studies projecting the possibilities of IMR consistency tests from combining information from current and future detectors. Especially, Reference [88] found that the multi-band observations can improve the constraints by a factor of 1.7.

In Reference [30], from the list of events satisfying the detection criteria, based on the false alarm probability of each event and the requirement of enough SNR in both inspiral and post-inspiral regimes, posteriors distributions on $\delta M_f$ and $\delta {\chi }_f$ are obtained assuming uniform priors on these parameters. In terms of the two-dimensional GR quantile, $Q_{GR}$ is the fraction of the posteriors enclosed by the iso-probability contours that contain the GR value. Reference [30] reports GW190814 as the most consistent event with the quantile $Q_{GR}=99.9\%$.

From the hierarchical combining method, the hyperparameters describing the Gaussian distribution are estimated to be $(\mu ,\sigma )=(0.{02}_{-0.09}^{+0.11},<0.17)$ for $\delta M_f$ and $(\mu ,\sigma )=(-0.{06}_{-0.16}^{+0.15},<0.34)$ for $\delta {\chi }_f$ within the 90% confidence interval in Reference [30]. The details can be found in Figure 4. Assuming that the fractional deviations take the same value for all events, a less-conservative combined 90% confidence interval of $\frac{\Delta M_f}{\overline{M_f}}=0.{04}_{-0.06}^{+0.08}$ and $\frac{\Delta {\chi }_f}{\overline{{\chi }_f}}=-0.{09}_{-0.08}^{+0.11}$ obtained in Reference [30]. This analysis employed IMRPhenomPv2 or IMRPhenomPv3HM waveform models depending upon the information about the higher-mode content present in the binary signal [90] and so far the analysis finds all events to be consistent with GR.

2.2. Parametrized Tests of GR Based on Generation and Propagation of GWs

The parametrized tests are designed to capture any deviations from GR in the generation and propagation of GWs. This is achieved by introducing model-independent parametric variations in the gravitational waveform models and constrain those from the observed data. If GR is the correct theory, parametric deviations vanish, and the statistical bounds can be used to put constraints on the alternative theory models.

2.2.1. Constraining the Parametrized Deviations from General Relativistic Inspiral-Merger-Ringdown Coefficients

Any generic deviation from GR may modify the binary dynamics and its time evolution. This leads to measurable modifications to the equation of motion through the energy and angular momentum of the source and the flux. However, the inspiral-merger-ringdown dynamics are uniquely determined and well studied in GR through various techniques such as post-Newtonian theory, numerical relativity, and BH perturbation theory once we fix the intrinsic parameters of the binary system [58,91,92,93,94,95,96,97,98,99,100].

The inspiral coefficients are modeled analytically using post-Newtonian (PN) theory, which finds perturbative solutions to the binary evolution in terms of a velocity parameter, $v/c$, in the slow-motion limit ($v<<c$, v is the PN parameter and c is the velocity of light). Parametrized tests based on inspiral coefficients are investigated in detail [98,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116] and also demonstrated the applicability of the test using Bayesian framework [9,20,117,118,119] in the past. Moreover, the possibility of constraining these parameters employing multiband observations has been studied in References [87,88,89,120]. For an inspiralling compact binary system, the GW waveform can be schematically represented in the frequency domain as,

$$\tilde{h}\left(f\right)=A\left(f\right)e^{i\varphi \left(f\right)},$$

(8)

where $A\left(f\right)$ denotes the amplitude and $\varphi \left(f\right)$ is the phase of the signal. Now, we introduce parametric deviations of the form, ${\varphi }_i\left(f\right)\to (1+{\phi }_i){\varphi }_i$. If GR is the correct theory, ${\widehat{\varphi }}_i$ vanishes for the Nth PN order, where $i=N/2$ denotes the Nth PN order. The parametrized tests for post-Newtonian coefficients (pPN analysis) measure these deviation parameters and are one of the main tests of GR analyses performed for the detected binary signals so far.

For a generic binary system, one needs to put bounds on the inspiral-merger-ringdown parametrized deviation coefficients separately. In this case, a relative deviation is introduced to each coefficient appearing through the inspiral-merger-ringdown regimes as,

$$p_i\to (1+\delta p_i)p_i.$$

(9)

The set of free parameters, $\delta p_i$ include the inspiral coefficients $\left\{{\phi }_i\right\}$ and post-inspiral coefficients $\{{\alpha }_i,{\beta }_i\}$ [62,106,121]. It is not plausible to represent the post-inspiral coefficients analytically and they are obtained by numerical fits.

The IMRPhenomPv2 waveform model [58,59,122] describes a precessing binary system in frequency domain with inspiral coeffiecients determined by PN theory and intermediate, and merger–ringdown regions by finding appropriate numerical fits. The transition frequency between the inspiral to merger-ringdown is defined as $GM(1+z)f_c^{PAR}/c^3=0.018$ (M is the total mass of the binary system and z is the redshift to the binary). The results for parametrized tests for post-Newtonian coefficients reported in Reference [30] relied on the IMRPhenomPv2 waveform model which allows for parametrized deviations of the phenomenological coefficients describing the inspiral, intermediate ${\beta }_i=\{{\beta }_1,{\beta }_2\}$, and merger-ringdown ${\alpha }_i=\{{\alpha }_2,{\alpha }_3,{\alpha }_4\}$ regions.

There is also another equally accepted method to perform the same analysis, which is not based on any particular waveform model but rather on theory-agnostic modifications applied to the inspiral coefficients of any waveform model. These low-frequency modifications (inspiral-only modifications $\left\{{\phi }_i\right\}$) are tapered to zero at high frequencies. That is, as the frequency reaches post-inspiral regions these modifications vanish and the signal agrees to a BBH signal described in GR. This test is carried out employing SEOBNRv4_ROM waveform model in References [61,118,123]. One of the main differences between the two approaches is that the deviations are applied to only the non-spinning coefficients for the first method. Still, for the second method, the aligned-spin waveform coefficients are modified. Both approaches provide consistent results when we compare them.

Denoting $\left\{{\phi }_i\right\}$ as the deviation from $N=i/2$ PN order, the list of inspiral coefficients we can put bounds from the data are,

$$\{{\phi }_{-2},{\phi }_0,{\phi }_1,{\phi }_2,{\phi }_3,{\phi }_4,\delta {\phi }_5,{\phi }_6,{\phi }_7\}.$$

(10)

This means the coefficients up to 3.5PN order ($i=7$) are available, including the logarithmic terms at 2.5 and 3PN orders. Due to its degeneracy with the coalescence phase, one cannot constrain the coefficient at 2.5 PN (term having no logarithmic dependence). Notice also that ${\phi }_{-2}$ is zero in GR and, so as ${\phi }_1$, it is introduced to account for specific alternative theories of gravity, especially those that predict dipolar radiation. For the two coefficients, ${\phi }_1$ and ${\phi }_{-2}$, we have the absolute deviations while all other parameters provide relative deviations from the respective GR coefficients. The absolute deviations are the differences between the true/actual value and the measured value, whereas the relative deviations are the ratio of absolute deviations to the true/actual value.

Among all the coefficients listed above, the best combined bound is obtained for the Newtonian coefficient in Reference [30], $|{\phi }_0|\le 4.4\times {10}^{-2}$ (neglecting the −1PN coefficient, ${\phi }_{-2}$). Notice that this bound is weaker than the bounds from the double pulsar measurement by a factor of ∼3. Reference [30] reports that the posterior on $\delta \widehat{p_i}$ is consistent with the GR prediction within the 90% confidence interval for all the events considered. From the hierarchical analysis, the tightest (loosest) bound obtained is ${\phi }_{-2}=0.{97}_{-4.07}^{+4.62}\times {10}^{-3}$ (${\phi }_{6\ell }=0.{42}_{-1.50}^{+1.67}$). For both the circumstances, the GR hypothesis is preferred with quantiles $Q_{GR}=68\%$ ($Q_{GR}=69\%$, which is close to the median values for both the cases. Figure 5 shows the combined posterior distributions obtained from the GWTC-2 events considering both the hierarchical method and a restricted method, where the deviation parameters do not allow for variance between events. Results from both the studies described above are found to be consistent with GR.

2.2.2. Tests of BBH Nature from Spin-Induced Quadrupole Moment Measurements

Spin-induced multipole moments arise due to the spinning motion of the compact object and take unique values for BHs given mass and spin. For the BBH signals, these effects are included in the post-Newtonian modeling, and they appear along with the spin-spin terms and can be schematically represented as,

$$Q=-\kappa {\chi }^2{m}^3,$$

(11)

here $\kappa$ is the spin-induced quadrupole moment coefficient, m is the mass, and $\chi$ is the dimensionless spin parameter. For BHs, ${\kappa }_{BH}=1$ from the no-hair conjecture [124,125,126] and for any other compact objects, the value may vary depending upon the properties of the star. Through numerical relativity simulations of slowly spinning neutron stars, it is found that the value of $\kappa$ varies between ${\kappa }_{NS}=$ 2 and 14 [127,128,129]. On the other hand, for more exotic stars like boson stars, the value of $\kappa$ can be even larger and found to vary between ∼10 and 100 [40,130].

It has been shown that one can introduce parametrized deviations of the form, $\kappa =(1+\delta \kappa )$ and put bound on $\delta \kappa$ [131,132,133,134,135]. An inspiralling binary system is parametrized by two such parameters, corresponding to both the binary components, $\delta {\kappa }_1$ and $\delta {\kappa }_2$. Simultaneous measurement of both these parameters will end up giving weak constraints on either parameter, hence it is proposed to measure $\delta {\kappa }_s=0.5(\delta {\kappa }_1+\delta {\kappa }_2)$ keeping $\delta {\kappa }_a=0.5(\delta {\kappa }_1-\delta {\kappa }_2)=0$. This is a safe assumption if we are testing the BBH nature of the detected signal [136,137]. The analysis performed on the first, second, and third observing runs of LIGO-Virgo detectors provided good constraints on events with non-zero spins which include GW151226, GW190412, GW190720_000836, and GW190728_064510. Employing IMRPhenomPv2 waveform model [58,59,122] and assuming prior distribution on $\delta {\kappa }_s$ ranging uniform between [−500, 500] [30]. These events, where the posteriors are very different from the prior knowledge, are highlighted in Figure 6.

With the restricted assumption that $\delta {\kappa }_s$ take the same value for all events, Reference [30] reported a combined bound on $\delta {\kappa }_s$ within the $90\%$ confidence interval as $\delta {\kappa }_s=-15.2_{-19.0}^{+16.9}$. Using the hierarchical analysis, the hyperparameters are constrained to $\mu =-24.6_{-35.3}^{+30.7}$ and $\sigma <52.7$ with $\delta {\kappa }_s=-23.2_{-62.4}^{+52.2}$ and which is again consistent with the null ($\mu =\sigma =0$) hypothesis at $90\%$ confidence interval. The hypothesis stating that the population contains all BBHs is favored by the population containing all the non-BBH hypothesis by a combined Bayes factor of 11.7. The analysis found that the data are consistent with the BBH hypothesis.

2.2.3. Tests of Gravity from GW Propagation

GW propagation in GR is non-dispersive and described by the dispersion relation,

$$E^2=p^2{c}^2.$$

(12)

Equivalently, the velocity of propagation of GWs in GR is independent of the frequency of the radiation. As a consequence, the graviton is massless with a corresponding infinite Compton wavelength. There are alternative theories of gravity predicting GWs with dispersion where the local Lorentz invariance is not respected [138].

For a generic theory of gravity, the GR dispersion relation may require modifications and the following equation can account for such propagation effects [139,140,141,142],

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

(13)

We can re-parametrize the gravitational waveforms so that they also account for the propagation effects given in Equation (13). In Equation (13), $\mathcal{A}$ is the dispersion amplitude and has dimension of ${\left[Energy\right]}^{2-\alpha }$, and $\alpha$ is a dimensionless constant. These paramtetrized modifications can be constrained from the data and these bounds can be translated into constraints on different alternative gravity models. For example, $\alpha =0$ and $\mathcal{A}>0$ correspond to massive graviton theories [141], $\alpha =2.5$ corresponds to multifractional spacetime [142], and $\alpha =3$ corresponds to double special relativity [140], etc.

As shown in References [67,139], one can use GW observations to get constraints on the modified dispersion parameters. The first bound on the Compton wavelength (which has a finite value for any massive graviton theory) from GW observations of a BBH signal is ${\lambda }_g>{10}^{13}\mathrm{km}$ [9] and this has been extended to more generic cases in the subsequent analyses [20,22,74]. This bound on the Compton wavelength translates to a graviton mass, $m_g\le 5\times {10}^{-23}\mathrm{eV}/{\mathrm{c}}^2$, and this is a stronger bound compared to the solar system constraints [74]. From the GWTC-2 data [30], a factor of 2.7 improvements is observed on this bound, and the graviton mass bound correspondingly changes to $m_g\le 1.76\times {10}^{-23}{\mathrm{eV}/\mathrm{c}}^2$ with $90\%$ credibility.

Note that the past studies have investigated the possibility of dispersion of GWs described by GR due to specific physical effects. For example, the nonlinear interaction between charged particles and GWs may lead to dispersion of GWs when GWs pass through astrophysical plasma in the presence of magnetic fields [143,144,145,146,147,148,149,150,151,152].

2.3. Tests Based on the Merger Remnant Properties

The merger remnant of a BBH system emits GWs to settle down to the stationary state, and this distinct signal from stellar-mass BBHs can be measured using the current GW detectors. Consequently, many tests for remnant nature have been proposed and performed on the detected GW events. We briefly discuss GW tests based on merger remnant properties here.

2.3.1. No-Hair Theorem Based Tests from the Quasi-Normal Mode Ringdown Radiation (BH Spectroscopy)

According to GR, the remnant formed after a BBH coalescence is a perturbed Kerr BH, and this BH attains the stationary state by emitting GWs. This damped sinusoidal signal (BH ringdown radiation) is characterized by quasi-normal-modes (QNMs) with frequency f and damping time $\tau$. Both the damping time and frequency of this oscillation are determined by the mass and spin of the Kerr BH formed after the merger ($M_f$ and ${\chi }_f$) [125,153,154,155,156,157]. In GR, the ringdown waveform is a superposition of damped sinusoids and takes the form,

$$\begin{array}{cc}\hfill h_{+}\left(t\right)-i{h}_{\times }\left(t\right)=\sum _{\ell =2}^{+\infty }\sum _{m=-\ell }^{\ell }\sum _{n=0}^{+\infty }& {\mathcal{A}}_{\ell mn}exp\left[-\frac{t-t_0}{(1+z){\tau }_{\ell mn}}\right]exp\left[\frac{2\pi i{f}_{\ell mn}(t-t_0)}{1+z}\right]{}_{-2}{S}_{\ell mn}(\theta ,\varphi ,{\chi }_{\mathrm{f}}),\hfill \end{array}$$

(14)

where z is the cosmological redshift, and the $(l,m,n)$ indices label the QNMs $(\ell ,m)$ are the angular multipoles, whereas n is the order of modes given $(\ell ,m)$. All the $f_{\ell mn}$ and ${\tau }_{\ell mn}$ are determined by the final mass and spin of the binary system (this is called the final state conjecture). For a perturbed Kerr BH, the damping time and frequency of each quasi-normal-mode can be calculated from BH perturbation theory as a function of its final mass and spin [158,159,160,161]. Assuming that the mergers we observe are BBHs, from the independent $M_f$ and ${\chi }_f$ post-merger measurements we can test the final state conjecture (commonly known as the BH spectroscopy) [55,56,83,162,163,164,165,166,167,168,169]. The complex amplitude ${\mathcal{A}}_{\ell mn}$ is a measure of the mode excitation and the phase of these modes at a reference time [170,171,172].

PYRING is a toolkit to perform BH spectroscopy which is completely implemented in the time domain [55,56]. As both templates and the likelihood are modeled in the time domain, spectral leakage is reduced [173]. Mainly assumed template models for this study are, $Ker{r}_{220}$ ($\ell =\left|m\right|=2,n=0$ contributions of Equation (14)), $Ker{r}_{221}$ ($\ell =\left|m\right|=2$, $n=0,1$ contributions of Equation (14)), and $Ker{r}_{HM}$ (all fundamental prograde modes with $\ell \le 4,n=0$ contributions of Equation (14) and also taking into account mode-mixing [172]). The frequencies and damping times are predicted in terms of final mass and spin for all these cases. The remnant quantities, $M_f$ and ${\chi }_f$, are estimated assuming uniform priors on these parameters. We do not consider the higher overtones (n > 1) as those are not expected to provide constraints with the current sensitivity of detectors.

Another equally established technique for the QNM analysis is the parametrized-SEOBNRv4HM (pSEOB) analysis, which employs a parametrized version of the EOB waveform model [174] accounting for aligned spins and higher modes [30,166,175]. The pSEOB analysis differs from PYRING in that it measures the ringdown frequency within a complete IMR waveform model framework using the full SNR of the signal. It is not dependent on a ringdown time definition. In this framework, one parameterizes the frequency and damping time of the $\ell =m=2$ by introducing a fractional deviation from the nominal GR prediction and constrains these fractional deviation parameters directly from the data. That is, $f_{220}=f_{220}^{GR}(1+\delta {\widehat{f}}_{220})$ and ${\tau }_{220}={\tau }_{220}^{GR}(1+\delta {\widehat{\tau }}_{220})$ [174], where $f_{220}^{GR}$ and ${\tau }_{220}^{GR}$ are frequency and damping time if GR is the correct theory of gravity. From both these approaches, by performing detailed analyses on the GWTC-2 events, no indication of the presence of non-BH behavior was reported [30].

2.3.2. Testing the Nature of Merger Remnant from the Measurement of Late Ringdown Echoes

One can ask the question as to if the merger remnant is not a BH but instead an exotic compact object (ECO) with a light-ring and reflective surface, instead of an event horizon as in the case of BHs [176,177,178,179]. For these hypothetical cases, the GWs can be trapped in between the effective potential at the centre and the reflective surface, leading to the emission of GWs as a train of repeating pulses known as GW echoes. BHs produce no echo signal as there is no possibility of an out-going boundary condition at the BH event horizon. GW echoes are unique probes of any non-BH compact object formation (especially exotic objects like gravastars, fuss balls) after the binary coalescence [47,176,180,181]. For an illustration of GW echoes originating from a binary merger, we point out Figure 7 (also see Figure 2 of Reference [176] for the original GW150914 signal on top of the best-fit echo template from a template-based echo analysis).

In the template-based framework, the echo signal is modeled with five extra parameters, characterizing the echo: the relative amplitude of the echoes, the damping factor between each echo, the start time of ringdown, the time of the first echo concerning the merger, and the time delay between each echo ($\Delta t_{echo}$ in Figure 7). Reference [176] studies this method and subsequent discussions on GW150914 in detail. For the GWTC-2 events [30], assuming uniform priors on each of these echo parameters and employing the IMRPhenomPv2 waveform model (except for the case GW190521, where Sur7dq4, is a surrogate model for precessing BBH system directly interpolates the numerical relativity waveforms, is used), a Bayesian analysis is performed to investigate the evidence for echoes. Bayes factor $B_{IMR}^{IMRE}$ (comparing the two hypotheses $IMRE$ data best fits an echo model, and data best fits a model without echoes $\left(IMR\right)$) are computed for each event. The data did not show evidence for echoes [30], except for the event GW190915_235702, which showed the highest value for $B_{IMR}^{IMRE}=0.17$ indicating a negligible evidence for echoes [30]. Reference [30] reports that the posteriors on the echoes parameters returned their prior distribution, pointing to a null detection of GW echoes.

2.4. Constraints on the Polarization States of GWs

Generic metric theories of gravity predict six independent degrees of freedom for the metric tensor, which can be identified as polarization states of GWs. More than the two tensor (spin 2) degrees of freedom allowed by GR, there is a possibility of two vectors (spin 1) degrees of freedom and two scalars (spin 0) degrees of freedom [141] in such cases. Out of these six modes, three of them are transverse, and three are longitudinal. The first constraints on the polarization states of the GWs from observations are detailed in [9]. Still, the results were uninformative as the data from the two LIGO instruments are not enough to constrain the extra polarization mode. This is possible only if there are data available from one another detector with another orientation. Hence the polarization test was first demonstrated with the three-detector event GW170814, and the improved results compared to Reference [9] are available in Reference [22]. In Reference [22] also the tensor modes hypothesis was favored over scalar and vector modes as shown in Reference [9].

A detailed analysis was performed in Reference [30] considering all GW events observed till the first half of the third observing run of LIGO/Virgo detectors. Reference [30] reports the highest (lowest) Bayes factor for GW190720_000836 (GW190503_185404) with ${log}_{10}{B}_V^T=0.139$, $B_S^T=0.138)$ (${log}_{10}{B}_V^T=0.074$ and $B_S^T=-0.072)$, here $B_V^T$ and $B_S^T$ represent the Bayes factors for full tensor versus full-vector and full-scalar hypotheses respectively. The $B_S^T$ is slightly larger than $B_V^T$ and this is explained by the intrinsic geometry of the LIGO-Virgo antenna patterns [182]. One should also notice that any of these results do not account for the possibility of mixed polarization. This topic has to be explored in the future when more detectors become operational.

3. Summary

This review article provides a brief overview of the tests of GR performed during the first three observing runs of the LIGO-Virgo detectors, including tests of consistency with GR (Section 2.1), parameterized tests (Section 2.2), tests based on the merger remnant properties (Section 2.3), and tests for GW polarizations (Section 2.4). Along with some technical details about each test, we also provide a short discussion pointing to the prospects of these various tests.

In this article, we only focused on signals which are consistent with BBHs. The detection of the first binary neutron star merger event opened up different possibilities of testing GR from combined electromagnetic and GW observations [23,118]. Many tests we detailed here may have overlaps or redundant information which is not accounted for here. Though the instrumental noise mainly dominates the current measurement uncertainties, we cannot exclude the possibility of any systematic bias arising due to un-modeled effects present in the waveform models. Another critical point is that the tests discussed here are all model-independent tests (null tests). In other words, we are not assuming any alternative gravity theory models here, and every test is capturing the deviation from GR in a model-agnostic way.