## Measuring the speed of gravity and the cosmic expansion with time delays between gravity and light from binary neutron stars

Leonardo Iampieri , 1, 2, ∗ Simone Mastrogiovanni , 2 and Francesco Pannarale 1, 2

1 Dipartimento di Fisica, Universit` di Roma 'Sapienza,' Piazzale A. Moro 5, I-00185, Roma, Italy a 2 INFN, Sezione di Roma, I-00185 Roma, Italy

(Dated: March 7, 2025)

The first observation of a gravitational wave (GW) and a short gamma-ray burst (sGRB) emitted by the same binary neutron star (BNS) merger o ffi cially opened the field of GW multimessenger astronomy. In this paper, we define and address lagging sirens , a new class of multimessenger BNSs for which associated GWs and sGRBs are observed without the identification of their host galaxy. We propose a new methodology to use the observed time delay of these sources to constrain the speed of gravity that is, the propagation speed of gravitational waves, the Hubble constant and the prompt time delay distribution between GWs and sGRBs, even though a direct redshift estimation from the host galaxy is unavailable. Our method exploits the intrinsic relation between GWs and sGRBs observed and prompt time delays to obtain a statistical redshift measure for the cosmological sources. We show that this technique can be used to infer the Hubble constant at the 10% level of precision with future-generation GW detectors such as the Einstein Telescope and only 100 observations of this kind. The novel procedure that we propose has systematics that di ff er completely from the ones of previous GW methods for cosmology. Additionally, we demonstrate for the first time that the speed of gravity and the distribution of the prompt time delays between GWs and sGRBs can be inferred conjointly with less than 10 sources even with current GW detector sensitivities.

## I. INTRODUCTION

Gravitational waves (GWs) from compact binary coalescences (CBCs) can be used to probe the expansion of the Universe. Given the current open tensions in cosmology [1], it is crucial to measure the Universe expansion parameters, particularly the Hubble constant H 0 , with independent classes of sources. CBCs detected via GWs are the only cosmological sources for which a direct measurement of the luminosity distance without the need for an astrophysical calibration is possible [2, 3]. For this reason, GW sources are referred to as 'standard sirens.' Unfortunately, GWs do not directly provide the source redshift which is the second ingredient required to measure the cosmic expansion parameters. Current scientific literature proposes and uses several methodologies to assign redshifts to GW events and infer the cosmic expansion parameters.

for dark sirens cosmology is the source-frame mass method [12-24]. This method consists of exploiting the di ff erence between the detector mass md and source mass ms of CBCs, which are related by a redshift factor md = (1 + z ms ) . Thus, by knowing the source-frame mass spectrum, it is possible to infer statistically the source redshift. Since the source mass spectrum is not known accurately, current studies conjointly fit flexible phenomenological models of the mass spectrum along with cosmological parameters that can in some limits approximate the true binary distributions [25]. The source frame mass method was used with the latest GW detections to obtain a value of H 0 = 68 + 12 -7 kms -1 Mpc -1 [26]. The low precision in the inferred value of H 0 is due to the broadness of the source-frame mass spectrum, which prevents us from obtaining precise redshift estimations for the GW event sources used to measure H 0 .

In terms of analysis, the most trivial scenario is the one of bright sirens , where an electromagnetic (EM) counterpart is observed conjointly with the GW emission. This is the case of the binary neutron star (BNS) merger GW170817 [4, 5]. For these cases, it could be possible to identify the host galaxy and obtain a spectroscopic estimate of the galaxy redshift. The observation of GW170817 and its EM counterpart was exploited to obtain H 0 = 70 + 19 -8 kms -1 Mpc -1 [6, 7]. The error budget on the H 0 measure is entirely dominated by the luminosity distance error from the GW detection [8], as the redshift is accurately measured.

So far, GW170817 is the only GW event with an unambiguously associated EM counterpart. Out of almost 100 GW detections achieved to date [9-11], most are due to binary black hole mergers observed without an EM counterpart ( dark sirens ). Among the several techniques proposed

∗ leonardo.iampieri@uniroma1.it

Following the same idea of the source-frame mass method, this paper proposes a new method for GW cosmology: the prompt time delay method. The methodology consists of obtaining an implicit source redshift information for GW sources with observed short gamma-ray bursts (sGRBs) but no identified host galaxy. We refer to this type of source as lagging sirens . Currently, there are no GW detections of this kind. However, we certainly expect them for future GW detectors such as Einstein Telescope (ET) [27, 28]. ET will detect BNSs up to redshift 2-3 [29]. Although the farthest GRBs, and their afterglows, could be detected even at these redshifts [30], it might not be trivial to identify and obtain a spectroscopic measurement of its host galaxy.

The paper is organized as follows. In Sec. II, we introduce the relations between prompt and observed GW-sGRB time delays considering also possible deviations between the speed of light and gravity. We further introduce the statistical framework used to infer the Hubble constant H 0 , speed of gravity ˆ α 0 and prompt time delay distribution. In Sec. III we simulate catalogs of GW and sGRB observations consider-

ing two scenarios for the sensitivity of GW detectors, namely, the sensitivity expected to be achieved in 2027-2030 with the fifth LIGO-Virgo-KAGRA Observing Run (O5) [31-34] and an ET-like sensitivity expected to be achieved in 2035 + . In Sec. IV we present the prospects of measuring H 0 , ˆ α 0 and the prompt time delay distribution. Finally, in Sec. V we draw our conclusions.

## II. COSMOLOGICAL AND STATISTICAL MODELS

This section provides the background knowledge required to build the prompt time delay method. Specifically, in Sec. II A, we introduce the relations between cosmological expansion, redshift and prompt time delay distributions, while in Sec. II B we review the Bayesian statistical framework used in this work.

## A. Theoretical background

We model the cosmic expansion with a flat Λ CDM model described by a Hubble constant H 0 , dark matter energy density Ω m and dark energy density ΩΛ = 1 -Ω m . In this model, the luminosity distance D L and redshift z of cosmological sources are related by

$$D _ { \mathrm L } ( z ) = \frac { c ( 1 + z ) } { H _ { 0 } } \int _ { 0 } ^ { z } \frac { d z ^ { \prime } } { \sqrt { \Omega _ { m } ( 1 + z ^ { \prime } ) ^ { 3 } + ( 1 - \Omega _ { m } ) } } \,, \quad ( 1 ) \quad \begin{pmatrix} \tilde { \alpha } _ { 0 } & \ell \\ \text{ of th} \end{pmatrix}$$

where c is the speed of light.

We define the prompt time delay ∆ t s as the elapsed time between the GW luminosity peak at the CBC merger and the emission of the sGRB. A positive ∆ t s indicates that the sGRB is emitted after the CBC merger, while a negative value indicates that it is emitted prior to the merger. The observed time delay between GW and sGRB at the merger is [35]

$$\Delta t _ { \text{d} } = ( 1 + z ) \Delta t _ { \text{s} } + \frac { \hat { \alpha } _ { 0 } } { 2 } \text{T} _ { l } \,, \quad \quad ( 2 ) \quad \text{wh} \\ \quad \text{the}$$

where T l is the lookback time

$$T _ { l } = \frac { 1 } { H _ { 0 } } \int _ { 0 } ^ { z } d z ^ { \prime } \frac { 1 } { \sqrt { \Omega _ { m } ( 1 + z ^ { \prime } ) ^ { 3 } + ( 1 - \Omega _ { m } ) } } \,, \quad ( 3 ) \quad \text{ is no}$$

and ˆ 0 α is a scalar factor proportional to the di ff erence between the GW's speed v GW and the speed of light relative to the speed of light itself, namely

$$\hat { \alpha } _ { 0 } = 2 \frac { v _ { G W } - c } { c } \,. \quad \quad \quad \ ( 4 ) \quad \ s p _ { i o i } ^ { \ } \,.$$

In writing Eq. (2), we implicitly assumed that the speed of gravity does not depend on the GW frequency and that its deviation from the speed of light is small [35]. Moreover, we assumed that the speed of gravity does not change during cosmic time, although this possibility is not excluded in principle for some theories of gravity [36, 37] or agnostic tests on the speed of gravity [38].

FIG. 1. Observed time delays (vertical axis) as a function of redshift (horizontal axis) for di ff erent values of the speed of gravity parameter and prompt time delay. The di ff erent lines indicate various choices for the prompt time delay and speed of gravity.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00362v2-figure-1318aa90-1e0f-4aa3-ba3c-926edb15ce7a.jpg)

Equation (2) can be used to predict the observed time delays between the GW and sGRB signals from the same event. For instance, if the GW and sGRB were prompted together ( ∆ t s = 0 s), and GWs were slightly slower than the light ( ˆ α 0 &lt; 0), then the GW would be detected after the arrival of the sGRB ( ∆ t d &lt; 0 s). If instead the sGRB were prompted after the GW ( ∆ t s &gt; 0 s) and the speed of gravity were equal to c ( ˆ α 0 = 0), then the GW would surely be detected prior to the sGRB ( ∆ t d &gt; 0 s). Figure 1 shows the observed time delays as a function of redshift for various values of the speed of gravity and prompt time delay. The larger the redshift, the longer the observed time delay between the GW and the sGRB is. From Fig. 1, we can already see that if the observed time delay were accurately measured 1 and the prompt time delay and speed of gravity were known, one could get implicit redshift information for the GW-sGRB sources. Unfortunately, while we can argue that the speed of gravity must be equal to the speed of light if we assume classical general relativity, the distribution of prompt time delays between GWs and sGRBs is not a constant and it is currently unknown. For this reason, we introduce a statistical framework able to infer it.

## B. Statistical background

Let us assume that we have a set of common population parameters ⃗ λ , such as the Hubble constant or the propagation speed of GWs, that we want to infer from N obs GW-sGRB joint observations collected from the data { ⃗ , . . . , d 1 ⃗ dN obs } . We further assume that each astrophysical source is described by

1 This is a reasonable assumption also supported by the measure of GW170817 with its sGRBs for which an observed time delay of ∆ t d = 1 74 . ± 0 05 s was observed [4]. .

a set of parameters ⃗ θ , the distributions of which are described by a CBC rate model

$$\frac { d N } { d t _ { d } d \vec { \theta } } ( \vec { \theta } | \vec { \lambda } ) \,,$$

where t d is the detector time. Our aim is to infer the parameters ⃗ λ by conjointly fitting the rate distribution of the parameters ⃗ θ . The hierarchical likelihood of obtaining these detections is [39]

$$\mathcal { L } ( \{ \vec { d } _ { i } \} | \vec { \lambda } ) = \prod _ { i = 1 } ^ { N _ { o b s } } \frac { \int \mathcal { L } ( \vec { d$$

where L ( ⃗ d i | ⃗ ⃗ θ λ , ) is the measurement likelihood, which quantifies how well we can measure ⃗ θ , conditioning on ⃗ λ , from a single event detection ⃗ d i , while p det ( ⃗ ⃗ θ λ , ) is a detection probability used to model the finite sensitivity of GW and sGRB detectors.

For this study, the relevant ⃗ θ parameters are the luminosity distance D L and the observed GW-sGRB time delay ∆ t d . Thus, the CBC rate model written in terms of these two quantities is

$$\frac { d N } { d t _ { d } d \vec { \theta } } \equiv \frac { d N } { d t _ { d } d D _ { L } d \Delta t _ { d } } \,. \quad \quad \quad$$

Since we want to exploit the prompt time delay distribution ∆ t s to obtain an implicit redshift estimation z , we need to rewrite the CBC rate in terms of these quantities. This can be done by performing a change of coordinates from luminosity distance and observed time delay to redshift and prompt time delay, namely

$$\frac { d N } { d t _ { d } d D _ { L } d \Delta t _ { d } } \equiv \frac { d N } { d t _ { s } d z d \Delta t _ { s } } \frac { 1 } { 1 + z$$

where the (1 + z ) -1 factor arises from the transformation between source frame t s and detector frame time t d , d N d sd d t z ∆ t s is the source frame rate per redshift z and per prompt time delay ∆ t s , while det JD → S is the determinant of the Jacobian associated to the change of variables between the detector and source frames (including also GW-sGRB time delays), that is,

$$\det J _ { D \rightarrow S } = \left | \frac { \partial D _ { L } } { \partial z } \frac { \partial \Delta t _ { d } } { \partial \Delta t _ { s } } \right | = \left | \frac$$

where

$$\frac { \frac { \partial D _ { \L } } { \partial z } = \frac { D _ { \L } ( z ) } { 1 + z } + \frac { c ( 1 + z ) } { H _ { 0 } } \frac { 1 } { \$$

We additionally parameterize the source frame rate per redshift z and per prompt time delay as

$$\frac { d N } { d t _ { s } d z d \Delta t _ { s } } & = \frac { d N } { d t _ { s } d V _ { c } d \Delta t _ { s } } \frac { d V _ { c } } { d z } = \\ & = R _ { 0 } \psi ( z ; \vec { \lambda } ) p _ { \text{pop} } ( \Delta t _ { s } | \vec { \lambda } ) \frac { d V _ { c } } { d z }, \\ \therefore R _ { 0 } \text{ is the BNS merger rate per comoving volume at red-}$$

where R 0 is the BNS merger rate per comoving volume at redshift z = 0, ψ ( z ; ⃗ λ ) is a phenomenological parametrization for the BNS merger rate as a function of redshift such that ψ ( z = 0; ⃗ λ ) = 1, p pop is a probability density function describing the prompt time delay distribution, and

$$\frac { d V _ { c } } { d z } = 4 \pi \left [ \frac { c } { H _ { 0 } } \right ] ^ { 3 } \left [ \int _ { 0 } ^ { z } \frac { d z ^ { \prime } } { E ( z ^ { \prime } ) } \right ] ^ { 2 } \text{ \quad \ \ } ( 1 2 )$$

is the di ff erential of the comoving volume.

The overall detector rate model written in terms of source rate quantities can be rewritten as

$$\underset { \Delta t _ { d } } { \text{nos} } \quad \frac { d N } { d t _ { d } d D _ { L } d \Delta t _ { d } } = R _ { 0 } \psi ( z ; \vec { \lambda } ) p _ { \text{pop} } ( \Delta t _ { s } | \vec { \lambda } ) \frac { d V _ { c } } { d z } \frac { 1 } { ( 1 + z ) ^ { 2 } \left | \frac { \partial D _ { L } } { \partial z } \right | } \,. \ ( 1 3 )$$

By substituting Eq. (13) in Eq. (6), we find the full expression of the hierarchical likelihood

$$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{	\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{ $\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{  $\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text	$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text
$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text#\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{#\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$
$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\test{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\test$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text}$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text'$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text<$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text.$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text`$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$$

