Tests of Gravitational-Wave Birefringence with the Open Gravitational-Wave Catalog

We report the results of testing gravitational-wave birefringence using the largest population of gravitational-wave events currently available. Gravitational-wave birefringence, which can arise from the effective field theory extension of general relativity, occurs when the parity symmetry is broken, causing the left- and right-handed polarizations to propagate following different equations of motion. We perform Bayesian inference on the 94 events reported by the 4th-Open Gravitational-wave Catalog (4-OGC) using a parity-violating waveform. We find no evidence for a violation of general relativity in the vast majority of events. However, the most massive event, GW190521, and the second most massive event, GW191109, show intriguing non-zero results for gravitational-wave birefringence. We find that the probability of association between GW190521 and the possible electromagnetic (EM) counterpart reported by Zwicky Transient Facility (ZTF) is increased when assuming birefringence. Excluding GW190521 and GW191109, the parity-violating energy scale is constrained to $M_\mathrm{PV}>0.05$ GeV at $90\%$ credible interval, which is an improvement over previous results from twelve events by a factor of five. We discuss the implications of our results on modified gravity and possible alternative explanations such as waveform systematics. More detections of massive binary black hole mergers from the upcoming LIGO/Virgo/KAGRA run will shed light on the true origin of the apparent birefringence.

This work tests GW birefringence using the currently largest population of gravitational-wave events from 4-OGC. Birefringence of GWs emerges when the parity symmetry, which is the invariance of physical laws regarding the inversion of spatial coordinates, is bro-ken between the left-and right-handed GW polarizations. While parity symmetry is conserved in GR, theories where parity is violated have been proposed such as Chern-Simons gravity [23][24][25][26][27], Hořava-Lifshitz gravity [28][29][30], ghost-free scalar-tensor gravity [31], and the symmetric teleparallel equivalent of GR [32] to account for dark matter and dark energy. Parity violation also arises at high energy scales in quantum gravity theories such as loop quantum gravity and string theory [23].
We utilize an effective field theory (EFT) extension of the linearized Einstein-Hilbert action to study how deviations from GR affect GW propagation. EFT is a flexible framework that includes all action terms that purposely preserve or violate certain symmetries. The leading order higher derivative modification of the linearized action comes from terms of mass dimension five that violate parity [33,34]. Parity violation leads to an asymmetry in the propagation speeds and amplitudes of the left-and right-hand polarization of GW, which, in turn, leads to phase and amplitude birefringence, respectively. Given the relationship between parity violation and Lorentz violation [35], our tests have implications for constraining the standard model extension (SME) of gravity [36,37], which is the most general EFT extension of linearized GR that violates Lorentz symmetry.
We test for GW birefringence in ninety-four GW events and find no evidence of birefringence in ninetytwo events. Intriguingly, we find two outliers, GW190521 and GW191109 010717 (hereafter GW191109), the most and second most massive binary black holes in 4-OGC [11,49,50], favor the birefringence waveform over the GR waveform with Bayes factors 11.0 and 22.9, respectively (ln B nonGR GR = 2.4 and 3.1). Excluding these two outliers, we place constraints on the cutoff energy scale M PV > 0.05 GeV, which is more stringent than Ref. [22] by a factor of five due to the larger data set and the more advanced waveform. The result can be mapped to the SME coefficient ς (5) < 9 × 10 −16 m, which characterizes isotropic birefringence with mass dimension d = 5.
The nonzero results may be caused by non-Gaussian noise fluctuation, systematic errors in waveform templates, or physical effects that are outside the standard assumptions for quasi-circular binary black hole mergers, e.g., the presence of orbital eccentricity [51,52], a hyperbolic encounter [53], or even entirely new physics [54]. We assess how the background noise affects the M pv measurement for both GW190521 and GW191109 by injecting 100 GR signals into the detector noise around the observed coalescence time for each event. The results suggest that our detection of GR deviation for GW190521 and GW191109 is not likely an artifact of the noise. Lastly, we consider the possible EM counterpart of GW190521, which ZTF observed and comes from an active galactic nucleus [55,56], by comparing the sky location posteriors for GW190521 with and without birefringence. Prior analyses of GW190521 do not strongly favor association with the EM counterpart. Reference [57] reported the Bayes factor of ln B −4 which disfavors a sky location fixed to that of the electromagnetic counterpart. However, if we assume GW birefringence, we find evidence in favor of the association with the EM counterpart with ln B of 6.6.

