Bounding the speed of gravity with gravitational wave observations

The time delay between gravitational wave signals arriving at widely separated detectors can be used to place upper and lower bounds on the speed of gravitational wave propagation. Using a Bayesian approach that combines the first three gravitational wave detections reported by the LIGO collaboration we constrain the gravitational waves propagation speed c_gw to the 90% credible interval 0.55 c<c_gw<1.42 c, where c is the speed of light in vacuum. These bounds will improve as more detections are made and as more detectors join the worldwide network. Of order twenty detections by the two LIGO detectors will constrain the speed of gravity to within 20% of the speed of light, while just five detections by the LIGO-Virgo-Kagra network will constrain the speed of gravity to within 1% of the speed of light.

The first detections of gravitational waves from merging black hole binaries [1][2][3] have been used to test many fundamental properties of gravity [3][4][5][6], and have been used to place the first observational upper limit on the speed of gravitational wave propagation [7]. In this letter we set a more stringent upper limit on the gravitational waves propagation speed c gw by combining all the detections announced to-date, and by applying a full Bayesian analysis. We also provide the first direct lower bound on the propagation speed: c gw > 0.55 c at 95% confidence. While there are strong theoretical arguments that demand c gw ≥ c to prevent gravitational Cherenkov radiation [8], the LIGO detections provide the first direct observational constraints.
Gravitational waves generically propagate at a speed different from c and with frequency dependence dispersion relations in theories of modified gravity, see e.g. Refs. [6,7,[9][10][11][12]. Thus, a precise determination of c gw is a test of gravitation complementary to other observations. To quantify what 'precise' tests mean for General Relativity, let us recall that some post-Newtonian parameters are known to O(10 −4 ) [13] while cosmological or other astrophysical observations typically constrain modifications to General Relativity at the O(10 −2 ) level [14,15].
A convenient parametrization for theories preserving rotation invariance is to write the dispersion relation as where m g refers to the mass of the graviton, c gw is what we call 'speed' of gravitational waves, and the rest of operators are wavenumber-dependent modifications suppressed by a high-energy scale Λ (for a parametrization in scenarios breaking rotation invariance, see e.g. Ref. [10]). Both m g and Λ can be constrained by the absence of dispersion of the waves traveling cosmological distances. The scale Λ is already constrained to be very large [9], making it very difficult to constrain the operator a. For the graviton mass the LIGO Scientific & Virgo Collaborations put the strong bound m g < 7.7 × 10 −23 eV/c 2 [3]. However, the parameter c gw can not be tested by dispersion measurements and other methods are required [7]. Measuring c gw : In the following we focus on possible ways to directly measure c gw . Since the signals measured by LIGO are dominated by the signal-to-noise accumulated in a narrow band between 50 Hz -200 Hz, our time delay bounds can be interpreted as constraints on the speed of gravity at a frequency f ∼ 100 Hz. Since the LIGO bounds constrain dispersion effects to be small over hundreds of Mpc, they can safely be ignored on the terrestrial distance scales we are considering. Note that the inference that the observed signals come from hundreds of Mpc away relies on waveform models derived from General Relativity, and may not apply to a theory that predicts c gw = c.
The most obvious way to measure the speed of gravitational wave propagation is to observe the same astrophysical source using both gravity and light. However, for the three gravitational wave detections that have been announced thus far no unambiguous electromagnetic counterparts have been detected, and a different approach must be taken to constrain c gw . The finite distance between the Hanford and Livingston gravitational wave detectors can be used to set an absolute upper limit on the propagation speed [7] since the observed gravitational wave signals did not arrive simultaneously in the two detectors. Here we show that a proper statistical arXiv:1707.06101v2 [gr-qc] 5 Oct 2017 treatment that folds in the probability distribution of the time delays as function of c gw also allows us to set lower bounds on the propagation velocity. It should be noted that when the first confirmed electromagnetic counterpart to a gravitational wave signal is finally observed, the bounds on the difference in propagation velocities, |c gw −c| will be many orders of magnitude more stringent than what we can ever hope to set using gravitational wave signals alone [9,[16][17][18]. Precisely for the same reason, this identification may never happen if the speed difference is not very small, since that would mean a significant time offset. For other possible stringent modelindependent bounds not relying on the detection of a counterpart see [11,19].
Constraints on c gw from LIGO detections: The LIGO gravitational wave detectors at the Hanford and Livingston sites are separated by a light-travel time of t 0 = 10.012 ms. The time delay for light along a propagation direction that makes an angle θ with the line connecting the two sites is ∆t EM = t 0 cos θ. For sources distributed isotropically on the sky, there are equal numbers of sources per solid angle element d cos θ dφ, thus the time delays for electromagnetic signals are uniformly distributed with p( The gravitational wave time delay is given by ∆t = (c/c gw )∆t EM , thus the probability of observing a time delay ∆t between gravitational wave signals arriving at the two sites for sources uniformly distributed on the sky is given by the likelihood While the sources may be uniformly distributed on the sky, the antenna patterns of the detectors make it more likely to detect systems above or below the plane of detectors. Assuming roughly equal sensitivity for the detectors, the observational bias scales as F 3 , where F (θ, φ) 2 = D=H,L F D +,× (θ, φ) 2 is the polarization averaged network antenna pattern [20]. The resulting distribution of electromagnetic time delays is then well fit by .65 ms. We use this modified distribution to define the likelihood p(∆t|c gw ). For multiple events the full likelihood is the product of the per-event likelihoods. The posterior distribution for c gw follows from Bayes' theorem: p(c gw |∆t) = p(∆t|c gw )p(c gw )/p(∆t). We consider two possibilities for the prior on the speed of gravity, p(c gw ): flat in c gw and flat in ln c gw in the interval c gw ∈ [c L , c U ]. For the results shown here we set c U = 100 c, and either c L = c/100, or c L = c. The latter limit takes into account the Cherenkov radiation constraint [8]. For three or more events the choice of prior has very little impact on the upper limit. To account for the measurement error in ∆t we use a Markov Chain Monte Carlo to marginalize over the errors in the arrival times. The first detections of black hole mergers by LIGO provide measurements of ∆t that were quoted in terms of central values and 90% credible intervals. Since the full posterior distributions for ∆t were not provided, we assume that the distributions can be approximated as normal distributions with mean µ and standard deviation σ with values: GW150914 (µ = 6.9 ms, σ = 0.30 ms) [5]; GW151226 (µ = 1.1 ms, σ = 0.18 ms) [5]; GW170104 (µ = 3.0 ms, σ = 0.30 ms) [3] (for a discussion about this assumption, see e.g. Ref. [21]). The upper bound on c gw quoted in Ref. [7] was found by taking the minimum time delay from GW150914 as ∆t = µ − 2σ = 6.3 ms, and demanding that c gw < c t 0 /∆t = 1.6 c. Note that this value is lower than the bound of 1.7 c quoted in Ref. [7] as they interpreted the error in ∆t quoted in Ref. [1] as one-sigma errors, when in fact they were the bounds on the 90% credible interval. We compute the posterior distribution for the gravitational wave propagation speed, p(c gw |∆t) using a Markov Chain Monte Carlo algorithm that marginalizes over the uncertainties in the time delays by drawing new values of ∆t from the assumed posterior distributions at each iteration of the Markov chain. Figure 1 shows the posterior distributions for c gw using each of the detections separately. Individually the three events yield 95% upper bounds on the propagation velocity for the linear and (log) uniform priors of c gw < 1.37 c (1.26 c) for GW150914; c gw < 10.1 c (8.57 c) for GW151226; and c gw < 3.19 c (2.94 c) for GW170104. Each event also yields a 95% lower bound on the propagation velocity, but these limits are not very interesting since all of the distributions have some support at c gw 0. Note that GW151226 produces the weakest upper bound on the propagation velocity even though it has the most accurately measured time delay. This is because the strongest upper limits come from events with the longest time delay, and even allowing for the uncertainties in the time delay measurements for GW150914 and GW170104, both are constrained to have delays that are much longer than for GW151226. Figure 2 shows the posterior distribution for c gw found by combining all three LIGO detections together for uniform priors in c gw or ln c gw . For the wider prior range with c L = c/100 the combination of the three detections yield an interesting lower bound on the propagation speed. The 90% credible interval for the linear and log priors are 0.55 c < c gw < 1.42 c and 0.41c < c gw < 1.39 c respectively. The upper limit is only weakly dependent on the choice of prior distribution. Notice that the form of the posterior distributions can be qualitatively understood analytically by using Bayes' theorem and ignoring the complications from antenna patterns and noise. For a uniform prior on c gw this gives p(c gw |∆t 1 , ∆t 2 , ∆t 3 ) = 4c 3 gw /c 4 * for c gw < c * , where c * = ct 0 /∆t 1 and ∆t 1 is the longest of the observed time delays. For the uniform-in-log prior the posterior distribution is p(c gw |∆t 1 , ∆t 2 , ∆t 3 ) = 3c 2 gw /c 3 * for c gw < c * . This explains the growing of the posterior in Figure 2, while the smoothing of the curve is related to the Gaussian noise. While an analytical understanding of the posterior distribution is possible, it is necessary to use numerical sampling methods to compute the full shape that includes marginalization over observational errors, possible orientation biases, and add new detectors and detections.
One possible limitation of using the published LIGO results, rather than analyzing the raw data, is that the standard LIGO searches exclude signals with Hanford-Livingston time delays greater 15 ms [22], so potentially missing some signals if c gw < 0.66c. On the other hand, loud single detector triggers consistent with a binary merger would not go un-noticed, as evidence by GW170104, which was missed by the standard search due to the Hanford detector being incorrectly flagged as out of observing mode, and only found later in an analysis of single detector triggers from the Livingston detector [3]. It is highly unlikely that pairs of triggers with time delays greater than 15 ms would be overlooked, especially if they shared similar parameters and occurred within minutes of each other.
Forecasts for more detections and more detectors: It is interesting to consider how the bounds will improve with additional detections and detectors, using just the gravitational time delays (as mentioned earlier, combined electromagnetic and gravitational observations will dramatically improve the measurements). The upper bound is mostly set by detections with large time delays, while the lower bound is set by having many signals with a wide range of delays. For the two-detector LIGO network, and assuming c gw = c as predicted by General Relativity, we should see one event with a time delay ∆t > 8.6 ms with ten detections, and one event with a time delay ∆t > 9.6 ms with 100 detections. Just those single events would yield 99% upper limits better than c gw < 1.2 c and c gw < 1.07 c respectively, assuming that ∆t is measured to a level of accuracy typical of the first three detections (the accuracy with which ∆t can be measured depends on the signal-to-noise ratio and the number of cycles completed in-band, among other things). Performing multiple Monte-Carlo simulations under the assumption that General Relativity is the correct theory of gravity indicates that with 100 detections by the 2detector LIGO network we will be able to constrain c gw to within a few percent of the speed of light for both the upper and lower bounds.
Far better constraints can be achieved with far fewer events by using a larger network of detectors. In the next few years the LIGO Hanford (H) and Livingston (L) detectors will be joined by the Virgo (V) detector in Italy, the Kagra (K) detector in Japan [23], and later, another LIGO detector in India [24]. With an N detector network there are N (N − 1)/2 time delays. From the geopositions of the detector sites [25] the maximum electromagnetic time delays between these sites are: HL = 10.012 ms, HV = 27.288 ms, HK = 25.158 ms, LV = 26.448 ms, LK = 32.455 ms, V K = 29.202 ms. (LIGO India would further improve the bound, but will not be available at least 2025). Upper bounds on the speed of gravity are dominated by events with sky locations that come close to maximizing the electromagnetic time delay between a pair of detectors. For the HL network just 5% of events are within 95% of the maximum time delay, while for the HLVK network 25% of events are. Thus, on average, it only takes a few events to produce tight limits using the larger network of detectors. Complete information about the inter-site time delays is contained in the joint probability distribution of the N − 1 electromagnetic time delays between one reference detector and the other detectors in the network. Figure 3 shows slices through the joint electromagnetic time delay distribution for the HLVK network using Hanford as the reference site assuming a uniform distribution of sources. Here we did not correct for the observational bias, since the network antenna pattern for a four detector network is fairly uniform. Applying the change of variable ∆t = (c/c gw )∆t EM to this distribution as we did for the two-detector HL case yields the joint likelihood p(∆t HL , ∆t HV , ∆t HK |c gw ). Using simulated detections of events measured to a precision of σ = 0.3 ms in each detector, we find that with just 3 detections the HLVK network will typically be able to constrain c gw to the 99% credible region c gw /c = 1.00 ± 0.02. The constraints improve to better than 1% of the speed of light with 5 detections.
FIG. 3. Slices through the joint electromagnetic time delay distribution for the HLVK network using Hanford as the reference site. Darker colors in the two-dimensional slices indicate higher density. Note the cup-like structure of the distributions (higher at the edges than in the center).
Summary: Combining the time delay measurements between detector sites for multiple gravitational wave events can be used to place interesting constraints on the speed of gravity. The LIGO detections made to-date already constrain the speed of gravity to within 50% of the speed of light. Additional LIGO detections in the next few years should improve the bound to of order 10%. The bounds will improve rapidly as more detectors join the worldwide network, with just a half-dozen detections by the Hanford-Livingston-Virgo-Kagra network constraining deviations to better than 1%. These bounds will allow to test General Relativity to the level of other standard tests, as those coming from the damping of orbits in binary systems or cosmology.