To evaluate the hierarchical likelihood, we use icarogw [40], a software able to calculate an approximate value of Eq. (14) starting from a set of posterior samples on luminosity distances and on observed time delays, and detectable injections used to evaluate the selection bias.

## C. Study case with GW170817 and GRB 170817A

To display an application of this method, we recalculate the constraints on the speed of gravity obtained in [4] from the joint observation of GW170817 and GRB 170817A. We consider two cases: the first in which we only consider the luminosity distance estimated from the GW signal and the observed time delay between GW170817 and GRB 170817A and the second in which we also include the redshift information from the event host galaxy. To include the redshift information from the host galaxy, the hierarchical likelihood in Eq. (14) can be modified as follows

$$\mathcal { L } ( | \vec { d } _ { i } | | \vec { x } ) \in \prod _ { i = 1 } ^ { N _ { \infty } } \frac { \int \mathcal { L } ( \vec { d } _ { i } | D _ { L }, \Delta t _ { d }, \vec { x } ) \mathcal { L } ( \vec { d } _ { l } ^ { \infty } | | z, \sigma _ { z } ) \frac { \vec { d } _ { z } } { \vec { d } _ { z } } \frac { \phi ( z ; \vec { x } ) p _ { p p } ( \Delta t _ { l }, \vec { l } ) } { ( 1 + z ) ^ { 2 } | ^ { \frac { \vec { d } p } { \sigma } } | ^ { \frac { \vec { d } p } { \sigma } } | } } { \int p _ { \det } ( D _ { L }, \Delta t _ { d }, \vec { x } ) \frac { \vec { d } _ { z } } { \vec { d } _ { z } } \frac { \phi ( z ; \vec { x } ) p _ { p p } ( \Delta t _ { l }, \vec { l } ) } { ( 1 + z ) ^ { 2 } | ^ { \frac { \vec { d } p } { \sigma } } | ^ { \frac { \vec { d } p } { \sigma } } | } } \, d D _ { L } d \Delta t _ { d } \,$$

