Eccentricity-induced systematic error on parametrized tests of general relativity: Hierarchical Bayesian inference applied to a binary black hole population

One approach to testing general relativity (GR) introduces free parameters in the post-Newtonian (PN) expansion of the gravitational-wave (GW) phase. If systematic errors on these testing GR (TGR) parameters exceed the statistical errors, this may signal a false violation of GR. Here, we consider systematic errors produced by unmodeled binary eccentricity. Since the eccentricity of GW events in ground-based detectors is expected to be small or negligible, the use of quasicircular waveform models for testing GR may be safe when analyzing a small number of events. However, as the catalog size of GW detections increases, more stringent bounds on GR deviations can be placed by combining information from multiple events. In that case, even small systematic biases may become significant. We apply the approach of hierarchical Bayesian inference to model the posterior probability distributions of the TGR parameters inferred from a population of eccentric binary black holes (BBHs). We assume each TGR parameter value varies across the BBH population according to a Gaussian distribution. The means and standard deviations that parametrize these Gaussians are related to the statistical and systematic errors measured for each BBH event. We compute the posterior distributions for these Gaussian hyperparameters, characterizing the BBH population via three possible eccentricity distributions. This is done for LIGO and Cosmic Explorer (CE, a proposed third-generation detector). We find that systematic biases from unmodeled eccentricity can signal false GR violations for both detectors when considering constraints set by a catalog of events. We also compute the projected bounds on the 10 TGR parameters when eccentricity is included as a parameter in the waveform model. This is done via multiplying the individual likelihoods for each event in the catalog, and also by combining them via hierarchical inference. We find that the first four dimensionless TGR deformation parameters can be bounded at 90% confidence to $\delta {\widehat{\phi }}_i\lesssim {10}^{-2}$ for LIGO and $\lesssim {10}^{-3}$ for CE [where $i=(0,1,2,3)$; GR predicts zero for all values of the $\delta {\widehat{\phi }}_i$]. The most stringent bound applies to the $-1\mathrm{PN}$ (dipole) parameter: it is constrained to $|\delta {\widehat{\phi }}_{-2}|\lesssim {10}^{-5}$ ($\lesssim 4\times {10}^{-7}$) by LIGO (CE). In comparison to the circular orbit case, the combined bounds on the TGR parameters worsen by a modest factor of $\lesssim 2$ when eccentricity is included in the waveform. (The dipole parameter bound degrades by a factor of $\sim 3–4$ in this case).

Figure 1

Schematic diagram depicting the expected behavior of the population hyperparameters ${\mu }_i$ and ${\sigma }_i$ for the unbiased (upper panels) and biased (lower panels) cases. The blue curves represent the one-dimensional marginalized posterior probability distributions for the hyperparameters ${\mu }_i$ and ${\sigma }_i$. For simplicity, we represent the posteriors via Gaussian distributions. In the absence of any bias (upper panel), the peaks of the probability distributions are exactly at the GR values (zeros) for both ${\mu }_i$ and ${\sigma }_i$. In this case, the width of the distributions depends only on the statistical width $({\tilde{\sigma }}_i^j)$ of the individual $\delta {\widehat{\phi }}_i$ posteriors [$p(\delta {\widehat{\phi }}_i^j)$] in the population. The use of an inaccurate waveform model shifts the peak of the individual posteriors $p(\delta {\widehat{\phi }}_i^j)$ away from zero. The shift in the ${\mu }_i$ posterior is given by the mean of the systematic bias (${\mu }_{\mathrm{sys}}$) in the $\delta {\widehat{\phi }}_i^j$ posteriors (assuming the statistical errors on the $\delta {\widehat{\phi }}_i^j$ are the same for all events, ${\tilde{\sigma }}_i^j={\tilde{\sigma }}_{0,i}$). The shift in the ${\sigma }_i^2$ posterior is proportional to the variance in the systematic biases (${\sigma }_{\mathrm{sys}}^2$). The widths of the ${\mu }_i$ and ${\sigma }_i$ posteriors in the biased case depend on ${\sigma }_{\mathrm{sys}}$ and ${\tilde{\sigma }}_{0,i}$ [see Eqs. (8) and (9)].

Figure 2

Eccentricity distributions for the BBH populations considered in this study. The Powerlaw distribution assumes a PDF with distribution $e_0^{-\alpha }$ with power law index $\alpha =1.2$. The Loguniform ($\alpha =-1$) distribution has the largest number of sources with higher eccentricities $e_0>{10}^{-3}$ (measurable in CE), followed by the Zevin [124] and Powerlaw distributions.

Figure 3

Posteriors on the hyperparameters ${\mu }_i$ (left panels) and ${\sigma }_i$ (right panels) for the leading-order TGR parameters $\delta {\widehat{\phi }}_i$ as measured with a single LIGO detector. Panels refer to the 0PN, 1PN, 1.5PN, and 2PN deviation parameters (top to bottom). The three colors represent the three eccentricity distributions considered in Fig. 2. The vertical black line indicates the GR value (zero). The dashed (colored) lines represent the 90% credible interval for each posterior. Most posteriors exclude the GR value, indicating a (false) deviation from GR. For each eccentricity distribution we combine data from 303 sources with $\mathrm{SNR}\ge 12$ and $M\le 60M_{\odot }$.

Figure 4

Posteriors on hyperparameters ${\mu }_i$ and ${\sigma }_i$ for the higher-order TGR parameters $\delta {\widehat{\phi }}_i$ (2.5PN log, 3PN, 3PN log, 3.5PN) as measured by LIGO. The bottom row shows the -1PN order (dipole) parameter. Conventions are as in Fig. 3.

Figure 5

Posteriors on the hyperparameters ${\mu }_i$ (left panels) and ${\sigma }_i$ (right panels) for the leading-order $\delta {\widehat{\phi }}_i$ parameters as measured with a single CE detector. Conventions are as in Fig. 3.

Figure 6

Posteriors on hyperparameters ${\mu }_i$ and ${\sigma }_i$ for the higher-order and -1PN order TGR parameters $\delta {\widehat{\phi }}_i$ as measured by a single CE detector. Conventions are as in Fig. 4.

Figure 7

90% upper bound on the testing GR parameters $\delta {\widehat{\phi }}_i$ from 303 BBHs in our simulated population. The bounds are obtained by assuming that the $\delta {\widehat{\phi }}_i$ take the GR values (zero) for all events. The green triangles represent the projected bounds for LIGO at design sensitivity, while the orange circles show the projected bounds for the 3G CE detector. The filled symbols correspond to bounds obtained with eccentric waveforms (including $e_0$ as a parameter). Unfilled symbols show bounds in the circular case. The eccentricity samples for ${10}^5$ BBHs are drawn from the Loguniform distribution with $e_0\in ({10}^{-7},0.2)$. The blue diamonds represent the 90% upper bound from the GWTC-3 [79] analysis.

Figure 8

Joint probability distributions for the hyperparameters ${\mu }_i$ and ${\sigma }_i$ that characterize the 10 TGR parameters. The contours represent the 90% credible region set by measurements with a LIGO-like detector operating at design sensitivity. We combine the 303 BBH events as discussed in Sec. 5. The solid contours consider the case of eccentric binaries; dashed contours represent the circular case.

Figure 9

Joint probability distributions for the ${\mu }_i$ and ${\sigma }_i$ hyperparameters in the case of a single CE detector. Conventions are the same as Fig. 8.