II. WAVEFORM TEMPLATES FOR GW BIREFRINGENCE
We briefly overview the construction of waveform templates to test GW birefringence following Ref. [34]. From the perspective of EFT, the leading order modifications to the linearized Einstein-Hilbert action are terms with three derivatives and mass dimension five: ijkḣ il ∂ jḣkl and ijk ∂ 2 h il ∂ j h kl [33,34,58], where i, j... = 1, 2, 3 refer to spatial coordinates, ∂ j denotes spatial derivatives, dot denotes derivatives with respect to time, ∂ 2 is the Laplacian, ijk is the antisymmetric symbol, and h ij is the tensor perturbation of metric. Notably, both terms violate parity. Therefore, EFT suggests that the leading order modification to GW propagation arises from parity-violation. We do not consider the more complicated anisotropic GW birefringence, for which the effects can be found in Refs. [38,48,59]. Combining the above higher derivative terms with the linearized Einstein-Hilbert action gives where a is the cosmic scale factor, M PV is the energy scale at which higher order modification starts to be relevant, and c 1 and c 2 are two undetermined functions of cosmic time; the speed of light and the reduced Planck's constant are set to c = = 1. Equation (1) is the generic form of the action; c 1 , c 2 and M PV can be mapped to the corresponding model parameters in a specific modified gravity theories [23][24][25][26][27][28][29][30][31][32], as explicitly demonstrated in ref. [34]. The equation of motion for the GW circular polarization modes h A , where A = R or L for the right-and left-hand modes, is where H is the conformal Hubble parameter, k is the wavenumber, and a prime denotes the derivative with respect to the cosmic conformal time τ , µ A and ν A are the phase and amplitude birefringence parameters. They have the exact forms ρ A = ±1 for left-and right-handed polarizations represents broken parity of GWs during propagation. We focus on velocity birefringence because the modification to GW strain from amplitude birefringence is negligible [34]. The GR solution can be found by setting Solving Eq. (2) gives the left-and right-handed circular polarization modes with parity violation. They are related to the GR waveform by The plus (h + ) and cross (h × ) modes of GW are given by . For binary black holes, we use the IMRPhenomX-PHM [60] GR waveform approximant, which includes sub-dominant harmonics and effects of precession; the IMRPhenomD NRTidal [61][62][63][64] and IMRPhenomNSBH [65], waveforms, which account for tidal deformability, are used for binary neutron stars and neutron-star-blackhole mergers respectively.
The phase modification to the GR waveform takes the following form where H 0 is the Hubble constant, z is the cosmic redshift, the frequency term f 2 corresponds to a modification at 5.5 post-Newtonian order, Ω M is the matter density, and Ω Λ is the dark energy density. We adopt a ΛCDM Cosmology with parameters Ω M = 0.3075, Ω Λ = 0.691, and H 0 = 67.8 km s −1 Mpc −1 [66] to convert luminosity distance to redshift for Eq. (5). As most GW detections are from the local Universe, we make the simplifying assumption that c 1 , c 2 are constants and c 1 − c 2 is of order unity. This is done by attributing its order of magnitude to M PV . Also note that we do not consider the special case where c 1 = c 2 exactly, and thus µ A = 0, in this work. This is the case for dynamical Chern-Simons gravity [67][68][69] and the constraints on amplitude birefringence in this case can be found in Refs. [41,45].