FIG. 2. Posterior distributions for the speed of gravity from the joint observation of GW170817 and GRB 170817A. The green line represents the posterior distribution using only the measurement of luminosity Distance ( D L) and Time Delay ( ∆ t d ). The blue line represents the posterior distribution when including the redshift information from the event host galaxy.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00362v2-figure-ec26a113-26dc-48b3-9854-bad33f36555e.jpg)

where L ( ⃗ d EM i | z , σ z ) encapsulates the measure of the redshift from the source host galaxy. Following [6], we assume this is a Gaussian distribution with mean z = 0 0100 and . σ = 0 0005. . Note also that in Eq. (15), the detection probability at the denominator should also contain the detection probability of the EM counterpart. However, we neglect this component as during the O2 run, the detection range of electromagnetic counterparts was significantly higher than the one of the GWsources. To estimate the source luminosity distance, and hence the redshift, we use the low-spin parameter estimation samples released 2 with [7]. The estimated value of the luminosity distance from the posterior samples is D L = 41 + 6 -12 Mpc [7]. We use a flat Λ CDMmodelwith cosmological parameters from [41] ( H 0 = 67 66 km s . -1 Mpc -1 , Ω m = 0 310) to convert . between luminosity distance and redshift. We also assume that the time delay at the detector between GW170817 and GRB170817A is ∆ t d = 1 74 . + 0 05 . -0 05 . s, measured with a Gaussian posterior. More importantly, we fix the distribution of source prompt time delays to be uniform between [0,10] s, meaning that the GRB can be emitted up to 10 seconds after the merger. This choice is made following [4], where to derive the constraints on the speed of gravity it was assumed that the GRB could be emitted at most 10 seconds after the merger.

Figure 2 shows the posterior distribution of the fractional di ff erence between the speed of gravity and the speed of light. The boundaries of the posterior can be constrained as a proxy for the constraints on the speed of gravity. When just using GW170817 and GRB 170817A, we obtain

$$- 3 \times 1 0 ^ { - 1 5 } \lesssim \frac { v _ { \text{GW} } - c } { c } \, \lesssim + 7 \times 1 0 ^ { - 1 6 }, \quad ( 1 6 ) \quad \text{the $l$} \\ ( 1, 3 )$$

which matches the constraints obtained in [4] by combining only GW and Fermi-INTEGRAL data. Notably, in [4], the au-

2 https://dcc.ligo.org/LIGO-P1800061/public

thors excluded inputs from external analyses, such as the identification of NGC 4993 via EM observations at lower wavelengths, and avoided assumptions about cosmology required to convert the luminosity distance into redshift. Instead, they assumed a conservative (in terms of speed of gravity constraints) lower limit of 26 Mpc for the luminosity distance. When we also include the information on the redshift of NGC 4993, we obtain

$$- 1 \times 1 0 ^ { - 1 5 } \lesssim \frac { v _ { \text{GW} } - c } { c } \, \lesssim + 3 \times 1 0 ^ { - 1 6 }. \quad \, ( 1 7 )$$

We can see that the posterior generated including the redshift estimation of NGC 4993 has shorter tails than the one generated without redshift. This is due to the fact that, the redshift of GW170817 is more precise than the luminosity distance inferred from the GW. As a consequence, when we fix the cosmology, the luminosity distance of the source is almost entirely given by the galaxy host redshift, with a negligible uncertainty on its value.

## III. MOCKCATALOGS OF GW AND SGRB DETECTIONS

We simulated mock catalogs of GW and sGRB detections to forecast the potential of the method we propose to constrain prompt time delay distributions, speed of gravity and H 0 . In Sec. III A we describe how we simulate joint detections of GW and sGRB, while in Sec. III B we provide a general overview of the observational properties in terms of luminosity distance and observed time delay of the GW-sGRB detections.

## A. Catalogs generation

In Fig. 3 we depict the simulation procedure used to generate the GW-sGRB mock catalogs.

a. Binary neutron star simulation: We choose a flat Λ CDM cosmology with parameters matching those derived from the latest observations of the cosmic microwave background by Planck [41] ( H 0 = 67 66 km s . -1 Mpc -1 , Ω m = 0 310) and we also assume that GWs propagate at the speed . of light ( ˆ α 0 = 0). We simulate BNS sources with an isotropic distribution over the sky in terms of right ascension α and declination δ . We also simulate an isotropic distribution for the binary orbital inclination angle ι and a uniform distribution of polarization angles ψ . The source masses m 1 and m 2 of the binary neutron stars are drawn uniformly in the interval (1 , 3) M ⊙ , with the condition m 2 ≤ m 1 . This distribution is consistent with the latest rate models inferred from the GW detections during the third LIGO-Virgo-KAGRA Observing Run (O3) [42]. The spin magnitudes of the two neutron stars are assumed to be vanishing. This is consistent with expectations for BNS, which are low-spin objects [7]. The distribution of the BNS sources in redshift is modeled according to a

FIG. 3. Flowchart depicting the simulation procedure to generate the GWs and sGRBs mock catalogs.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00362v2-figure-02c64d83-71e8-4092-867a-d06625086b8e.jpg)

Madau-Dickinson-like star formation rate [43], namely

$$\frac { d N } { d V _ { c } d t _ { s } } \in [ 1 + ( 1 + z _ { p } ) ^ { - \gamma - k } ] \frac { ( 1 + z ) ^ { \gamma } } { 1 + \left ( \frac { 1 + z } { 1 + z _ { p } } \right ) ^ { \gamma + k } } \quad ( 1 8 ) \quad \text{the i} \\ \text{sens} \text{i}$$

where the model parameters ( γ , k , zp ) are set to (2.7, 6.0, 2.0). The parameters chosen for the BNS rate function are in agreement with the binary black hole merger rate inferred from O3 [11, 44] as it is currently not possible to constrain the BNS rate as a function of redshift. GW events are simulated up to a maximum redshift that depends on the GW detector network configuration. We chose to simulate BNSs up to a maximum redshift z max = 0 2 for O5, and . z max = 2 for ET. These values are high enough to embrace all sources detectable by the GW detectors with a joint sGRB observation.

- b. Gavitational wave detection: We consider a signal to be detected if its observed detector network signal-to-noise ratio (SNR) is greater than 12. This is a conservative threshold for the detection of GW events [34]. The network SNR is defined as the quadrature sum of the individual detector SNRs

$$\rho _ { \text{net} } ^ { 2 } = \sum _ { i \in \text{IFO} } \rho _ { \text{opt,} i } ^ { 2 } \,. \quad \quad \ \ ( 1 9 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

In this work, we model the GW SNR ρ opt , i at the detector i at the 0th post-Newtonian order [45], namely

$$\rho _ { \text{opt}, i } ^ { 2 } & = \frac { 4 A ^ { 2 } } { D _ { L } ^ { 2 } } \left [ F _ { +, i } ^ { 2 } ( \alpha, \delta, \psi, t ) ( 1 + \cos ^ { 2 } \iota ) + 4 F _ { \times, i } ^ { 2 } ( \alpha, \delta, \psi, t ) \cos ^ { 2 } \iota \right ] I _ { i }, \\ \text{where the amount} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$} \, \text{$1$}$$

where the amplitude A is given by

$$A = \sqrt { \frac { 5 } { 9 6 } } \left ( \frac { G \mathcal { M } _ { \text{obs} } } { c ^ { 3 } } \right ) ^ { 5 / 6 } c \pi ^ { - 2 / 3 } \,, \quad \quad ( 2 1 )$$

