Tests of gravitational wave propagation with LIGO-Virgo catalog

The detection of gravitational waves (GWs) by LIGO [1] and Virgo [2] collaborations heralds a new era in both astrophysics and fundamental physics [4, 5, 3, 6, 7, 8]. During the third observing run (O3), Advanced LIGO, Advanced Virgo and KAGRA added 79 GW events, as detailed in GWTC-2.1 [7] and GWTC-3 [8]. These discoveries deepen our understanding of the universe and provide powerful tools to check the rules of gravitational theories such as general relativity (GR), and help us better understand the basic structure of spacetime.

According to GR, GWs are theorized to propagate at the speed of light in a vacuum, and their velocity $v_{g}$ is predicted to be invariant across all frequencies. This characteristic suggests that the propagation speed of GWs is unaffected by frequency, ensuring uniformity in the velocity across different wavelengths. The GWs detected by Advanced LIGO and Advanced Virgo can be used to make statistical inferences about $v_{g}$[9, 10, 11, 12], and the speed of GWs can be constrained within the range of $0.99^{+0.01}_{-0.02}c$ at a 90% credible interval[11]. In the framework of the quantum gravity theory, gravitons are postulated to be carriers for transmitting gravitational interactions. GWs travel at the speed of light across all frequencies which supports the idea that gravitons are massless.

However, certain modified gravitational theories typically anticipate violations of Lorentz invariance (LI) [13, 14, 15, 16, 17, 18, 19, 20], indicating that the speed of GWs might not always align with the speed of light and could vary based on their frequencies during propagation [21, 22, 23, 24]. These theories imply that gravitons may possess a rest mass that is not equal to zero. We use modified dispersion relations to depict these theories[26],and we focus our discussion on a specific scenario, which is the massive gravity[25]. It is hypothesized that low-frequency GWs propagate slower than high-frequency GWs. Consequently, during the inspiral phase of compact binary stars, gravitons emitted earlier have a lower propagation speed than those emitted later. This phenomenon induces phase alterations in GWs during propagation, resulting in received waveforms that differ from those predicted by GR. If such dephasing of waveforms comparing with GR’s templates is not detected, this could imply a constraint on the graviton mass, thereby limiting the potentially extent of LI violation. The graviton mass has been measured through various experimental datasets[27, 28, 29, 30, 31, 32].

Similar to the impact of the massive gravity, if the weak equivalence principle (WEP) is violated, the speed of GWs will no longer be the speed of light but will depend on the frequencies. This frequency dependence allows us to limit the extent of WEP violation. The violation of WEP results in different Shapiro delays [33] for different massless particles traversing an identical gravitational field. The Shapiro delay is proportional to $\gamma+1$, where $\gamma$ is the parametrized post-Newtonian parameter ($\gamma=1$ in general relativity) [34, 35], and different massless particles will have different values of $\gamma$ during free fall in a gravitational field if WEP is violated. It is often considered that two different massless particles are emitted from the same astrophysical source, with respective values of $\gamma_{1}$ and $\gamma_{2}$, then $|\Delta\gamma|$ ($|\gamma_{1}-\gamma_{2}|$) can be used to quantitatively represent the violation of the WEP. If we know the intrinsic time delay of these two particles, we can limit the difference in the Shapiro time delay between the two particles $\Delta t_{\rm gra}$ [36] which is proportional to $|\Delta\gamma|$, by measuring the observed time-delays $\Delta t_{\rm obs}$. Therefore, we can establish an upper limit for $|\Delta\gamma|$ by observing the $\Delta t_{\rm obs}$ of various signals emitted from the same source. Several constraints on $|\Delta\gamma|$ have been derived from various astrophysical events, including the emissions from the supernova event SN1987A [36, 37], gammy-ray bursts [38, 40, 39, 41, 42, 43] and fast radio bursts (FRBs) [44, 45, 46, 47]. In this study, we use GWs to test WEP, which had previously been employed [48, 49, 50, 51]. The utilization of GW signals can eliminate the intrinsic time delay when signals are emitted, which will enable us to estimate the up-limit of $|\Delta\gamma|$ more accurately.