III. BAYESIAN INFERENCE
We use Bayesian parameter estimation and model selection to test GW birefringence. Given data d(t), which is a sum of the detector noise n(t) and a possible GW signal h(t, θ) characterized by parameters θ, Bayes theorem states that where P ( θ|d, H) is the posterior probability distribution for parameters θ, P ( θ|H) is the prior distribution containing any a priori information, P (d| θ, H) is the likelihood for obtaining the data given model parameters, P (d|H) is a normalization factor called evidence, and H is the hypothesis for modeling the data. This work considers two competing hypotheses: H GR where GWs are accurately described by GR and H nonGR where GWs have birefringence and are described by Eq. (4). The Bayes factor, or the ratio of the evidences of two hypotheses, quantifies the degree that data favor one hypothesis over another.
Using PyCBC Inference [70], we numerically sample over all GW intrinsic (mass m 1,2 , spin s 1,2 , and, in the case of neutron stars, tidal deformability Λ 1,2 ) and extrinsic parameters (luminosity distance d L , inclination angle ι, polarization angle Ψ, right ascension α, declination δ, coalescence time t c , and phase φ c ) as well as the parity violation parameter M −1 PV for H 1 . The priors for the intrinsic and extrinsic parameters are the same as those in the 4-OGC [10] which are uniform for mass, spin amplitude, polarization, coalescence phase, and time. The distance prior is uniform in comoving volume. The spin orientation, sky location, and spatial orientation priors are isotropic. The prior for M −1 PV is uniform in [0,200] GeV −1 except for the two outlier events which use [0,1000].
Assuming M −1 PV is a universal quantity, the M −1 PV posteriors from the individual GW events can be combined to obtain an overall constraint, where d i denotes data of the i-th GW event.

IV. RESULTS
We find that ninety-two out of ninety-four events are consistent with GR; the M −1 PV posteriors are shown in Fig. 1, with the most constraining events individually indicated in the legend. In general, the tightest constraints are all from events with relatively low mass (∼ 10 M ) binary black hole mergers e.g., GW190707 093326 has component masses of 12.9M and 7.7M [10]. This result is unsurprising because the birefringence effect is more significant at higher frequencies, where low mass binaries spend more time (see Eq. (5)).
The overall constraint is obtained by multiplying the posterior distributions of the individual events together using Eq. (8). We find the 90% upper limit of M −1 PV to be 19 GeV −1 , which corresponds to M PV > 0.05 GeV. This result is more stringent than previous results, based on twelve GW evens, by a factor of five [22]. Note that Ref. [22] used h = 1 not = 1. Taking this into account, their result is M PV > 0.01 GeV. This improvement is due to the increased number of events analyzed and the use of improved GW waveforms [60] with longer duration and higher harmonic modes.

A. GW190521 and GW191109
We find non-zero results for birefringence in GW190521 and GW191109 with M −1 PV = 400 +460 −230 GeV −1 and 220 +150 −100 GeV −1 (90% credible interval). To investigate any systematics causing the apparent deviation from GR, we perform two more birefringence  tests for the outlier events using a different phenomenological template model IMRPhenomPv3HM [72] and a numerical relativity surrogate model NRSur7dq4 [73]. We find consistent support for non-zero M −1 PV in the posteriors for these runs. We further check the data quality around GW190521 and GW191109 by calculating the background noise power spectral density (PSD) variation (see [74] for definition), which measures the noise non-stationarity. For GW190521, the PSD variation in a one-hour interval around the event deviates from Gaussian stationary noise by 0.1 (except for a glitch in LIGO Hanford 400 s after the event), showing no significant deviation from other ordinary times. However, the LIGO detector data for GW191109 contains non-Gaussian and non-stationary transient noise artifacts or glitches [6]. Due to this, we performed our analysis using data with the glitch removed released by LIGO/Virgo [6].
To further quantify whether detector noise could be the cause of the observed non-zero result, we simulate 100 GR signals with parameters drawn from the GW190521 and GW191109 posteriors [11] and inject them into the LIGO/Virgo detector noise nearby the GW190521 and GW191109 triggers, respectively (see the appendix for technical details). We find only 1(4) events, out of 100 injections, have ln B nonGR GR larger than what we found for GW190521(GW191109). We thus conclude that the false alarm rate for detecting birefringence is 1(4) in 100 observations.
Lastly, we consider the possible EM counterpart for GW190521 reported by the ZTF [55,56]. Interestingly, we find that including birefringence significantly improves the chance of association. Fig. 3 shows the right ascension (α), declination (δ), and luminosity distance posteriors from the birefringence analysis. The red lines mark the independent measurements by ZTF (α = 3.36 rad, δ = 0.61 rad, and redshift=0.438). The Bayes factor supports coincidence with ln B overlap = 6.6 in favor of the association. The GR analyses did not favor association. For instance, Ref. [57] reports ln B overlap ∼ −4.4.