with M obs the redshifted chirp mass. F + , i and F × , i are the detector antenna patterns for the detector i , for a GW signal located at right ascension α and declination δ , with polarization angle ψ , and arrival at a time t . In our simulations, we fixed the antenna pattern functions at the arrival time t of the GW. This approximation is reasonable for the current groundbased detector, as the BNS signals enter the sensitivity band around 10 Hz, with only a few minutes left prior to the merger. Since the antenna patterns vary with the sidereal day, they can be considered constant for this case. However, for ET, BNS signals will enter the sensitivity band around ∼ O (1) Hz and be visible for a few days. In order to ease the computational load of our simulations, we impose a lower limit to the ET sensitivity at 10 Hz. In other words, we are artificially reducing the GW detection range that ET would have. However, as we will see later, this is not a problem for the simulation as the detection range for GRBs with Fermi-like satellites is significantly smaller than the restricted ET detection range.

The quantity I i is defined as

$$I _ { i } = \int _ { 1 0 \, { \tt H z } } ^ { f _ { i s o o } } \frac { f ^ { - 7 / 3 } } { S _ { n, i } ( f ) } d f \,,$$

where f isco ≈ 4 4 kHz ( . M ⊙ / M tot ) 3 is the orbital frequency of the innermost circular orbit of the BNS system, and S n i , ( f ) is the detector's power spectral density. For this, we use the sensitivities expected to be reached by Advanced LIGO and Virgo during O5 4 [34] and by ET [29]. Finally, we assume 100% duty cycles for all detectors.

- c. Short gamma-ray burst detection: We approach the simulation of sGRB signals following [46]. We assume that every BNS merger produces an sGRB prompted by a jet of relativistic matter ejected from the merger. We assume that the jet is narrowly collimated with an opening angle θ B around the direction perpendicular to the binary orbital plane. We model the luminosity of the sGRB as

$$L ( \iota ; \theta _ { B }, L _ { \max } ) = L _ { \max } e ^ { - \frac { \iota ^ { 2 } } { 2 \theta _ { B } ^ { 2 } } } \, \text{ \quad \ \ } ( 2 3 )$$

with θ B = 4 5 rad, compatible with the average value of the . opening angle as found by [47], and L max is drawn from a log-normal distribution with a mean of 5 × 10 51 erg s -1 and

3 M tot is the total mass of the BNS system at the observer.

4 https: // dcc.ligo.org LIGO-T2000012-v1 public / /

width of 0.56 dex [46]. The simple parametrization in Eq. (23) predicts a small luminosity if the binary is observed edge-on ι = π/ 2 ≫ θ B and a luminosity equal to L max when the binary is observed face-on ( ι = 0 rad). This model is consistent with most of the observed sGRBs with known redshift and jet opening angle [4, 48]. To compute the sGRB peak energy in the source frame Ep , we model a relation between Ep and L ( ; ι θ B , L max) in analogy to the one of long bursts [49]. Following [47], we write this correlation in the form

$$\log _ { 1 0 } \left ( \frac { E _ { p } } { 6 7 0 \, \text{keV} } \right ) = q + m \log _ { 1 0 } \left ( \frac { L } { 1 0 ^ { 5 2 } \frac { \text{erg} } { s } } \right ) \,, \quad ( 2 4 ) \quad \text{where} \quad.$$

the parameters q and m are set to 0.034 and 0.82, respectively. For the detection of the sGRB, we assume a Fermi gamma-ray burst monitor like experiment [50]: if the sGRB peak flux surpasses the detection threshold of 10 ph s cm 2 and Ep falls within the 10-1000 keV energy band, we consider it to be detected. The spectra of sGRBs typically exhibit an observer frame peak energy distribution centered between 0.5 and 1 MeV. In this work, we assume that the entire luminosity is radiated at the energy of the peak flux Ep .

The next step is to draw a value of the GW-sGRB prompt time delay ∆ t s for each BNS merger to simulate the time of arrival of the GW and sGRB signals. Unfortunately, there is currently no robust model for this quantity, although the literature generally agrees that these should have values between 10 -2 and 100 seconds [51]. Given these uncertainties, we use a simple Gaussian model to describe the distribution of prompt time time delays. Specifically, we consider three Gaussian distributions, each with a mean µ ∆ t of 1.7 seconds and standard deviations σ ∆ t of 1.7, 0.17, and 0.017 seconds, respectively. We chose three di ff erent values of σ ∆ t as we expect that, if the distribution of prompt time delay were almost a constant in nature, it would be easier to obtain a statistical estimation of the sources' redshift. On the other hand the value of µ ∆ t is set to be in agreement with the GW170817 GRB 170817A joint / observations [4].

d. Error budgets on luminosity distance and observed time delays: To represent measurement uncertainties inherent in actual observations, we generate simulated error values for the sources' luminosity distances and observed GWsGRB time delays. The measurement of the luminosity distance is provided from the GW detection. In principle, to obtain reliable distance estimates, one should use a full Bayesian analysis using simulated data and implementing state-of-theart waveform models for the GW signals. However, this approach is highly computationally demanding. For this reason, we adopt a GW likelihood [52] model able to describe statistically the typical error budgets for the GW luminosity distance. With this model, we associate a set of data xi with an 'observed luminosity distance' D L obs , namely

$$\mathcal { L } _ { \text{noise} } ( D _ { L } ^ { \text{obs} } | D _ { L } ^ { \text{true} } ) = \mathcal { N } ( D _ { L } ^ { \text{obs} } | \mu = D _ { L } ^ { \text{true} }, \sigma = 0. 2 \cdot D _ { L } ^ { \text{true} } ) \,, \ ( 2 5 ) \quad \text{figs.\ } \mathbb { A } _ { \text{one} }$$

where D L true is the true physical value of the luminosity distance. Additionally, we also assume that the observed SNR is not the true SNR in Eq. (20) but it is modified by possible noise fluctuations. To mimic this e ff ect, we generate the observed SNR ρ obs from a Gaussian distribution with a mean equal to the intrinsic SNR ρ true and a standard deviation of 1 [53] and apply the SNR selection criterion on it. The likelihood for the observed SNR is

$$\mathcal { L } _ { \text{noise} } ( \rho _ { \text{obs} } | \rho _ { \text{true} } ) = \mathcal { N } ( \rho _ { \text{obs} } | \mu = \rho _ { \text{true} }, \sigma = 1 ) \,. \quad ( 2 6 )$$

For the observed time delays ∆ t d , we sample the noisy measurement from a Gaussian distribution,

$$\mathcal { L } _ { \text{noise} } ( \Delta t _ { d } ^ { \text{obs} } | \Delta t _ { d } ^ { \text{true} } ) = \mathcal { N } ( \Delta t _ { d } ^ { \text{obs} } | \mu = \Delta t _ { d } ^ { \text{true} }, \sigma = 0. 0 5 \, s ) \,, \ \ ( 2 7 )$$

where the standard deviation is calibrated on the precision achieved for the timing of GW170817 and GRB 170817A [5]. The overall likelihood model is given by the product of Eqs. (25)-(27), the posterior estimation values of the luminosity distance and time delays are drawn sampling from this likelihood.

## B. Catalogs description

We describe the population distribution of the mock catalogs, focusing on the luminosity distance D L, the inclination angle ι , and the detected time delay ∆ t d . Figures 4 and 5 show the distributions obtained for D L, ∆ t d an cos ι for the cases of O5 and ET. For these example cases, we simulate 10 6 BNS mergers with a standard deviation of the Gaussian prompt time delay distribution set to σ ∆ t = 0 17 s. .

The distributions of detected BNS mergers feature two peaks at cos ι ∼ -1 and cos ι ∼ 1. This is the result of the dependence of the GW SNR and the sGRB intrinsic luminosity on the inclination angle ι [see Eqs. (20) and (23)], which causes the detection probability to be highest when the BNS merger orbital angular momentum is aligned or antialigned with respect to our line of sight (face-on face-away). / While GW detections span all cos ι values, no sGRBs are detected within the range cos ι ∈ ( -0 70 0 70). . , . This is a consequence of the exponential dependency of the sGRB luminosity on the inclination angle [see Eq. (23)].