Building on the previously discussed framework, we further explore the phenomena of massive gravity and WEP violation, both of them modify the speed of GWs resulting in frequency dependence. Specifically, if gravitons possess mass, the speed of gravitons emitted during the early inspiral phase of GW emission is expected to be lower than that emitted later. Consequently, the impact of this effect on the GW waveform is expected to become more pronounced as the luminosity distance of GW events from Earth increases. This is consistent with the behavior under WEP violation, where the influence of gravitational fields leads to cumulative effects as the luminosity distance of the GW increases. However, it is essential to note that WEP violation manifests through changes in the speed of GWs due to the influence of gravitational fields.

In this study, our analysis focuses on the gravitational potential generated by the Milky Way, under the assumption that its entire mass is concentrated at the center. Owing to the unclear details of the host galaxies and galaxies along the propagation paths of these events, we neglect the time delays caused by these galaxies. This approach renders the constraints in our study very conservative, as the impacts from the host galaxy and other galaxies may possess similar or even greater magnitude compared to those from the Milky way. Therefore, the actual values of $|\Delta\gamma|$ might be several times smaller than our results.

The remainder of this paper is organized as follows. Section 2 introduces our methodology for analyzing the specific dispersion of GWs. In Section 3, we present our results for 22 GW events. Conclusions and discussions are presented in Section 4.

We use the common phenomenological modification to GR introduced in [26], which was previously applied to LIGO and Virgo data in [52, 28, 29, 30]:

where c is the speed of light, E and p are the energy and momentum of the GWs, respectively, $A_{\alpha}$ and $\alpha$ are phenomenological parameters. With different $A_{\alpha}$ and $\alpha$, we can use Eq. (1) to denote different modified dispersion relations in different alternative gravitational theories. In this study, we only consider the scenarios where $\alpha=0$ and $A_{0}>0$, for which the graviton mass is given by $m_{g}=A_{0}^{1/2}/c^{2}$, following the framework of massive gravity[25].

The modified gravitational theory will lead to a difference in the dispersion relation of GWs, which will lead to the propagation of GWs at a speed different from that of light. Therefore, a modified waveform template that includes such effects is generally used to describe the waveform of GWs, which is expressed in the frequency domain as follows

where $\tilde{A}(f)$ denotes the complex amplitude, $\Psi_{\mathrm{GR}}(f)$ denotes the complex phase predicted by GR, and $\delta\Psi(f)$ is the modification term. In the scenarios where $\alpha=0$ and $A_{0}>0$, it is given by

where $z$ is the cosmological redshift, $h$ is the Planck constant, $f$ is the frequency of the GWs, and $D_{0}$ is the distance parameter given by

where $H_{0}=67.90\text{ km s}^{-1}\text{ Mpc}^{-1}$ is the Hubble constant, $\Omega_{\text{m}}=0.3065$ and $\Omega_{\Lambda}=0.6935$ are the matter and dark energy density parameters, respectively. These are the TT+lowP+lensing+ext values from [56].

The Shapiro time delays generated during the propagation of GWs will also induce frequency-dependent variations in their speed owing to the violation of WEP. Here, we assume that $\delta\Psi(f)$ is caused by WEP deviation. Considering a GW event, GWs emitted at $t_{\rm{e}}$ and $t^{\prime}_{\rm{e}}$ with different frequencies are received by the detectors at the corresponding arrival times $t_{\rm{a}}$ and $t^{\prime}_{\rm{a}}$. We can ignore the cosmological inflation effect if the difference in emitting time ($\Delta t_{\rm{e}}=t_{\rm{e}}-t^{\prime}_{\rm{e}}$) is so small that the difference in the arrival times of the corresponding GWs ($\Delta t_{\rm{a}}=t_{\rm{a}}-t^{\prime}_{\rm{a}}$) is

In this study, we only consider the violation of the WEP caused by the Milky Way, and $\Delta t_{\rm gra}$ would be [49]

where the Milky Way mass $M_{\mathrm{MW}}\approx 6\times 10^{11}~{}M_{\odot}$ and the distance from the sun to the center of the Milky Way $r_{G}\approx 8~{}\mathrm{kpc}$, $d$ denotes the distance from the GW event to the Milky Way center, $b$ represents the impact parameter of the GW paths relative to the center of the Milky Way. We use a transform formula [57] to convert the celestial coordinates to $b$. $s_{n}=+1$ denotes the GW event positioned in the direction of the Milky Way center, and $s_{n}=-1$ denotes the GW event positioned away from the Milky Way center. However, Minazzoli et al. have raised concerns about the application of Eq. (6) for calculating Shapiro delay on cosmological scales [58], as these concerns extend beyond the current study’s scope, we exclude them from our analysis. Then, we obtain $\delta\Psi(f)$