V. DISCUSSION AND CONCLUSION
We test for GW propagation birefringence using stateof-the-art waveform templates and 4-OGC, the most extensive GW catalog currently available. Combining the results from all of the events except the outliers GW190521 and GW191109, we constrain the lower energy scale cutoff to M PV > 0.05 GeV, which is an improvement over previous constraints by a factor of 5. The constraint on M PV allows us to limit the SME isotropic birefringence parameter with mass dimension d = 5 to ς (5) < 9 × 10 −16 m. These results show that GW astronomy is a promising future avenue by which to study gravity at high energies.
We surprisingly find evidence in support of GW birefringence for GW190521 and GW191109, which happen to be the most and second most massive events in 4-OGC. Furthermore, we find that including birefringences increases the likelihood of association between GW190521 and its possible EM counterpart.
We find no disparity between three state-of-the-art waveform approximants and no significant issues with the data quality for GW190521. While there is a glitch in the GW191109 data, it was removed by LIGO/ Virgo, and our analysis suggests the noise fluctuation is unlikely to have caused the non-zero M −1 PV result. However, it is well documented that GW190521 is an exceptional event that may not fit well into the simple quasi-circular binary black hole merger picture. For instance, Refs. [51,52] show that GW190521 is consistent with the merger of a binary black hole system with eccentric orbit, Ref. [53] gives a Bayes factor that highly favors a hyperbolic encounter over a quasi-circular merger, and Ref. [54] shows that GW190521 could be genuinely new physics, such as a Proca star collision. The accuracy of current GR waveform approximants is limited at the merger stage. This is quite relevant for GW190521 and GW110919 as most of the data is in the merger band. Our work provides further evidence for non-standard physical effects in GW data, which the available GR waveform approximants cannot currently account for. Even if the apparent deviation from GR is from new physics, our GW birefringence model can not provide a universal explanation. A possible extension might include a parity violation that depends on the masses of the binary.
As the advanced LIGO and Virgo detectors are upgraded and the KAGRA detector joins the network, we expect more high mass detections similar to GW190521 and GW191109, which may provide further insight into the physics behind the observed behavior of these outliers.
We release all posterior files and the scripts necessary to reproduce this work in https://github.com/ gwastro/4ogc-birefringence. We investigate the data quality to determine if nonstationary or non-Gaussian noise could be responsible for the apparent deviation from GR for GW190521 and GW191109. We generate 100 GR waveforms with the IMRPhenomXPHM [60] template and inject them into the LIGO and Virgo data near the coalescence times of GW190521 and GW191109. The waveform source parameters (component masses and spins, sky location, source orbital orientation, GW polarization angle, and coalescence phase) are sampled from the GW190521 and GW191109 posteriors released by 4-OGC [11]. For GW190521, the waveforms are injected into the time interval [-20, 20] seconds around the trigger time of GW190521. However, we exclude the region [-4,4] seconds around the trigger time, where the event is predominant over the noise. Both LIGO detectors contain transient noise artifacts at the merger time of GW191109 [6]. We have analyzed GW191109 using a deglitched version of data released by LIGO and Virgo [6]. However, a cleaned version of the data in the region of GW191109 is not readily available. To avoid the glitches, which were removed for the GW191109 analysis, we inject our signals into a segment of data [-100, -30] seconds from the trigger.
The injections are then analyzed the same as the actual event using the method presented in the main text. For the parameter estimation, all priors are the same as those used to analyze the GW190911 and GW190521 respectively.
The injections' M −1 PV posteriors are shown in the first row of Fig. 4. For comparison, the figure also shows the GW190521 and GW191109 posteriors. We note that some results for GW191109-type injections have a spikey or multi-modal shape; we attribute this to the noisy data around GW191109.
We find that, in rare cases, the background noise fluctuation can produce non-zero peaks for M −1 PV . To assess the significance of our non-zero M −1 PV results for GW191109 and GW190521, we extract the Bayes factor from the posteriors using the Savage-Dickey density ratio; the second row of Fig. 4 shows the results. Using ln B nonGR GR as a statistic, only 1 in 100 and 4 in 100 simulations exceed the Bayes factor from GW190521 and GW191109, respectively. The false alarm rate or p-value to reject the null hypothesis that GR is correct is thus 1(4) out of 100 realizations.
We take the mean of M −1 PV of all injections to determine the background distribution for the null hypothesis and note that the background is qualitatively different from the actual data results. Overall we do not find strong evidence that background noise artifacts caused the nonzero results for testing GR.