We detect GWs up to luminosity distances DL ∼ 900 Mpc for O5 and DL ∼ 15000 Mpc for ET. For the O5 case, sGRBs will be detectable for almost all GW signals nearly face-on ( cos | ι | &gt; 0 70). . However, for ET, the GW detection reach will greatly surpass the sGRB detection reach and we will be able to detect sGRBs only for face-on binaries below DL ∼ 3500 Mpc. Consequently, about ∼ 44% of BNS detections in O5 are followed by an sGRB observation, and this rate drops to just ∼ 4% for ET (even though in terms of absolute numbers GW-sGRB joint detections in ET will be more abundant) [54]. An increase in luminosity distance detection range correlates with a shift and a broadening in the observed time delay ∆ t d distribution as may be clearly seen by comparing Figs. 4 and 5. This is a consequence of their shared dependency on the merger redshift, z [refer to Eqs. (1) and (2)]. The DL -∆ t d distribution could then be used to extract implicit redshift information from a set of detected time delays and luminosity distances. However, for BNS merger events detected in O5, the observed time delay distribution does not evolve

FIG. 4. O5 catalog: normalized distribution of the luminosity distance D L, cosine of the inclination angle cos ι , and observed time delay ∆ t d for Astrophysical events (blue line), GW detections (green line), and joint GW-sGRB detections (red line). The plots report the 68 3% and 90% credible intervals. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00362v2-figure-c33166ca-3fdd-435c-a39e-beb5bc8095d2.jpg)

FIG. 5. ET catalog: normalized distribution of the luminosity distance D L, cosine of the inclination angle cos ι , and observed time delay ∆ t d for Astrophysical events (blue line), GW detections (green line), and joint GW-sGRB detections (red line). The plots report the 68 3% and 90% credible intervals. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00362v2-figure-f36f505b-7b9b-4d8e-b1c6-f802ec31ddde.jpg)

significantly with luminosity distance. In such scenario, the redshift e ff ect is very small and, as we will show later, it will be di ffi cult to measure it.

## IV. RESULTS

Using the catalogs of GW and sGRB detections, we forecast the precision with which the prompt time delay distribution, speed of gravity and H 0 can be constrained thanks to time-delay cosmography. To accomplish this, we use icarogw a Python tool developed to infer the population properties of CBCs observed through GWs [40]. More precisely, we employ a Markov Chain Monte Carlo algorithm 5 , to sample the population parameters. We consider three inference scenarios:

- I. Speed of gravity and prompt time delays . We fix H 0 to the Planck value, i.e., H 0 = 67 66 km s . -1 Mpc -1 [41], and jointly infer the speed of gravity and prompt time delay distribution with a set of 3 parameters ( ˆ α , µ 0 ∆ t , σ ∆ t ).
- II. Cosmology and prompt time delays . We set ˆ 0 α = 0 (GR value), and jointly infer the cosmology ( H 0 ) and prompt time delay distributions ( µ ∆ t , σ ∆ t ).
- III. Fully agnostic . We jointly infer all four parameters ( ˆ α , 0 H 0 , µ ∆ t , σ ∆ t ).

In each scenario, we consider sets of N obs = [4 , 10 100] GW-, sGRB detections. The choice to start with N obs = 4 detections is motivated by the fact that in scenario III, where we jointly infer 4 parameters, it would be impossible to jointly constrain all parameters with less than 4 observations. In each scenario, we also consider multiple catalogs constructed for the O5 and ET sensitivities, and with various standard deviations of the prompt time delay ∆ t s distributions, namely σ ∆ t = [1 7 . , 0 17 0 017] s. . , .

The priors that we set for the population parameters are summarized in Table I. For the prompt time delay distribution, we use a uniform prior between [0 , 5] s for the mean µ ∆ t and a flat-in-log prior between [10 -3 , 10] s for the standard deviation σ ∆ t . The choice to use a flat-in-log prior within this range stems from the di ff erent σ ∆ t values used for the ∆ t s distribution models: the flat-in-log prior guarantees a uniform distribution across all these orders of magnitude. For the Hubble constant, we choose a uniform prior distribution between [30 140] km s , -1 Mpc -1 . Finally, for | ˆ α 0 | we use a uniform prior with a maximum value dependent on the detection and simulation scenario, spanning between ∼ 10 -19 and ∼ 10 -15 . The di ff erent prior range is chosen so that the posterior on ˆ α 0 is fully contained in the prior range. The values of the remaining population parameters are fixed in all scenarios. In the remainder of this work, any references to 'uncertainty' and 'precision' correspond to the median with symmetric 68.3% credible interval.

5 The emcee sampling algorithm [55] from B ilby [56, 57]

TABLE I. List of priors for the population parameters.

| Parameter                         | True Value                                                         | Prior distribution                                         | Distribution range                                         |
|-----------------------------------|--------------------------------------------------------------------|------------------------------------------------------------|------------------------------------------------------------|
| µ ∆ t σ ∆ t H 0 ˆ α 0 Ω M γ k z p | 1 . 7 s [1.7,0.17,0.017] s 67 . 66 kms - 1 Mpc - 1 0 0.310 2.7 6 2 | Uniform logUniform Uniform Uniform Fixed Fixed Fixed Fixed | [0,5] s [10 - 3 ,10] s [30 , 140] kms - 1 Mpc - 1 Variable |

FIG. 6. Joint and marginalized posterior distributions for ˆ α 0 , µ ∆ t and σ ∆ t in the case of 100 events detected by ET and the true value of σ ∆ t = 0 17 s. . The orange lines represent the true value of the parameters. The plots report the 68 3% and 90% credible intervals. . The marginal 1-dimensional distribution reports the 68 3% credible . intervals.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00362v2-figure-b08c148f-8a2f-4449-9544-dad2a4bba8fc.jpg)

## A. Scenario I: Speed of gravity and prompt time delays

Scenario I corresponds to the case in which we are interested in constraining only the speed of gravity and the prompt time delay distribution. We find that future observations will allow us to jointly constrain the speed of gravity and the prompt time delay distribution.

As an example, in Fig. 6, we present the joint posterior distribution for 100 events detected by ET. The interplay between source-frame time delay parameters and the speed of gravity is evident in the strong correlation between ˆ α 0 and µ ∆ t exhibited by the posterior distribution. Considering the expres- sion of the detected time delay in Eq. (2), a lower ˆ α 0 must be compensated by a higher source-frame time delay to maintain the observed time delay value fixed at its detected value. Figure 6 clearly shows that the GW-sGRB prompt time delay distribution and speed of gravity can be measured jointly. In Fig. 7 we summarize the precisions achievable for the population parameters ( ˆ 0 α , µ ∆ t , σ ∆ t ). In our simulations, we suc-

FIG. 7. Precision at the 68.3% credible interval for scenario I on the population parameters when we combine more and more joint GW-sGRB detections. Each row of subplots represents a deduced population parameter (ˆ 0 α , µ ∆ t , σ ∆ t ). Each subplot shows 6 curves; one for each prompt time delay model ( σ ∆ t = [1 7 . , 0 17 0 017] s) . , . and configuration of the detector network (O5, ET).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00362v2-figure-5990ad32-8164-4eb4-82d9-9872a4d53c78.jpg)

cessfully constrained the speed of gravity across all observation counts and prompt time delay models with an uncertainty ∆ ˆ α 0 ranging from ∼ 2 · 10 -16 down to ∼ 10 -19 . These constraints are up to 10 4 times better than those derived from the GW170817 measurement [4]. We also notice that the preci-