where $\Delta f=f-f^{\prime}$ and $f,~{}f^{\prime}$ are two different frequencies of GWs in a single event. It is intuitive to assume that $\gamma$ is proportional to the particle energy. Considering $E=hf$, we assume that $\Delta\gamma\propto\Delta f$.

In this study, we employ the BILBY Bayesian parameter estimation software [54] and the DYNESTY nested sampling package [55] to perform parameter estimation. We use the parametrized waveform model described in Eq. (2), which incorporates modification terms into the IMRPhenomXP waveforms [53]. These modifications are designed to estimate parameters $m_{g}$ and $\Delta\gamma$ using the selected GW events in GWTC-2.1 and GWTC-3. The GW strain data are obtained from the Gravitational Wave Open Science Center (GWOSC), which is publicly accessible at [https://gwosc.org/]. For the parameter priors, we adopt BILBY’s default settings for binary black holes (BBHs) [54] except for $m_{g}$ and $\Delta\gamma$. For $\Delta\gamma$, we select a modified logarithmic prior that encompasses both negative and non-negative values (refer to Eqs.(6) in [51]). Furthermore, inspired by this approach, we similarly define the prior distribution for $m_{g}$ as

where $\alpha$ is a uniform distribution parameter. Based on previous studies, we set the prior range for $m_{g}$ as [0, $10^{-20}eV/c^{2}$] and for $\Delta\gamma$ as [$-10^{-7}$, $10^{-7}$].

If deviations are present, indicating non-zero values for $m_{g}$ or $\Delta\gamma$. Specifically, in the context of massive gravity, we hypothesize that an increase in the luminosity distance of the GW event correlates with a more pronounced impact on the GW signal. For WEP violation, we hypothesize that the impact on the GW signal intensifies as the propagation path of the GWs approaches closer to the center of the Milky Way, and also increases with increasing luminosity distance of the GW event.

However, the further the luminosity distance of a GW event, the lower its Signal-to-Noise Ratio (SNR) tends to be. In the presence of the deviations, discrepancies from the standard GW waveforms become more obvious with increasing luminosity distance for both models. But the data precision of these GW events decreases correspondingly. Therefore, we can selectively analyze events with similar SNRs but significant differences in luminosity distance to examine their estimated values and error ranges more accurately. If the deviation exist, it is anticipated that events with greater luminosity distances and similar SNRs will result in more precise estimations of $m_{g}$ and $\Delta\gamma$, with smaller margins of error.

We also can ignore the influence of SNR and use the Bayes factor to directly compare the advantages and disadvantages of the modified models and the standard model using real GW data.

where $\mathcal{Z}_{A}$,$\mathcal{Z}_{B}$ is the Bayesian evidence of models A and B. To enhance the data’s intuitiveness, we use the natural logarithm of the Bayes factor.

Massive gravity and the violation of WEP induce frequency-dependent variations in the propagation speed of GWs. Assuming that the impact of GW waveforms caused by WEP violation is solely attributable to the Milky Way, it follows that as the angle between the GW event and the Milky Way center as observed from Earth gradually decreases, its impact should become more pronounced. Conversely, the impact of GW waveforms caused by massive gravity should be independent of the angles of GW events and would uniformly affect the propagation of GWs through spacetime across all directions. This impact is expected to exhibit cumulative effects over the propagation distance. In other words, as the propagation time of GW increases, so does the luminosity distance of the GW events, the deviation of GW waveforms will become increasingly evident. We can elucidate this phenomenon using Eq. (2), where a larger luminosity distance corresponds to a greater variation in $\delta\Psi$. This variation enables the differentiation of the GW waveform from the one predicted by GR. Thus, this induces an increasing impact on the GW waveform as the luminosity distance of the GW event increases.

Thus, as the luminosity distance of GW events increases, the Bayes factor between the massive gravity and the standard model derived from GWs data will increase. We also hypothesize that the potential impact of the WEP violation would be magnified as the luminosity distance of GW events increases and the angle decreases. Consequently, the Bayes factor comparing WEP violation model to the standard model based on GW data should exhibit higher values, indicating more pronounced WEP violation effects. Conversely, lower values indicate weaker effects.

Our data is based on GW events of BBH in GWTC-2.1 and GWTC-3. We select GW events that exhibit a SNR greater than 10. We calculate the angle between each GW event and the center of the Milky Way as observed from Earth. Using the median of the posterior distribution of angles as the reference value, we approximately uniformly select GW events across the entire range of angles. In this selection, events within 0 to 30 degrees and 150 to 180 degrees are relatively rare. For events with similar angles, we select the one with higher SNR to do the analysis. Following this methodology, we then obtain 22 suitable GW events.

Fig. 1 and Fig. 2 illustrate the posterior distributions of $m_{g}$ and $\Delta\gamma$ for a subset of events (90% credible interval), respectively. The prior distribution of $m_{g}$ follows the distribution presented in Eq. (8). The prior distribution of $\Delta\gamma$ follows the distribution presented in Eqs.(6) in [51], which provides continuous coverage for both negative and non-negative $\Delta\gamma$ value. The detailed posterior information of $\Delta\gamma$ and $m_{g}$ for each event is tabulated in Table  1. Furthermore, the potential data-quality issues inherent in GW observations have been discussed in studies[7, 8]. It is crucial to acknowledge that our method could not measure the exact values of $m_{g}$ and $|\Delta\gamma|$. This limitation arises from the assumption that all observation errors are attributable to the deviation arise from massive gravity and WEP violation. Consequently, the results presented here represent an upper limit for $m_{g}$ and $|\Delta\gamma|$, implying that the actual values of $m_{g}$ and $|\Delta\gamma|$ should be lower than the estimations obtained in this study.

In Fig. 3, contrary to expectations, GW events at a greater luminosity distance and higher SNR cannot provide more tighter constraints on $m_{g}$ than other events. Specifically, high-SNR event such as GW190521_074359, the observed outcome cannot demonstrate superior results in terms of estimated $m_{g}$ compared to other low-SNR events. In cases with similar SNR, events with greater luminosity distances, such as GW190519_153544 and GW190512_180714, their outcomes still fail to demonstrate superior results in terms of estimated $m_{g}$ compared to other events.

In Fig. 4, we do not observe the anticipated trend where GW events with smaller angles and higher SNRs would exhibit tighter constraints on $\Delta\gamma$ than other events. Specifically, high-SNR event GW190521_074359, its outcome cannot demonstrate superior results in terms of estimated $\Delta\gamma$ compared to other events. In cases with similar SNR, GW190512_180714 with small angle, the outcome still fails to demonstrate superior results in terms of the estimated $\Delta\gamma$ compared to other events.

For the analysis of massive gravity, we select four groups of GW events with similar SNRs differing by only 0.1, but exhibiting substantial variations in luminosity distance as shown in Fig. 5. The GW events are arranged from left to right in increasing order of SNR, with events on the left having a smaller luminosity distance than those on the right in each group. Due to the logarithmic scale used for the coordinates of $m_{g}$, it is essential to focus on the upper limit of the 90% credible interval of the error bars, rather than their total length. We observed that at the SNR of 12.5, the upper limit for the $m_{g}$ at the 90% credible interval of GW191222_033537 at a luminosity distance of $3000^{+1700}_{-1700}$ Mpc is tighter than that of GW190513_205428 at $2210^{+550}_{-500}$ Mpc. The differences in the upper limit for the $m_{g}$ at the 90% credible interval for other groups do not show significant variability.

For the analysis of WEP violation, we selected four groups of GW events with similar SNRs differing by only 0.1, but exhibiting substantial variation in the angles between the GW event and the center of the Milky Way, as shown in Fig. 6. The GW events are arranged from left to right in increasing order of SNR, with events on the left having a smaller angle than those on the right in each group. Here, it is important to consider both the upper and lower limits of the 90% credible interval of the error bars. At the SNR = 13.2, GW190915_235702 with an angle of ${94.6^{+15.1}_{-11.1}}$ degrees establishes a tighter limit at the 90% credible interval on the $|\Delta\gamma|$ than the GW event GW190602_175927 with an angle of ${122.1^{+19.8}_{-13.8}}$ degrees. The variations in the upper limit for the $|\Delta\gamma|$ at the 90% credible interval for the remaining groups do not show significant variability.

In fact, potential degeneracies between the deviations from massive gravity or WEP violation and other source parameters could introduce imprecision in the estimated values of $m_{g}$ and $\Delta\gamma$. However, we present upper limits, the results are considered acceptable. Fig. 7 presents the relationship between the chirp mass and $|\Delta\gamma|$ within our posterior distributions.

Additionally, we employ Bayes factors to evaluate the validity of the massive gravity and WEP violation models relative to the standard GW model. Fig. 8 presents the natural logarithm of Bayes factors between the WEP violation model and the standard model for selected GW events of BBH in GWTC-2.1 and GWTC-3. The natural logarithm of Bayes factors are predominantly below an absolute value of 2 and do not show any significant trends in the top-left or bottom-right corners of the distribution, which correspond to the scenarios of small angle with large luminosity distance and large angle with small luminosity distance. These outcomes indicate that significant deviations in Bayes factors are not evident, suggesting a negligible impact of WEP violation on GW propagation under the current experimental conditions and data precision in our study.

Fig. 9 presents the natural logarithm of Bayes factors between the massive gravity and the standard model for selected GW events of BBH in GWTC-2.1 and GWTC-3. Similar to the previous findings, the natural logarithm of Bayes factors are generally below an absolute value of 2, with no obvious deviation that correlate with variations in luminosity distance. These outcomes suggest the absence of significant deviations in Bayes factors relative to the luminosity distance under the current experimental conditions and data precision, which means that the massive gravity cannot provide a superior fit to the data compared to the standard model.

Further, we also calculate Bayes factors between the WEP violation model and massive gravity, as shown in Fig. 10. The Bayes Factors are below zero in approximately two-thirds of the events. However, the absolute values of these natural logarithm of Bayes factors are generally less than 2. These results indicate that there is insufficient evidence to prefer either model. Therefore, the GW data do not provide sufficient grounds to assert the superiority of one model over the other in terms of its fit to the GW observations.

Since the direct detection of GWs [59], the LIGO and Virgo collaboration has rigorously evaluated the consistency between observed GW signals [52, 29, 30] and the theoretical predictions derived from GR. Therefore, previous studies have already placed bounds on the graviton mass and tested potential deviations induced by violation of the WEP during the propagation of GWs, and no significant deviations from the predictions of GR have been revealed through the examination of the GW data. It should be noted, however, that under the current experimental conditions, it is still necessary to use some novel methods to examine the existence of graviton mass and to test WEP violation during the propagation of GWs.

In this study, we do not find a significant relationship between the parameter $m_{g}$ and luminosity distance and a relationship between the parameter $|\Delta\gamma|$ and both luminosity distance and sky location of the GW events. We do not find a significant preference for either the massive gravity or WEP violation models. Moreover, Bayes factors fail to provide obvious evidence favoring one model over the other. With the current level of experimental precision, it remains uncertain to determine whether these deviations exist.

In this study, we have exclusively focused on the scenarios where $\alpha=0$ and $A_{0}>0$ in modified dispersion relations. In an upcoming work, we will consider more cases of $A_{\alpha}$ with different $\alpha$, and will adopt a more comprehensive theory that incorporates the effects of WEP and LI violations, which will allow us to directly compare the impact of these two effects on GWs. In the method of detecting the properties of GWs, we will no longer just focus on phase variations in the frequency domain of GW waveform. However, we will investigate how violations of WEP and LI may influence on other properties of GWs, such as their polarization.

The fourth observing run (O4) of the LIGO-Virgo-KAGRA (LVK) GW detector network has started running, it is promising to get more GW events, including more binary neutron star (BNS) merger events. These events are anticipated to have a promising probability of presenting multi-messenger characteristics, which would permit a better determination of parameters such as right ascension, declination, and luminosity distance. Then, we can clearly analyze the differences in propagation speeds between GWs and electromagnetic waves and the impact of large mass gravitational sources on the propagation of GWs. So it is possible to better constrain the upper bound of $m_{g}$ and $\Delta\gamma$.

In the future, the implementation of space-borne interferometers will expand our ability to detect GWs from a wider array of sources and improve the SNRs of detected events [60, 61, 62, 63, 64]. We will also develop statistical methods and try more efficient ways to search for traces of deviations from GR in the multitude of GW events. The utilization of O4 data and forthcoming space-borne interferometers will allow us to detect various GW events, and conduct population analysis on GW events that may contain deviations from GR. Hence, it is worthwhile to discuss the capability of GW detectors to detect signals that exhibit such deviations in future studies.