sion on ˆ α 0 depends on the GW detection range and prompt time delay model ( σ ∆ t = [1 7 . , 0 17 0 017] s) as well as the . , . number of detected GW events N obs . In Fig.7, for a fixed detector sensitivity and prompt time delay distribution model, the precision in inferring the parameters can decrease from 4 to 10 observations. This behavior arises due to statistical fluctuations, as di ff erent event realizations lead to variations in the precision of the inferred parameters. The scatter is more significant with fewer observations. In fact, with a su ffi ciently large number of events, the uncertainty is expected to scale as 1 / √ N obs , reducing the impact of these fluctuations. We find that a 'sharper' prompt time delay ∆ t s distribution, i.e., σ ∆ t is 10 (100) times smaller than σ ∆ t = 1 7 s, decreases the uncer-. tainty of an ˆ α 0 measurement by a factor ∼ 0.1 (0.035). Meanwhile, we find that using more sensitive detectors will improve the constraints on ˆ 0. Specifically, ET will be able to constrain α ˆ α 0 3 times better than O5-like sensitivities. The measurements of the prompt time delay mean µ ∆ t and standard deviation σ ∆ t are constrained with an uncertainty of ∼ 0 6 . -0 015 s and . ∼ 0 5 . -0 006 s respectively. . Similar to ˆ α 0 , the precision of these constraints depends on both N obs and the prompt time delay model scenario. However, µ ∆ t could only be constrained with a minimum of 100, 10 and 4 observations for a prompt time delay distribution with σ = 1 7 s, . 0.17 s, and 0.017 s, respectively. To constrain σ ∆ t at least 10 observations were necessary for σ ∆ t = 1 7 s, . while 4 observations su ffi ced for both 0.17 s and 0.017 s. The uncertainty of the ( µ ∆ t , σ ∆ t ) measurements varies with the 'sharpness' of the prompt time delay distribution in the same qualitative way as in the ˆ α 0 case. From Fig. 7, we conclude that the constraints on the prompt time delay population parameters are less dependent on the GWdetector sensitivities than the ˆ 0 α constraints. This is consistent with the expected error budgets from Eq. (2). Assuming Gaussian errors for ∆ t d , and a known prompt time delay distribution, we can see that the expected standard deviation on ˆ0 α is inversely proportional to the lookback time T . l Repeating the same procedure for the prompt time delay but fixing ˆ 0 α = 0, provides a precision on ∆ t s inversely proportional to 1 + z . As the lookback time T l grows faster than 1 + z , for high redshift values (higher detection ranges), the precision on ˆ 0 α will improve faster than the ones on prompt time delay distributions.

## B. Scenario II: Cosmology and prompt time delays

Scenario II is the case in which we would like to exploit time delays to measure H 0 when the source redshift is unknown.

Figure 8 shows the precision on H 0 and the prompt time delay parameters that we could be able to achieve with this technique. In general, we observe that the constraining power on H 0 increases as the detector network sensitivity improves and the prompt time delay distribution becomes 'sharper' (lower σ ∆ t /µ ∆ t ). For a GW detector network with O5-like sensitivities, H 0 can be constrained with a precision &lt; 30% when 100 sources will be detected and the prompt time delays have a narrow distribution ( σ ∆ t = 0 017 s). . If the dis-

FIG. 8. Precision at the 68.3% credible interval for scenario II on the population parameters when we combine more and more joint GW-sGRB detections. Each row of subplots represents a deduced population parameter ( H 0 , µ ∆ t , σ ∆ t ). Each subplot shows 6 curves; one for each prompt time delay model ( σ ∆ t = [1 7 . , 0 17 0 017] s) . , . and configuration of the detector network (O5, ET).

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00362v2-figure-10d837d6-cf28-4a3b-ac70-e2e85fd12d4c.jpg)

tribution of prompt time delays is wider, with an O5-like detector network, we will not be able to constrain H 0 even with 100 of these sources. Instead, for ET, which will observe GW events at higher redshifts, we find that H 0 could be measured already at ∼ 30% precision with 100 sources and a distribution of prompt time delays that is mildly spread ( σ ∆ t = 0 17 s). . In the ET case, we obtain that with 100 of these sources, it will be possible to constrain H 0 at &lt; 10% precision, which would be enough to provide useful hints on the H 0 tension. In the best case scenario, for N obs = 100, ET and σ ∆ t = 0 017 s, we obtain a value for the Hubble constant . of H 0 = 63 69 . + 2 81 . -2 71 . kms -1 Mpc -1 . For the prompt time de-

FIG. 9. Joint and marginalized posterior distributions for H 0 , µ ∆ t and σ ∆ t in the case of 100 events detected by ET and true value of σ ∆ t = 0 17 s. . The orange lines represent the true value of the parameters. The plots report the 68 3% and 90% credible intervals. The marginal . 1-dimensional distribution reports the 68 3% credible intervals. .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00362v2-figure-3eb27564-9d0d-44fa-8179-5bdc0812ece2.jpg)

lay parameters, by comparing Figs. 8 and Fig. 7, we find that the qualitative behavior for these two parameters is similar to what was found in scenario I.

In Fig. 9, we show the joint posterior distribution for 100 events detected by ET and a mildly narrow prompt time delay distribution ( σ ∆ t = 0 17 s). . The posterior distribution shows a strong correlation between the measurement of H 0 and µ ∆ t . As can be seen from Eq. (2), since lower H 0 values lead to lower redshifts, µ ∆ t must increase to keep ∆ t d fixed to its detected value. We also observe a weak correlation for the pairs ( H 0 , σ ∆ t ) and ( µ ∆ t , σ ∆ t ). In fact, a larger H 0 results in a larger spread of the BNS redshifts. Then, due to Eq. (2), an increase in H 0 also widens the spread of ∆ t d beyond the range actually observed. To reconcile this, the standard deviation of the prompt time delays σ ∆ t needs to be reduced. Then H 0 is correlated with both µ ∆ t and σ ∆ t : this induces a correlation between them in the marginalized posterior distribution.

## C. Scenario III: Fully agnostic

In this scenario, we take an agnostic approach and assume that we wish to measure prompt time delay distributions, H 0 and also the speed of gravity.

In this case, we find that it will be di ffi cult to constrain the H 0 inside the prior range used for the analysis even when considering narrow prompt time delay distributions and a hundred sources. Figure 10 displays the uncertainties on the speed of gravity and prompt time delay parameters that we can constrain in this case. The parameters ˆ α 0 , µ ∆ t and σ ∆ t are con-

FIG. 10. Precision at the 68.3% credible interval for scenario III on the population parameters when we combine more and more joint GW-sGRB detections. Each row of subplots represents a deduced population parameter (ˆ 0 α , µ ∆ t , σ ∆ t ). Each subplot shows 6 curves; one for each prompt time delay model ( σ ∆ t = [1 7 . , 0 17 0 017] s) and . , . configuration of the detector network (O5, ET). We do not display H 0 as we are unable to constrain it for such inference scenario.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00362v2-figure-448e21a2-d0e0-4eb7-b2a0-928922929ba0.jpg)

strained with uncertainty of ∼ 2 · 10 -16 -10 -18 , ∼ 0 5 . -0 02 s . and ∼ 0 5 . -0 006 s, respectively. . In comparison to scenario I (cf. Fig. 7), these parameters are less constrained. This is a consequence of the fact that in this case we are also fitting for H 0 , which introduces additional degeneracies in the estimation of prompt time delay distributions.

For H 0 , we find that, for each detector and prompt time delay model, we will not be able to measure the Hubble constant even with 100 GW detections. This is a consequence of the fact that the determination of ˆ α 0 is degenerate with the prompt time delay distribution, which is a crucial ingredient to get the implicit redshift information required for H 0 . In other

FIG. 11. Joint and marginalized posterior distributions for ˆ α 0 , H 0 , µ ∆ t and σ ∆ t in the case of 100 events detected by ET and the true value of σ ∆ t = 0 17 s. . The orange lines represent the true value of the parameters. The plots report the 68 3% and 90% credible intervals. . The marginal 1-dimensional distribution reports the 68 3% credible . intervals.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00362v2-figure-898fbe59-6944-41fc-89df-a53b49d6ebbe.jpg)

words, H 0 and the GW propagation speed are correlated with each other. The two parameters are degenerate since di ff erent combinations of their values can result in the same observed time delay. We display the correlations among the population parameters in Fig. 11, where we show the joint posterior distribution for 100 events detected by ET.

## V. CONCLUSION

In this paper, we presented a new methodology to jointly estimate the speed of gravity, the universe expansion, and prompt time delay distribution using joint GW-sGRB observations without the need for a direct redshift measurement.

In Sec. III, we described the procedure followed to generate mock catalogs of joint GW-sGRB detections. We performed extensive simulations, considering two di ff erent detection networks (O5 and ET), three prompt time delay models ∆ t s ( σ ∆ t = [1 7 . , 0 17 0 017] s), and three di . , . ff erent observing scenarios with 4, 10, and 100 joint GW-sGRB detections. In Sec. IV, we used these catalogs to jointly infer the speed of gravity, the prompt time delay distribution, and the Hubble constant.

Wefound that constraining the Hubble constant through our method requires ∼ O (100) GWs-sGRB detections that will probably be observed with next-generation detectors. For instance, with 100 events detected by ET and standard deviation σ ∆ t = 0 17 s, we constrained the 68% credible interval of . H 0

to 74 + 15 -14 kms -1 Mpc -1 . We think this is not a viable method to measure the Hubble constant for O5 as we will have too few GWs-sGRBs events.

Constraints on the speed of gravity and prompt time delay distribution can be measured even with ∼ 10 events, however. With 10 detections observed during O5 and a true value of σ ∆ t = 0 17 s, we set a constraint on . ˆ α 0 of order O (10 -17 ) (100 times better than the current measurement from GW170817) while also constraining the mean and standard deviation of the ∆ t s distribution to 2 08 . + 0 23 . -0 24 . s and 0 21 . + 0 07 . -0 05 . s precision, respectively.

Finally, in the fully agnostic scenario (joint inference of all four parameters), we were unable to constrain H 0 for any of the GW-sGRB catalogs considered in this work.

## ACKNOWLEDGMENTS

The authors acknowledge the support of the computing facilities at INFN Rome of the Amaldi Research center funded by the MIUR program 'Dipartimento di Eccellenza' (CUP: B81I18001170001). This material is based upon work supported by NSF's LIGO Laboratory which is a major facility fully funded by the National Science Foundation.

- [1] E. Di Valentino, E. Saridakis, and A. Riess, Cosmological tensions in the birthplace of the heliocentric model, NatAs 6 , 1353 (2022), arXiv:2211.05248 [astro-ph.CO].
- [2] B. F. Schutz, Determining the Hubble Constant from Gravitational Wave Observations, Nature 323 , 310 (1986).
- [3] D. E. Holz and S. A. Hughes, Using Gravitational-Wave Standard Sirens, ApJ 629 , 15 (2005), arXiv:astro-ph 0504616 / [astro-ph].
- [4] B. P. Abbott et al. , Gravitational waves and gamma-rays from a binary neutron star merger: GW170817 and GRB 170817a, The Astrophysical Journal 848 , L13 (2017).
- [5] B. P. Abbott et al. , Multi-messenger Observations of a Binary Neutron Star Merger, ApJL 848 , L12 (2017), arXiv:1710.05833 [astro-ph.HE].
- [6] B. P. Abbott et al. (LIGO Scientific, Virgo, 1M2H, Dark Energy Camera GW-E, DES, DLT40, Las Cumbres Observatory, VINROUGE, MASTER), A gravitational-wave standard siren measurement of the Hubble constant, Nature 551 , 85 (2017), arXiv:1710.05835 [astro-ph.CO].
- [7] B. P. Abbott et al. , Properties of the Binary Neutron Star Merger GW170817, PhRvX 9 , 011001 (2019), arXiv:1805.11579 [grqc].
- [8] E. Chassande-Mottin, K. Leyde, S. Mastrogiovanni, and D. A. Steer, Gravitational wave observations, distance measurement uncertainties, and cosmology, PhRvD 100 , 083514 (2019), arXiv:1906.02670 [astro-ph.CO].
- [9] R. Abbott et al. (LIGO Scientific, Virgo), GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, Phys. Rev. X 11 , 021053 (2021), arXiv:2010.14527 [gr-qc].
- [10] R. Abbott et al. (LIGO Scientific, VIRGO), GWTC-2.1: Deep Extended Catalog of Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run, Phys. Rev. D 109 , 022001 (2024), arXiv:2108.01045 [gr-qc].
- [11] R. Abbott et al. (KAGRA, VIRGO, LIGO Scientific), Population of Merging Compact Binaries Inferred Using Gravitational Waves through GWTC-3, Phys. Rev. X 13 , 011048 (2023), arXiv:2111.03634 [astro-ph.HE].
- [12] S. R. Taylor, J. R. Gair, and I. Mandel, Hubble without the Hubble: Cosmology using advanced gravitational-wave detectors alone, Phys. Rev. D 85 , 023535 (2012), arXiv:1108.5161 [grqc].
- [13] D. Wysocki, J. Lange, and R. O'Shaughnessy, Reconstructing phenomenological distributions of compact binaries via gravitational wave observations, Phys. Rev. D 100 , 043012 (2019), arXiv:1805.06442 [gr-qc].
- [14] W. M. Farr, M. Fishbach, J. Ye, and D. Holz, A Future PercentLevel Measurement of the Hubble Expansion at Redshift 0.8 With Advanced LIGO, Astrophys. J. Lett. 883 , L42 (2019), arXiv:1908.09084 [astro-ph.CO].
- [15] S. Mastrogiovanni, K. Leyde, C. Karathanasis, E. ChassandeMottin, D. A. Steer, J. Gair, A. Ghosh, R. Gray, S. Mukherjee, and S. Rinaldi, On the importance of source population models for gravitational-wave cosmology, Phys. Rev. D 104 , 062009 (2021), arXiv:2103.14663 [gr-qc].
- [16] M. Mancarella, E. Genoud-Prachex, and M. Maggiore, Cosmology and modified gravitational wave propagation from binary black hole population models, Phys. Rev. D 105 , 064030 (2022), arXiv:2112.05728 [gr-qc].
- [17] S. Mukherjee, The redshift dependence of black hole mass dis-
- tribution: is it reliable for standard sirens cosmology?, Mon. Not. Roy. Astron. Soc. 515 , 5495 (2022), arXiv:2112.10256 [astro-ph.CO].
- [18] R. Gray, Gravitational wave cosmology: measuring the Hubble constant with dark standard sirens , Ph.D. thesis, University of Glasgow (2021).
- [19] R. Gray, C. Messenger, and J. Veitch, A pixelated approach to galaxy catalogue incompleteness: improving the dark siren measurement of the Hubble constant, Mon. Not. Roy. Astron. Soc. 512 , 1127 (2022), arXiv:2111.04629 [astro-ph.CO].
- [20] J. R. Gair et al. , The Hitchhiker's guide to the galaxy catalog approach for gravitational wave cosmology, The Astronomical Journal 166 , 22 (2023), arXiv:2212.08694 [gr-qc].
- [21] K. Leyde, S. Mastrogiovanni, D. A. Steer, E. ChassandeMottin, and C. Karathanasis, Current and future constraints on cosmology and modified gravitational wave friction from binary black holes, JCAP 09 , 012, arXiv:2202.00025 [gr-qc].
- [22] J. M. Ezquiaga and D. E. Holz, Spectral Sirens: Cosmology from the Full Mass Distribution of Compact Binaries, Phys. Rev. Lett. 129 , 061102 (2022), arXiv:2202.08240 [astroph.CO].
- [23] N. Borghi, M. Mancarella, M. Moresco, M. Tagliazucchi, F. Iacovelli, A. Cimatti, and M. Maggiore, Cosmology and Astrophysics with Standard Sirens and Galaxy Catalogs in View of Future Gravitational Wave Observations, The Astrophysical Journal 964 , 191 (2024), arXiv:2312.05302 [astro-ph.CO].
- [24] C. Karathanasis, S. Mukherjee, and S. Mastrogiovanni, Binary black holes population and cosmology in new lights: signature of PISN mass and formation channel in GWTC-3, MNRAS 523 , 4539 (2023), arXiv:2204.13495 [astro-ph.CO].
- [25] G. Pierra, S. Mastrogiovanni, S. Perri` es, and M. Mapelli, A Study of Systematics on the Cosmological Inference of the Hubble Constant from Gravitational Wave Standard Sirens, Phys. Rev. D 109 , 083504 (2024), arXiv:2312.11627 [astroph.CO].
- [26] R. Abbott et al. (LIGO Scientific, Virgo,, KAGRA, VIRGO), Constraints on the Cosmic Expansion History from GWTC-3, Astrophys. J. 949 , 76 (2023), arXiv:2111.03604 [astro-ph.CO].
- [27] M. Punturo, M. Abernathy, F. Acernese, B. Allen, N. Andersson, K. Arun, F. Barone, B. Barr, M. Barsuglia, M. Beker, et al. , The Einstein Telescope: a third-generation gravitational wave observatory, CQGra 27 , 194002 (2010).
- [28] M. Maggiore et al. , Science Case for the Einstein Telescope, JCAP 03 , 050, arXiv:1912.02622 [astro-ph.CO].
- [29] M. Branchesi, M. Maggiore, D. Alonso, C. Badger, B. Banerjee, F. Beirnaert, E. Belgacem, S. Bhagwat, G. Boileau, S. Borhanian, et al. , Science with the Einstein Telescope: a comparison of di ff erent designs, JCAP 2023 (7), 068, arXiv:2303.15923 [gr-qc].
- [30] L. A. Antonelli, P. D'Avanzo, R. Perna, L. Amati, S. Covino, S. Cutini, V. D'Elia, S. Gallozzi, A. Grazian, E. Palazzi, and et al., GRB 090426: the farthest short gamma-ray burst?, A&amp;A 507 , L45 (2009), arXiv:0911.0046 [astro-ph.HE].
- [31] T. L. S. Collaboration et al. , Advanced ligo, Classical and Quantum Gravity 32 , 074001 (2015).
- [32] F. Acernese et al. , Advanced virgo: a second-generation interferometric gravitational wave detector, Classical and Quantum Gravity 32 , 024001 (2014).
- [33] T. Akutsu et al. , Overview of KAGRA: Detector design and construction history, Progress of Theoretical and Experimental Physics 2021 ,

- 05A101 (2020), https: // academic.oup.com ptep article-/ / pdf 2021 5 05A101 37974994 ptaa125.pdf. / / / / /
- [34] B. P. Abbott et al. (KAGRA, LIGO Scientific, VIRGO), Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA, Living Rev. Rel. 21 , 3 (2018), arXiv:1304.0670 [gr-qc].
- [35] S. Mastrogiovanni, D. A. Steer, and M. Barsuglia, Probing modified gravity theories and cosmology using gravitationalwaves and associated electromagnetic counterparts, PhRvD 102 , 044009 (2020), arXiv:2004.01632 [gr-qc].
- [36] A. Enea Romano and M. Sakellariadou, Constraining the time evolution of the propagation speed of gravitational waves with multimessenger astronomy, arXiv , arXiv:2309.10903 (2023), arXiv:2309.10903 [gr-qc].
- [37] A. E. Romano, E ff ective speed of gravitational waves, PhLB 851 , 138572 (2024), arXiv:2211.05760 [gr-qc].
- [38] R. Ghosh, S. Nair, L. Pathak, S. Sarkar, and A. S. Sengupta, Does the speed of gravitational waves depend on the source velocity?, Phys. Rev. D 108 , 124017 (2023).
- [39] I. Mandel, W. M. Farr, and J. R. Gair, Extracting distribution parameters from multiple uncertain observations with selection biases, Monthly Notices of the Royal Astronomical Society 486 , 1086 (2019).
- [40] S. Mastrogiovanni, G. Pierra, S. Perri` es, D. Laghi, G. Caneva Santoro, A. Ghosh, R. Gray, C. Karathanasis, and K. Leyde, ICAROGW: A python package for inference of astrophysical population properties of noisy, heterogeneous, and incomplete observations, Astron. Astrophys. 682 , A167 (2024), arXiv:2305.17973 [astro-ph.CO].
- [41] N. Aghanim et al. (Planck), Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 , A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astroph.CO].
- [42] R. Abbott et al. (KAGRA, VIRGO, LIGO Scientific), Population of Merging Compact Binaries Inferred Using Gravitational Waves through GWTC-3, Phys. Rev. X 13 , 011048 (2023), arXiv:2111.03634 [astro-ph.HE].
- [43] P. Madau and M. Dickinson, Cosmic star-formation history, Annual Review of Astronomy and Astrophysics 52 , 415 (2014).
- [44] R. Abbott et al. (LIGO Scientific, Virgo), Population Properties of Compact Objects from the Second LIGO-Virgo Gravitational-Wave Transient Catalog, Astrophys. J. Lett. 913 , L7 (2021), arXiv:2010.14533 [astro-ph.HE].
- [45] C. Cutler and E . E. Flanagan, Gravitational waves from merg-´ ing compact binaries: How accurately can one extract the binary's parameters from the inspiral waveform?, Physical Re-

view D 49 , 2658 (1994).

- [46] O. S. Salafia, G. Ghisellini, A. Pescalli, G. Ghirlanda, and F. Nappo, Structure of gamma-ray burst jets: intrinsic versus apparent properties, Monthly Notices of the Royal Astronomical Society 450 , 3549 (2015).
- [47] G. Ghirlanda et al. , Short gamma-ray bursts at the dawn of the gravitational wave era, Astron. Astrophys. 594 , A84 (2016), arXiv:1607.07875 [astro-ph.HE].
- [48] W. Fong and E. Berger, The locations of short gamma-ray bursts as evidence for compact object binary progenitors, The Astrophysical Journal 776 , 18 (2013).
- [49] P. O. Francisco J. Virgili, Bing Zhang and E. Troja, Are all short-hard gamma-ray bursts produced from mergers of compact stellar objects?, The Astrophysical Journal 727 , 109 (2011), published 2011-01-11.
- [50] D. J. Thompson and C. A. Wilson-Hodge, Fermi gamma-ray space telescope, in Handbook of X-ray and Gamma-ray Astrophysics (Springer Nature Singapore, 2022) pp. 1-31.
- [51] B. Zhang, The delay time of gravitational wave - gammaray burst associations, Frontiers of Physics 14 , 10.1007 s11467-/ 019-0913-4 (2019).
- [52] J. R. Gair, A. Ghosh, R. Gray, D. E. Holz, S. Mastrogiovanni, S. Mukherjee, A. Palmese, N. Tamanini, T. Baker, F. Beirnaert, et al. , The hitchhiker's guide to the galaxy catalog approach for dark siren gravitational-wave cosmology, The Astronomical Journal 166 , 22 (2023).
- [53] M. Fishbach, D. E. Holz, and W. M. Farr, Does the Black Hole Merger Rate Evolve with Redshift?, Astrophys. J. Lett. 863 , L41 (2018), arXiv:1805.10270 [astro-ph.HE].
- [54] S. Ronchini, M. Branchesi, G. Oganesyan, B. Banerjee, U. Dupletsa, G. Ghirlanda, J. Harms, M. Mapelli, and F. Santoliquido, Perspectives for multimessenger astronomy with the next generation of gravitational-wave detectors and high-energy satellites, A&amp;A 665 , A97 (2022), arXiv:2204.01746 [astro-ph.HE].
- [55] D. Foreman-Mackey, D. W. Hogg, D. Lang, and J. Goodman, emcee: The MCMC Hammer, PASP 125 , 306 (2013), arXiv:1202.3665 [astro-ph.IM].
- [56] G. Ashton et al. , BILBY: A user-friendly Bayesian inference library for gravitational-wave astronomy, Astrophys. J. Suppl. 241 , 27 (2019), arXiv:1811.02042 [astro-ph.IM].
- [57] I. M. Romero-Shaw et al. , Bayesian inference for compact binary coalescences with bilby: validation and application to the first LIGO-Virgo gravitational-wave transient catalogue, Mon. Not. Roy. Astron. Soc. 499 , 3295 (2020), arXiv:2006.00714 [astro-ph.IM].
Leonardo Iampieri Simone Mastrogiovanni Francesco Pannarale
Measuring the speed of gravity and the cosmic expansion with time delays between gravity and light from binary neutron stars
The first observation of a gravitational wave (GW) and a short gamma-ray
burst (sGRB) emitted by the same binary neutron star (BNS) merger officially
opened the field of GW multimessenger astronomy. In this paper, we define and
address $\textit{lagging sirens}$, a new class of multimessenger BNSs for which
associated GWs and sGRBs are observed without the identification of their host
galaxy. We propose a new methodology to use the observed time delay of these
sources to constrain the speed of gravity that is, the propagation speed of
gravitational waves, the Hubble constant and the prompt time delay distribution
between GWs and sGRBs, even though a direct redshift estimation from the host
galaxy is unavailable. Our method exploits the intrinsic relation between GWs
and sGRBs observed and prompt time delays to obtain a statistical redshift
measure for the cosmological sources. We show that this technique can be used
to infer the Hubble constant at the $10\%$~level of precision with
future-generation GW detectors such as the Einstein Telescope and only 100
observations of this kind. The novel procedure that we propose has systematics
that differ completely from the ones of previous GW methods for cosmology.
Additionally, we demonstrate for the first time that the speed of gravity and
the distribution of the prompt time delays between GWs and sGRBs can be
inferred conjointly with less than 10 sources even with current GW detector
sensitivities.