Issue 
A&A
Volume 652, August 2021



Article Number  A1  
Number of page(s)  13  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/202040229  
Published online  30 July 2021 
The potential role of binary neutron star merger afterglows in multimessenger cosmology
^{1}
Université de Paris, CNRS, Astroparticule et Cosmologie, 75006 Paris, France
email: mastrosi@apc.in2p3.fr
^{2}
Sorbonne Université, CNRS, Institut d’Astrophysique de Paris, 75014 Paris, France
email: duque@iap.fr
Received:
23
December
2020
Accepted:
13
May
2021
Binary neutron star mergers offer a new and independent means of measuring the Hubble constant H_{0} by combining the gravitationalwave inferred source luminosity distance with its redshift obtained from electromagnetic followup. This method is limited by the intrinsic degeneracy between the system distance and orbital inclination in the gravitationalwave signal. Observing the afterglow counterpart to a merger can further constrain the inclination angle, allowing this degeneracy to be partially lifted and improving the measurement of H_{0}. In the case of the binary neutron star merger GW170817, afterglow lightcurve and imaging modeling thus allowed the H_{0} measurement to be improved by a factor of three. However, systematic access to afterglow data is far from guaranteed. In fact, though each one allows a leap in H_{0} precision, these afterglow counterparts should prove rare in forthcoming multimessenger campaigns. We combine models for emission and detection of gravitationalwave and electromagnetic radiation from binary neutron star mergers with realistic population models and estimates for afterglow inclination angle constraints. Using these models, we quantify how fast H_{0} will be narrowed down by successive multimessenger events with and without the afterglow. We find that because of its rareness and though it greatly refines angle estimates, the afterglow counterpart should not significantly contribute to the measurement of H_{0} in the long run.
Key words: gravitational waves / cosmological parameters / stars: neutron / Kuiper belt: general
© S. Mastrogiovanni et al. 2021
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
The detection of gravitational waves (GWs) from compact binary coalescence (Abbott et al. 2019) has opened a new window to study the Universe. Gravitationalwave sources are a new type of “standard candle”, usually referred to as “standard sirens” as it is possible to directly measure their luminosity distance (d_{L}) from the GW signal (Schutz 1986). Hence, if supplied with the source redshift information, GW detections can be used to measure cosmological parameters (Holz & Hughes 2005; Nissanke et al. 2013; Chen et al. 2018; Mortlock et al. 2019), such as the Hubble constant H_{0}. This possibility is of great interest given the current tension between the H_{0} measurement at early and late epochs of the Universe (Freedman 2017).
The first GW measurement of H_{0} was made possible by the multimessenger observation of the binary neutron star (BNS) merger GW170817 (Abbott et al. 2017a) and its associated kilonova, which enabled the identification of the host galaxy and its redshift, leading to a new and independent measurement of km s^{−1} Mpc^{−1} (Abbott et al. 2017b). Currently, GW170817 is the only GW event observed with an electromagnetic counterpart that allowed this kind of cosmological measure. In the absence of an electromagnetic counterpart, one can use the redshifts of all cataloged galaxies with positions consistent with the 3D GW skymap to measure H_{0} (the socalled “dark siren” method; Fishbach et al. 2019; SoaresSantos et al. 2019; Gray et al. 2020), leverage tidal effects in the GW waveform to estimate the source redshift (Messenger & Read 2012; Del Pozzo et al. 2017), or exploit the power spectrum of GWs and galaxy distributions (Mukherjee et al. 2021b, 2020).
These methods encounter challenges due to galaxy catalog incompleteness and the difficulty in detecting tidal effects in the GW signal, as well as uncertainty in the equation of state of neutron star matter. Overall, an electromagnetic counterpart detection remains the best prospect for obtaining precise H_{0} measurements with GWs.
The H_{0} measurement is limited by the large uncertainty on the GW luminosity distance. This arises because the GW emission is not isotropic and the luminosity distance determination is degenerate with the binary orbital inclination to the lineofsight, hereafter referred to as ι (Chen et al. 2019). To improve the H_{0} measurement, it is crucial to break this degeneracy.
One solution is to measure ι from anisotropic electromagnetic signals emitted after the merger. The most stringent constraint comes from the photometry and very long baseline interferometry (VLBI) imaging of the afterglow. This is nonthermal radiation from the shock formed as the relativistic merger outflow decelerates in the circummerger environment. For GW170817, these observations provided inclination angle measurements precise to ∼12 deg (e.g., Guidorzi et al. 2017; Troja et al. 2019; Hajela et al. 2019) and ∼5 deg, respectively (Mooley et al. 2018a; Coughlin et al. 2019a; Ghirlanda et al. 2019; Ascenzi et al. 2021), and refined the H_{0} measurement to km s^{−1} Mpc^{−1} (Hotokezaka et al. 2019).
Further afterglow counterparts could drastically improve H_{0} measurements and play a leading role in multimessenger cosmology. However, they are faint and difficult to detect for distant or very inclined binaries and should therefore prove rare in the future (Saleem et al. 2018; Gottlieb et al. 2019; Duque et al. 2019). Furthermore, using electromagnetic measurements of ι entails selection effects that, if not correctly understood, can significantly bias the subsequent measurement of H_{0} and must be carefully taken into account (Chen 2020).
In this prospective study, we consider both the likeliness of the afterglow to be observed and its capability to constrain ι to realistically quantify its benefit for the measurement of the Hubble constant with future multimessenger events. We will assume that the peculiar motion of the galaxies is accurately known, as their mismatch could introduce a systematic bias.
Afterglowderived information on ι is much more precise than that from GW data only. Nonetheless, we find that, for advanced GW detectors in their O3 and O4+ configurations, detecting afterglows should be so rare that, statistically, afterglows will not accelerate the narrowingdown of H_{0}.
This paper is organized as follows. In Sect. 2 we describe the possible multimessenger observation scenarios of a GW event with electromagnetic counterparts following which a measurement of H_{0} can be made. In Sect. 3 we evaluate and quantify the likeliness of the different scenarios by combining models for the emission of GWs and electromagnetic radiation from BNSs, models for the detection of these signals, and population models for the sources. In Sect. 4, we briefly describe the prospects of joint GW and electromagnetic detections in these different scenarios. In Sect. 5, we evaluate the potential of electromagnetic counterparts in contributing to the measurement of H_{0} by combining their capability to measure the inclination angle ι with their likelihood of being observed in the scenarios described above.
2. Inferring H_{0} with gravitational waves and electromagnetic counterparts
In this section, we review the possible multimessenger observation scenarios expected after a GW event from a BNS merger.
Until the advent of deep radio surveys such as the Square Kilometer Array (Dewdney et al. 2009), the kilonova will continue to be the only electromagnetic counterpart to GW events that can lead to an identification of the system host galaxy and redshift. For a kilonova to be detected during followup of a GW event, it must be in reach of the followup telescopes in terms of both magnitude and sky position, and the source must be localized sufficiently well to be discovered in time before it fades. Magnitude and localization are, however, not the only conditions for detection. As recent searches for kilonova signals both during GW followup campaigns (Kasliwal et al. 2020) and in archival data (Andreoni et al. 2020) have shown, difficulties can arise in recognizing kilonovae among a myriad of optical transients, even with quality spectroscopic or color evolution observations.
Detection of the kilonova counterpart and thereby acquiring the system’s redshift through its host galaxy is the minimal scenario required for a multimessenger measurement of H_{0}. We refer to this scenario as “Level 1”. It is subject to detection criteria concerning the GW and kilonova signals, described in Sect. 3.
In this scenario, the information on d_{L} is provided solely by the GW data, without any contribution from electromagnetic counterparts. Indeed, any direct distance information from an electromagnetic counterpart would require using the cosmic distance ladder, which is of course excluded in the perspective of measuring H_{0}.
In principle, the kilonova signal could indirectly contribute through the measurement of ι from colorevolution considerations (Kashyap et al. 2019; Dhawan et al. 2020). However, currently, these signals crucially lack modeling and observing history. Inclination angle measurements from kilonova data are very model dependent and lack robustness (Doctor 2020; Heinzel et al. 2021). Only when the kilonova sample grows will the potential impact of kilonovaderived angle constraints be appreciated. Therefore, we excluded any contribution to ι from kilonova data.
As expected from the observation of short gammaray bursts and evidenced in the case of GW170817, relativistic jets are launched from BNS mergers (Mooley et al. 2018a; Ghirlanda et al. 2019). The interaction of this jet with the circummerger environment leads to longlived^{1} nonthermal emission in the radio to the Xray bands: the afterglow counterpart.
The afterglow photometry can provide an independent measurement of inclination angle. In such measurements from offaxis jet afterglow data fitting, there is unavoidable degeneracy with the jet opening angle (Nakar & Piran 2021). However, when combined with prior estimates for short gammaray burst jet opening angles (see references in Sect. 3.2.3), afterglow data can still lead to inclination angle information. It can therefore indirectly inform on d_{L} by breaking the d_{L}/ι degeneracy in the GW data. We refer to the scenario where Level 1 is realized and a ι measurement from the jet afterglow photometry is made as “Level 2”. It is subject to the realization of the GW and kilonova detection criteria and to detection criteria on the jet afterglow light curve, which we detail in the next section.
Additionally, the relativistic nature of the jetted outflow can be revealed by VLBI observations (Mooley et al. 2018a; Ghirlanda et al. 2019) that evidence an apparent superluminal motion of the jet head. Detecting this centroid displacement is possible for events that are particularly close or bright or under specific inclination angle conditions (Duque et al. 2019; Dobie et al. 2020). Doing so further constrains ι and narrows down the measurement of H_{0}, as shown in the case of GW170817. We refer to the scenario where such a constraint on ι can be extracted from afterglow VLBI imaging, in addition to those of afterglow photometry, as “Level 3”. This final level of ι constraint is the most informative on ι and H_{0}, but also the most difficult to obtain.
Other electromagnetic counterparts could potentially provide further independent measurements of ι. These are the short gammaray burst and the rebrightening in the source’s multiwavelength signal due to the emergence of emission from the front shock of the decelerating mildly relativistic ejecta responsible for the kilonova emission, called the “kilonova afterglow” (Hotokezaka et al. 2018; Kathirgamaraju et al. 2019). However, the short gammaray burst should prove extremely rare in future events (Ghirlanda et al. 2016; Beniamini et al. 2019) and there is a lack of robust modeling for gammaray signals, especially for very inclined events. Furthermore, the singularity of GRB170817A with respect to other short gammaray bursts has cast more uncertainty on the modeling of gammaray emission from BNS mergers (Kasliwal et al. 2017; Nakar et al. 2018). Modeling of the kilonova afterglow is still uncertain, rendering any angle measurement difficult (Duque et al. 2020). To sum up, we consider neither the gammaray burst nor the kilonova afterglow as viable to measure ι and thus H_{0}.
3. Evaluation of detection probabilities and selection effects
Here, we evaluate the likelihoods of the different observing scenarios to occur, by computing the probabilities of detection of the GW and electromagnetic signals.
3.1. Gravitationalwave detection probability
In GW searches an event is detected if its detected signaltonoise ratio (S/N) exceeds a certain threshold. The detected S/N is a measure of the power of the GW signal registered at the GW detectors. The higher the S/N, the more chance there is to recognize the signal from detector noise. It is important to note that differs from the optimalfilter S/N ρ_{opt} that is calculated taking into account the average sensitivity of the detector network. Indeed, noise fluctuations are not included in optimalfilter S/N. Also, due to the same noise fluctuations, the GW event will be detected with values and that differ from the true d_{L} and cos ι. The detected , and the corresponding can be sampled from the GW likelihood model.
In this article, we used the Cutler and Flanagan (CF) approximant for the GW likelihood (Cutler & Flanagan 1994; Poisson & Will 1995; ChassandeMottin et al. 2019). This approximate can reproduce the GW likelihood in the and space under the assumption that the chirp mass of the signal is well estimated, as is always the case for BNS detections (Cutler & Flanagan 1994). By sampling from the CF approximate, we can obtain a value of the detected and , which can then be used to compute the detected S/N with the following (Cutler & Flanagan 1994; ChassandeMottin et al. 2019):
where , while ϵ_{d}, σ_{d} are variables that depend on the detector network and the skyposition of the GW source and is the GW polarization angle. The variable ρ_{fo} is the optimalfilter S/N that the binary would have had if it had been faceon,
with S_{n, aver}(f) the harmonic mean of the noise power spectrum densities (PSD) of the interferometric detectors composing the network, G and c universal gravitational and speedoflight constants and ℳ_{c} the system’s chirp mass in the detector frame. As per the low redshift range of the systems in this study (see below), we assimilate the chirp mass and last stable orbit frequency f_{LSO} in the detector frame to their values in the source frame. The integral low boundary f_{low} = 20 Hz is set to the lowfrequency cutoff for groundbased GW detectors.
To evaluate the GW detection probability, we perform a Monte Carlo simulation with synthetic signals. We divide the cos ι range into 200 bins and we simulate 10 000 BNS merging at fixed d_{L} and uniformly distributed over the celestial sphere. The BNS masses are generated from a Gaussian distribution with a mean of 1.35 M_{⊙} and a standard deviation of 0.15 M_{⊙} (Farrow et al. 2019). For each binary, we then draw a detected and from the CF likelihood and calculate following Eqs. (1) and (2). We then count the binaries with S/N exceeding the online matchfiltering threshold for detection of 14^{2} and compute the detection probability as the fraction of BNS events detected. The previous procedure is repeated in d_{L} with a step of 1 Mpc until only 1 of the 10 000 simulated systems is detected. Following Chen et al. (2021), we define the 0.02% “response distance” at which 0.2% of the simulated binaries (with isotropic distribution in the sky and orientation) will be detected by the network. We verified by running the simulation several times that the GW detection probability estimation is not prone to Monte Carlo statistical fluctuations.
For this paper, we use three PSDs for the HLV network, composed of the LIGO Hanford, LIGO Livingston, and Virgo detectors. The first PSD is indicative of the detectors’ sensitivity during O2 (Abbott et al. 2017c), the second is indicative of the sensitivity reached during the first 3 months of O3 (Abbott et al. 2020a), while the third one is a projection for the O4 run sensitivity (Biwer et al. 2019).
Figure 1 shows the GW detection probability as a function of the BNS d_{L} and cos ι marginalized over the GW polarization angle and skyposition. The systems we consider have redshifts ≤0.06 for both Planck (Planck Collaboration VI 2020) and SH0ES (Riess et al. 2019) values of H_{0}. Because of the low redshift range of our simulation, the shape of the GW detection probability functions are not affected by the nonlinear d_{L} − z relationship and are therefore the same for all the PSD hypotheses. As it can be seen from the figure, for all three GW sensitivity levels, faceon binaries are easier to detect as the GW emission is stronger perpendicularly to the orbital plane, and they can be observed at higher luminosity distances.
Fig. 1. GW detection probability as a function of the BNS luminosity distance and inclination for the detector network HLV using a S/N threshold of 14. The horizontal axis is scaled to the BNS 0.2% response distance , which is Mpc for O2; Mpc for O3; Mpc for O4. 
3.2. Electromagnetic counterpart detection probability
3.2.1. Level 1 observing scenario
In the perspective of upcoming highcadence and largefieldofview optical facilities such as the Zwicky Transient Facility (ZTF, Bellm et al. 2019) or Vera C. Rubin Observatory (LSST, Ivezic et al. 2008), we assume the detection of the kilonova is limited only by the kilonova magnitude and skyposition, and not by the size of the GWprovided skymap that the followup network must cover in its searches. Indeed, we consider that these survey facilities can cover all the sky available to them within the first nights of the search, that is, before the estimated time for significant dimming of kilonova signals.
While this can be considered a “bestcase” assumption, this level of performance was reached during the campaigns following GW events in the recent O3 run of the LIGOVirgo Collaboration by, for example, the GROWTH Collaboration (Kasliwal et al. 2020). In particular, for the only confirmed BNS event of the O3 run GW190425 (Coughlin et al. 2019b), the ZTF covered the ∼8000 deg^{2} of the skymap overlapping with their night sky.
Similarly to the evaluation of GW detection probability, we place systems uniformly in the sky and uniformly in cos ι and d_{L} out to 600 Mpc. This maximum distance ensures that a fraction of less than 10^{−4} of kilonovae were detected in the furthest distance bin. Our model for the inclinationangledependent peak magnitude of the kilonova signal is an empirical fit to model “W2” in Wollaeger et al. (2018, their Fig. 19). It appears that the peak magnitudes are approximately linear in cos ι up to a maximum angle after which the dependence on angle strongly decreases. We therefore set the peak absolute magnitude of these events as:
with numerical values reported in Table 1. We note that a similar dependence of kilonova magnitudes on ι was considered in Villar et al. (2017, their “asymmetric model”) and found in Kawaguchi et al. (2020). We select the detectable events as those satisfying the magnitude threshold criterion g < g_{lim} or r < r_{lim}, with lowercase letters denoting the apparent magnitudes.
Description and numerical values of various constants used in the electromagnetic emission and detection models.
Among all the events with large enough flux, followup can only detect those in its accessible sky. For observatories at latitudes of ∼33 deg such as the ZTF (northern hemisphere) and the LSST (southern hemisphere), this represents a seasonaveraged fraction p_{night} ∼ 52% of the whole sky (Bellm 2016). Thus, we deem detected the remaining events with a cointoss with probability p_{night} = 0.52.
Forming the ratio with the original number of events allows us to evaluate the detection probability , the probability of acquiring an event redshift as a function of its luminosity distance and inclination angle. This is plotted in Fig. 2. Here, the detection probabilities were normalized to p_{night}. As the p_{det} are even in ι, we represent only the cos ι > 0 range. Following what was developed in Sect. 2 and given the independence of the GW and kilonova detection processes, we define the probability of the Level 1 scenario occurring as .
Fig. 2. Probabilities of detection for the electromagnetic counterparts considered in the study: (top), (bottomleft) and (bottomright) under the two different hypotheses for the jet’s kinetic energy distribution denoted G16 and WP15. See Sect. 3.2.2 for details on these hypotheses. 
3.2.2. Level 2 observing scenario
For events where emission from the core jet dominates the afterglow radiation, the afterglow light curves are expected to display a single peak occurring when the jet has decelerated to a Lorentz factor of Γ ∼ 1/ι (e.g., Gill & Granot 2018; Mooley et al. 2018b, but see Beniamini et al. 2020 for a study of multipeaked afterglows). This peak can occur up to hundreds of days after the merger (Duque et al. 2019). Assuming instrument availability and longterm followup efforts, detecting the afterglow is simply a matter of flux sensitivity, once the source’s position is settled by the kilonova. However, detecting the afterglow at its peak does not suffice to make an inclination angle measurement. This requires an extended and wellsampled light curve on which to fit afterglow models, as was extensively done for GW170817 (e.g., Lamb & Kobayashi 2017; Resmi et al. 2018; Lazzati et al. 2018; Troja et al. 2019). We therefore define the criterion for ι measurement with the afterglow light curve as:
where F_{p} is the peak flux of the afterglow light curve, and F_{lim} is the limiting sensitivity of the followup facility. In this study, we consider the 3 GHz band and the Very Large Array (VLA) as the limiting radio facility, with F_{lim} = 15 μJy.
In addition to ι and d_{L}, the peak flux F_{p} of every event depends on the jet’s kinetic energy, the particle density in the surrounding medium, and on the microphysical parameters of the front shock formed by the decelerating jet, such as the spectral index of the shockaccelerated population of electrons, denoted p. The analytical form for F_{p} as a function of these parameters that we use can be found in Nakar et al. (2002). There is some uncertainty in the distributions of these parameters in the population of jets from BNS mergers, in particular for the jet kinetic energy. To establish , the probability of making the inclination measurement allowed by the afterglow photometric data in addition to acquiring the redshift, we use the same distribution of parameters as the population model of Duque et al. (2019). In particular, we use two distinct hypotheses regarding the distribution of the jets’ kinetic energies. These are labeled G16 and WP15 in the sequel and are derived by assuming a constant conversion factor between the luminosity of a short gammaray burst and the postburst jet kinetic energy, and starting from the short gammaray bursts luminosity functions found by Ghirlanda et al. (2016) and Wanderman & Piran (2015), respectively. Among published short gammaray burst luminosity functions, these two represent extremes in the steepness of the luminosity function, with G16 predicting many more bright bursts (and therefore many more energetic afterglows) than WP15.
By applying this population model to those events selected to establish and applying the detection criterion in Eq. (4), we calculate . It is plotted in Fig. 2. We define the likelihood that a Level 2 scenarios occurs as . As expected, the range of parameter space allowed by this scenario is much smaller than for a kilonovaonly event. In particular, the distant or inclined events are largely cut off because (for a shockaccelerated population spectral index of p = 2.2, Nakar et al. 2002), whereas r, g ∼ 1 − cos ι + log d_{L}.
We note that GW170817 does not exactly qualify for our Level 2 scenario. Indeed, it had log F_{p}/F_{lim} ∼ 0.9 < 1 (Mooley et al. 2018b). It seems that much of the uncertainty in the measurement of ι with GW170817 is held in the very early phases of the afterglow, where the fitting models most diverge (Ghirlanda et al. 2019). Had the afterglow been brighter–at the level of our Level 2 scenario–and these earlier points observed, a better measurement of ι would have certainly ensued. Nonetheless, as we detail in Sect. 5, we consider the case of GW170817 as representative of the ι measurements possible in the Level 2 scenario.
3.2.3. Level 3 observing scenario
The Level 3 scenario occurs if the motion of the jet can be measured by VLBI imaging, in addition to the measurements of afterglow photometry. Once again this assumes constant instrument availability and longterm followup. We estimate the angular displacement of the jet centroid as:
where δt_{VLBI} is the total time the afterglow remains detectable by the radiotelescope array, that is, the time its flux is above F_{lim}^{3}, and where dθ/dt_{max} is the proper motion of the remnant at the time of the afterglow peak, that is when Γ × ι ∼ 1. At this time, it is straightforward to estimate the remnant’s proper motion as:
with the apparent velocity of the remnant, and Γ_{p} = 1/ι the jet head Lorentz factor at the afterglow peak. By considering the source’s proper motion to be that at afterglow peak during the entire followup, Eq. (6) in fact overestimates δθ_{r}. Our Level 3 events are therefore likely treated in an optimistic manner, as we discuss in Sect. 6.
Also, Eq. (6) is not valid if the observer is within the jet’s opening, when ι < θ_{j} with θ_{j} the halfopening angle of the jet. In this case, no jet displacement is observed, δθ_{r} = 0. We consider θ_{j} = 0.1 ∼ 6 deg, in line with measurements on GW170817 (Gill & Granot 2018; Troja et al. 2019; Mooley et al. 2018b) and other short gammaray burst studies (Fong et al. 2015; Beniamini et al. 2019).
The afterglow flux will remain in reach of the radio network for a duration that depends on the details of its light curve and therefore on, for example, the jet structure, its expansion dynamics, and the surrounding medium density profile. To simplify, we assume all jets launched from mergers have the same structure as GW170817 and all mergers occur in a rarefied medium with a constantdensity profile as suggested by short gammaray burst observations. In this case, the slopes of the increasing and decreasing phases of the light curve are the same as for GW170817’s afterglow, regardless of inclination angle (Beniamini et al. 2020). Therefore, we empirically modeled the afterglow light curves as a broken powerlaw with slopes +0.80 and −2.2 (Mooley et al. 2018b) respectively before and after a peak occurring at a time T_{p} and a flux F_{p}. The time T_{p} depends on the same parameters as F_{p}, and relevant equations can be found in Nakar et al. (2002). The time the signal is above the radio threshold can thus be analytically estimated, and thereby the total source displacement. Then, detection of this displacement is simply conditioned by the VLBI array angular resolution: δθ_{r} > Δθ_{VLBI} ∼ 2 mas (Ghirlanda et al. 2019).
As for the other electromagnetic counterparts, we determine the probability of detecting the remnant proper motion and making the corresponding angle measurement in addition to the other measurements denoted and define . This is plotted in Fig. 2, where one can see further suppression of events.
3.3. Selection effects in multimessenger cosmology
In Appendix A, we make estimates of the selection effects impacting the measurements of H_{0} with GW and electromagnetic data.
In the circumstances of measuring H_{0} with data from GW170817, that is, a Level 3 scenario during the O2 run, the selection effects are less than 2.0% over the 60−80 km s^{−1} Mpc^{−1} range, according to Fig. A.1. This is well below the 14% precision claimed by studies making this measurement (Hotokezaka et al. 2019). We conclude that no selection effects significantly impacted the measurement of H_{0} with the GW170817 data.
However, the selection effect is not negligible compared to the 4% precision required to resolve the H_{0} tension (Freedman 2017; Feeney et al. 2018), particularly with events detected in O3 and O4type runs, where selection effects reach 2% and 4% in the 60−80 km s^{−1} Mpc^{−1} range, respectively. Then, careful consideration of selection effects becomes necessary.
4. Prospects of joint electromagnetic and gravitationalwave detections
In this section, we use the detection probability models to generate a population of events detected jointly in GW and the different electromagnetic counterparts. This allows us to study the dominant effects in the multimessenger detection process.
4.1. Simulation description
To evaluate the multimessenger capability to measure H_{0}, we simulate BNS events in a Universe described by flat ΛCDM cosmology and with H_{0} = 70 km s^{−1} Mpc^{−1} and Ω_{m} = 0.308. We generate 80 000 BNS mergers uniformly distributed in the sky up to a luminosity distance of for a given detector network with BNS 0.2% response distance (see Fig. 1). We assume a system formation rate uniform in redshift and the BNS masses are generated from a Gaussian distribution with a mean of 1.35 M_{⊙} and standard deviation 0.15 M_{⊙} (Farrow et al. 2019), cos ι is generated uniform on the unitary sphere and the GW polarization angle is distributed uniformly in the range ∈[0, π].
For each BNS merger we calculate the measured GW S/N as described in Sect. 3.1 by drawing a measured d_{L} and cos ι from the CF approximation. If the GW detected S/N exceeds a threshold of 14, then we assume that the event has been detected by the HLV network. Each GW event detected is passed to the electromagnetic followup chain described in Sect. 3.2. We use the different as described previously to decide which electromagnetic counterpart is detected according to the value of d_{L} and cos ι. At the end of the simulation, each detected GW event is associated with a flag describing the corresponding scenario: either not detected in the electromagnetic domain or detected with Levels 1, 2, or 3. We repeat this process for O2, O3 and O4type GW sensitivities and for the G16 and WP15 electromagnetic counterpart population models, for a total of six GW and population model combinations.
4.2. Rates of electromagnetic counterpart detections
In Table 2, we show the fraction of GW detections with different electromagnetic information levels predicted for O2, O3 and O4like runs. As the population models G16 and WP15 are extremes in terms of afterglow luminosities, one can consider the figures in Table 2 as confidence interval bounds for the corresponding fractions^{4}.
Average fraction of GW events observed with different electromagnetic counterpart levels.
In an O2type run, 52% of the GW events are expected to have a detectable kilonova counterpart. Since we assumed the optical instruments can cover a fraction p_{night} = 52% of the sky, this means that the kilonova magnitude was not a limiting effect for multimessenger events during O2. In other words, the GW detection probability dominates the multimessenger detection probability, as was predictable from Figs. 1 and 2.
Conversely, the fraction of events expected to have been detected in O2 with Levels 2 and 3 are between 4% and 12% and between 1% and 7%, respectively. This is from 10% to 25% and from 2% to 15% of kilonovae with detectable afterglow light curve and remnant proper motion, respectively. Therefore, as of O2, the sample of events with afterglow counterpart was limited by selection. As it can be seen from Figs. 1 and 2, this is mostly because electromagnetic signals are much more sensitive to binary inclination than GW signals.
GW170817 was a Level 3 event during O2 (Abbott et al. 2017a; Hotokezaka et al. 2019). In light of our study, the probability of this to occur was between 1% and 7%, and GW170817 was a very lucky event.
In O3type runs, the kilonova counterparts are still largely in reach of followup instruments. However, the probability of measuring cos ι with afterglow counterparts is strongly reduced compared to an O2type run, meaning the multimessenger detection process is strongly dominated by electromagnetic selection. Nonetheless, there is a nonnegligible probability of an afterglow counterpart to a GW event. In the first 6 months of O3, GW190425 (Abbott et al. 2020b) was the only plausible BNS event observed. For this event, there was no electromagnetic counterpart reported by large fieldofview facilities (Coughlin et al. 2020a). Let us note that, as opposed to the events considered in this simulation, GW190425 was detected by only two detectors (LIGO Livingston and Virgo), thus producing a large skylocalization area covering 10 200 deg^{2}.
From our simulations, we estimate that maintaining an O2like fraction of radiodetectable afterglows during an O3like GW run would require a factor of five increase in the radio sensitivity. That is, a F_{lim} on the order of 3 μJy, as projected for the Square Kilometer Array 1 “Mid” band (Dewdney et al. 2009).
In an O4like run, GW detectors detect a large fraction of BNS mergers at higher distances. In this case, the effects of magnitude limitations start to kick in, with only half of the kilonovae with sufficient flux. Taking into account skyposition limitations, this number is one out of 4. In O4type runs, it should prove extremely rare to obtain an electromagnetic cos ι measurement (Level 2 and Level 3).
The different selection biases introduced by the GW and electromagnetic detections processes when going from an O2like to an O4like run are depicted in Figs. 3 and 4. The two figures show the distributions in d_{L}, cos ι, and optimalfilter S/N detected in different levels. In an O2like run, the GW and electromagnetic detections roughly probe the same BNS mergers, with the exception of some regions in cos ι that forbid a detection by VLBI. On the other hand, in an O4like run, the electromagnetic and GW detection clearly corresponds to a different subpopulation of sources. In fact, the electromagnetic facilities can observe only closeby events–which thus have a high GW S/N–corresponding to a very small fraction of GWdetected events.
Fig. 3. Histograms for the distributions of inclination angle ι, luminosity distance d_{L} and optimalfilter GW signaltonoise ratio (see Eq. (1)) for all simulated BNS mergers (black line), those detected with GW only (red), and those detected with Level 1 (green line) or Level 3 (blue line) for an O2type run. 
In Table 3 we report the average number of years of observation required to detect BNS mergers with different levels of EM counterparts. We assume a BNS merger rate of Gpc^{−3} yr^{−1} as estimated after the first six months of O3 (The LIGO Scientific Collaboration 2020). We note that this is overestimated as this is based on a conservative choice for the detection threshold. For this simulation we selected confident BNS detections with S/N > 14 (∼8 in each detector) while current GW detectors can detect BNS with S/N ∼ 8 (Abbott et al. 2020a). Choosing a detection threshold of S/N ∼ 8 corresponds to a surveyed volume ∼6 times larger, and thus an observation time to detection reduced by the same factor. Table 2 reports the relative fractions of events detected with different levels of EM emission. From the results in Table 3 we conclude that the probability of detecting a BNS merger during a oneyear O2like run is small, and especially so for a loud event such as GW170817.
Average observation time (years) needed to detect one BNS merger with different EM counterpart levels.
5. Forecast on multimessenger H_{0} measurements
5.1. Method
We now set to quantify the benefit of electromagneticprovided information in measuring H_{0}. We use the statistical framework described in Appendix A with the posterior as in Eq. (A.1). The likelihood p(dH_{0}) encodes the statistical uncertainties of the GW and electromagnetic data x_{GW/EM} measurements:
where p_{pop}(z, cos ιH_{0}) is the assumed population distribution in redshift and inclination for the entire BNS population.
The function p_{GW}(x_{GW}d_{L}(H_{0}, z), cos ι) is the GW likelihood, which provides the GWdetected source parameters distribution for given true parameters. We take the CF approximation for consistency with our earlier computation of selection effects.
Likewise, p_{EM}(x_{EM}cos ι, z) is the likelihood for electromagnetic measurement of source parameters cos ι and z. This is informative on z only for a Level 1 scenario and on both z and cos ι for higherlevel scenarios. It is not informative on d_{L}. We decompose this electromagnetic likelihood as p_{EM}(x_{EM}z)p_{EM}(x_{EM}cos ι) by supposing the electromagnetic measurements of redshift and angle are independent. This is reasonable as the redshift information is deduced from the host galaxy alone while the ι information is expected to be provided by the jet itself.
We assume the redshift measurement is unbiased and set to a Gaussian distribution centered on the true event redshift with standard deviation 5 × 10^{−4} for all our scenarios. This is the same accuracy measured for GW170817’s redshift (Abbott et al. 2017b). For Level 1 scenarios, the electromagnetic counterpart is uninformative on cos ι and we set p_{EM}(x_{EM}cos ι) to a flat function.
For Level 2 scenarios, we assumed that one can obtain an unbiased ι constraint at the level of that deduced from GW170817’s afterglow light curve alone. That is, a Gaussian constraint with a 12 deg 1σ uncertainty (Troja et al. 2019; Hajela et al. 2019) for .
For Level 3 scenarios, we based our predicted constraints on those of GW170817 and set to an unbiased Gaussian constraint with a width of 4 deg (Mooley et al. 2018a; Ghirlanda et al. 2019; Hotokezaka et al. 2019) for all the events. We discuss the validity and impacts of these assumptions in Sect. 6.
Using Eqs. (A.1), (A.3) and (7), we simulated the H_{0} measurement process of 500 binary systems in all three observing scenarios and all three GW sensitivity hypotheses. For the individual measurements, we systematically assumed a prior on H_{0} uniform in [40, 120] km s^{−1} Mpc^{−1}.
We then combined the measurements of the first 100 events to emulate a thread of multimessenger events. We repeated the combining step after reordering the 500 events to reproduce different possible time orderings of the events. This allowed us to study the reconstruction of H_{0} by the multimessenger measurements, and in particular the speed of convergence.
5.2. Bulk comparison of observing scenarios
Breaking the d_{l} − cos ι degeneracy is fundamental for measuring H_{0}. Let us show this by inferring H_{0} using only Level 2 or 3 scenarios in comparison with Level 1 scenarios.
Figure 5 shows the H_{0} posteriors obtained by combining 10 BNS events in different observing scenarios. From the plot, we can see that, when the knowledge of cos ι is refined by electromagnetic observations, the estimation of H_{0} improves, as also noted in Chen et al. (2019). We can also observe that the H_{0} posterior reaches Gaussian convergence after under 10 events.
Fig. 5. Posterior on H_{0} obtained for ten BNS events observed with different level of electromagnetic information on ι in an O2type run. Left: Level 1 (no electromagnetic information on ι), center: Level 2 (electromagneticbased ι precision of 12 deg), right: Level 3 (electromagneticbased ι precision of 4 deg). 
The bottom panel of Fig. 6 shows the relative uncertainty ΔH_{0}/H_{0} with 1σ confidence intervals for the H_{0} estimation as more events are detected. The uncertainty corresponds to the different population realizations of the detected events. In Fig. 6, one can read that a single Level 3 event during an O2type run results on average in the uncertainty of 14% on the estimation of the Hubble constant, as observed for GW170817 (Hotokezaka et al. 2019).
Fig. 6. Convergence of H_{0} with number of detected events. Bottom: precision on the estimation of H_{0} with 68.3% confidence intervals as a function of the number of events detected at each EM level. Here, all events are assumed to be at a given electromagnetic information level, as denoted in different colors. We also indicate the tension in the Hubble constant and the precision on H_{0} obtained with GW170817 and counterparts. Top: expected number of years of continuous observation required to detect events in various multimessenger scenarios. 
It is clear that (i) the precision on H_{0} improves as more events are combined, (ii) the convergence is faster when cos ι is more constrained from the electromagnetic emission, (iii) above about ten events (Chen et al. 2018), the combined H_{0} posterior becomes Gaussian, with where N is the number of events, allowing us to define Θ as an effective singleevent H_{0} estimation standard deviation.
In the top panel of Fig. 6, we show the expected number of years of continuous observing required to detect the number of events read on the bottom horizontal axis. We assume an O2like sensitivity^{5} and the detection rates reported in Table 3. From the bottom panel of Fig. 6, the number of Level 2 events required to resolve the H_{0} tension assuming an O2like sensitivity is ∼25. According to the top panel, one would need 100 to 600 years of observation to collect these events, depending on the population model.
Counterparts to GW events detected with different GW sensitivities probe different regions of the distance and inclination parameter space (compare Figs. 3 and 4). Therefore, one cannot deduce the number of events or observing time required to resolve the H_{0} tension for an O4type run from those for an O2type run by simply comparing the GW detection rates. In the next section, we derive the number of events required to resolve the H_{0} tension assuming a longlived O4type run.
In Fig. 7 we show the values of the average effective singleevent standard deviation Θ of the different scenarios, as fit on the curves in Fig. 6. There is a clear boost in the H_{0} convergence speed when considering the information in cos ι from the electromagnetic counterparts. We find that a cos ι precision of 4 deg (Level 3) provides a convergence that is 1.4 times faster than one of 12 deg (Level 2), which is itself about 1.5 times faster when there is no angle information at all from the electromagnetic domain. Roughly, it means that the H_{0} accuracy reached combining 10 Level 2 events is equivalent to that reached by combining 5 Level 3 events.
Fig. 7. Effective singleevent standard deviation Θ in different observing scenarios, assuming all events are at a given electromagnetic information level, in an O2type GW run. 
We also find that detectors with better sensitivities will be able to better constrain H_{0}, even without cos ι measurements from electromagnetic counterpart. This is due to the higher redshifts of the events used to infer H_{0}, as we consider a constant uncertainty in redshift measurements. Considering events without electromagnetic contributions to cos ι for O4type runs, we find an average effective singleevent standard deviation of 14%, consistent with previous simulations (Chen et al. 2018; Gray et al. 2020).
We note that for a given GW sensitivity, there is no difference in this first approach between G16 and WP15 population models as these only impact the probabilities of detecting the electromagnetic counterparts. This is clear in Figs. 6 and 7.
5.3. Considering realistic detection rates
In a real observing run, not all the GW events with a redshift estimation will have cos ι measurements from electromagnetic counterparts. In a second approach, we estimate the H_{0} convergence by including the relative detection rates of the different electromagnetic counterparts. More precisely, we generated threads of events as in the bulk comparison of the previous section, but successively allowing for only Level 1 events, then up to Level 2, then up to Level 3, to quantify the acceleration of H_{0} convergence each EM level allows.
We then simulated the H_{0} measurement for each of them, and combined their measurements throughout the first 100 events. We repeat this process 500 times with varying event time orderings to simulate different realizations of the subset of 100 detections.
Figure 8 shows the H_{0} convergence as a function of the total number of detections, while Fig. 9 shows the singleevent standard deviation Θ for all scenarios. According to Table 2, about 7%−23% of O2 events with associated redshift would have had Level 2 information on ι, and about 1%−13% Level 3. As seen in Fig. 8, this detection fraction is enough to somewhat improve the H_{0} convergence. In other words, the convergence speed allowing for Level 2 and 3 events is larger than with kilonova events alone.
Fig. 8. Same as Fig. 6 (bottom panel), with the bottom axis counting the total number of GW events, regardless of the nature of the EM counterpart. Therefore, this considers the realistic detection rate of events with different electromagnetic counterparts. 
The situation changes drastically when we consider the multimessenger events in O3 and O4type observing runs. In this case, the fraction of Level 2 or Level 3 events are so small that on average they bring no additional improvement. Indeed, the precision on H_{0} obtained with a given number of GW events does not change whether we allow for afterglow counterparts or if we do not. In particular, for an O3like run, only allowing for Level 3 events and assuming the optimistic G16 population prescription could provide a slight acceleration in the H_{0} narrowingdown, while for all the other cases the improvement is negligible. In O4type runs, neither Level 2 or Level 3 events should statistically speed up the convergence of H_{0}.
During an O4type run, Level 2 or 3 events are too rare to significantly improve the H_{0} convergence and shorten the time needed resolve the H_{0} tension. Relying on Level 1 events only is just as fast. In Sect. 6, we argue that, in this case, discarding measurements of the inclination angles from afterglows prevent further biases in the H_{0} estimate.
For O4like runs, Table 3 shows that one Level 1 events are detected every 6 months on average. Figure 8 shows that, to resolve the H_{0} tension problem, 30 such events are required. Thus, with an O4like sensitivity, 15 years of data taking are necessary to collect the number of detections with measured redshift.
6. Discussion
In this article, we studied the prospects of measuring the Hubble constant with GW standard sirens coupled to inclination angle measurements from merger afterglow counterparts. We first studied the potential impact of selection effects in multimessenger cosmology and showed that these were negligible in the H_{0} measurement reported after GW170817. We illustrated how selection effects increase with the GW sensitivity. For events in future GW observing runs, selection effects will be important and should be taken into account.
We then studied the likelihood of detecting the electromagnetic counterparts required for multimessenger measurements of H_{0}. We showed that for O2 and O3type GW sensitivities, their optical magnitudes are not limiting in the detection of kilonova signals, and therefore in acquiring source redshifts. This is true provided GW interferometers are operating and followup facilities available and efficient.
We showed that, in O4like runs and beyond, the GW detection probability will extend further in distance compared to the electromagnetic detection probabilities, largely decreasing the likelihood of detecting electromagnetic counterparts. This will especially be the case when secondgeneration GW interferometers reach their design sensitivity level. Fortunately, the lower number of observable counterparts should be counterbalanced by the advent of large fieldofview, highcadence optical instruments such as the ZTF and LSST, which are not limited in terms of skymap coverage (but only by their limiting magnitude). Also, the advent of additional GW interferometers such as KAGRA (Kagra Collaboration 2019) will largely improve the median GW localization skymap down to ∼40 deg^{2} (Abbott et al. 2020a), paving the way to more effective followup by smaller fieldofview instruments and therefore bettersampled kilonova light curves and better leveraging of events.
In the future, thirdgeneration GW detectors such as the Einstein Telescope (ET, Punturo et al. 2014) and Cosmic Explorer (Reitze et al. 2019) will open up the detection range of even larger redshifts. Naturally, the lowredshift events will still be accessible, with a much better sky resolution that can be covered entirely by small fieldofview facilities. For larger redshifts, followup observations will be limited by the sky localization, with an average resolution of 200 deg^{2} for BNS at z = 0.1 with ET. This coverage issue combined with the unprecedented dimness of the counterparts at these redshifts may call for totally different followup strategies for the highcadence large fieldofview survey facilities. Supposing the source is identified, photometric and spectroscopic followup would still be limited to z = 0.5 and z = 0.3, respectively, for the largest optical telescopes such as the Extremely Large Telescope (Maggiore et al. 2020). Therefore, the access to those events at cosmic redshifts should rely on the observation of counterparts other than kilonova, such as the short gammaray burst, the detection of which could be facilitated by the GW signal being present in the ET band hours before merger.
We finally studied whether the observation of merger afterglow signals and subsequent measurements of cos ι will significantly accelerate the narrowingdown of H_{0} when combined with GW detections in the future. We considered only the afterglow signal–its photometry and imaging–as potential providers of ι measurements. We deemed the other counterparts such as the kilonova and gammaray burst unfit for such a measurement, for their still large modeling uncertainties.
We considered an optimistic measurement model in which all events with an afterglow counterpart contribute a ι measurement with an accuracy comparable to GW1710817, for both afterglow photometry and imaging. This is an optimistic assumption as the uncertainty on ι depends on the number of photometric points detected from the light curve, and thereby on the event distances and density of the followup. Also, only a subset of the followup campaigns is expected to provide such detailed multiwavelength photometric data. Taking the variability of the followup scope and data quality into account is a possible extension to this study.
Furthermore, our analysis assumes that cos ι estimated from EM is accurate, which is a simplifying assumption. For both the afterglow photometry and imaging analysis, the leading uncertainty in the electromagnetic modeling is the treatment of the jet lateral expansion (Ghirlanda et al. 2019). The jet expansion affects both the time of afterglow peak flux (e.g., Duque et al. 2019) and the dynamics of the VLBI image (Fernández et al. 2021), possibly biasing and widening the electromagnetic posterior on ι more than assumed here.
From the point of view of observations, most of the indeterminations in measuring ι from GW170817’s afterglow lie in the very early phases of the light curve and the late phases of the source displacement curve (see the posterior sampling in Ghirlanda et al. 2019). Acquiring early photometric data with deeper searches would have provided a better estimate of ι from this source. Recently, radio points were acquired from GW170817 to flux levels much deeper than the early radio monitoring of GW170817 (Balasubramanian et al. 2021), proving that this is in reach of current radio facilities. We thus advocate for such early deep searches.
Also, we assumed perfect kilonova detection and identification above a magnitude threshold over the accessible sky, and our expression of the jet proper motion in Eq. (6) is clearly an overestimate, leading to overpredict the number of sources with detectable proper displacement.
Finally, the effect of peculiar velocities can bias the estimation of H_{0}. Galaxies’ peculiar motions can be as high as v_{p} = 300 km s^{−1}, with an associated error of their measurement of the same order (Mukherjee et al. 2021a). This corresponds to a redshift correction (and additional uncertainty on H_{0}) of 10% for events such as GW170817, and it is fundamental to take into account. With future GW detectors, which will detect BNSs up to redshift ∼0.1, this type of correction will generally be negligible though it will remain important in the case of closeby, highS/N events.
For all these reasons, our hypotheses are optimistic. Even so, we have found that, for all GW runs after O3, events with afterglow counterparts should prove so rare that, statistically, using the afterglow counterparts when available for the multimessenger measurement of H_{0} on individual events will not bring any acceleration to the measurement of H_{0}.
We found that, for the electromagnetic measurement of cos ι to significantly increase the H_{0} convergence, the GW and electromagnetic detection probabilities should be comparable, or the understanding of electromagnetic emission from compact binary mergers should drastically improve, to the point where a degreelevel precision on ι is accessible from a typical afterglow light curve. Even then, care should be taken with systematic effects in angle measurements, for example, from VLBI imaging, to not create a second tension on H_{0}.
The current state of kilonova modeling does not allow for robust measurements of the inclination angle. This may change in the future as the model uncertainties will probably reduce after more signals are observed. Given the abovementioned rates of kilonova associations with GW, we estimate that kilonovae could accelerate the narrowingdown of H_{0} if a ∼10% calibration can be reached between kilonova data (such as light curves or color evolution) and inclination angle. In this respect, the advent of wide fieldofview, highcadence optical facilities is an asset as they are expected to collect a large sample of kilonovae detected both serendipitously and as GW counterparts. With tens of wellsampled kilonova light curves, the 10% accuracy level for the light curve versus inclination relation may be within reach, especially if these sources are standardizable (Kashyap et al. 2019; Coughlin et al. 2020b,c). This perspective would truly allow multimessenger cosmology to develop.
Our results should not be misunderstood: If the opportunity of making an electromagnetic measurement of cos ι occurs, then it should obviously be made as the improvement on H_{0} from such an event is important. All that we found is that, statistically, being able to make such measurements will not significantly speed up the narrowingdown of H_{0} and the resolution of the Hubble tension on the long run. This is because of the rareness of electromagnetic counterparts and their still toopoor constraints on cos ι.
We proved that the electromagneticprovided cos ι measurements will likely not drive the H_{0} narrowingdown. Therefore, direct biases to H_{0} through electromagnetic mismeasurements of cos ι should not be feared. However, as the detection probabilities of the electromagnetic counterparts should dominate the selection effect for GWEM standard sirens, incorrectly modeled dependence of the kilonova signal on the inclination can lead to H_{0} biases through uncontrolled selection effects, as discussed in Chen et al. (2018), Chen (2020). Correct modeling of the kilonova signal to control the selection effects in followup campaigns should be a point of care for future endeavors in multimessenger cosmology.
Once these selection biases are dealt with, the limiting uncertainty in multimessenger cosmology should be the GW data calibration. This uncertainty is at the level of ∼1% (Karki et al. 2016), below the Hubble constant tension, and therefore should not impede the resolution of the tension by multimessenger cosmology when combining a low number of events. However, a systematic effort on calibration uncertainties when combining a large number of events should be performed.
As the afterglow counterparts should not accelerate the measurement of H_{0}, we can state that the number of multimessenger events necessary to resolve the H_{0} tension is still that given by Chen et al. (2021) and Mortlock et al. (2019), that is, 20−50. This represents about fifteen years of continuous O4level GW observation.
7. Conclusion
The afterglow counterparts of binary neutron star mergers represent viable means to measure the inclination angle of sources, and thereby to improve the standardsiren measurement of the Hubble constant. Afterglows could therefore play the role of narrowing down H_{0} and possibly resolve the Hubble tension with fewer events than by leveraging only the gravitationalwave data and source redshift. To quantify how much faster afterglowenhanced H_{0} measurements could solve the Hubble tension, we carried out a realistic population model considering that every future afterglow counterpart could provide a constraint on the source inclination angle at the same level as GW170817. We found that, while each afterglow allows for a jump in H_{0} precision, events with afterglow counterparts should prove very rare, to the point that allowing for afterglowenhanced measurements should not statistically make any difference in the number of events required. Once models have improved, kilonova light curves could be viable for inclination angle measurements and, as these should be much more frequently acquired, kilonovae could play the leading role in multimessenger cosmology. Whether for kilonova or afterglow counterparts, one must treat selection effects with care so as to not produce yet another tension because upcoming gravitationalwave observing runs will probe distances where selections effects are important.
The afterglow of GW170817 was still detected in the Xray band 1000 days postmerger (Hajela et al. 2020).
Usually a S/N threshold of 8 is assumed for the detection by a single detector Abbott et al. (2020a). As the CF approximant is valid in the high S/N regime, here we assume a threshold of 14. This is equivalent, on average, to a S/N of about 8 in each detector.
We note that the results on afterglow counterpart fractions presented in this section are consistent with the predictions of Duque et al. (2019), though the treatment of the GW and jet motion in the present study is more refined.
We consider an O2like sensitivity for Fig. 6 for computational reasons. Using O3 or O4like sensitivities decreases the number of Level 3 events in our simulation to the point that emulating the H_{0} measurement with even ten such events was computationally infeasible.
Acknowledgments
We thank H. Chen for useful comments during the development of this work. RD acknowledges the Centre National d’Études Spatiales (CNES) for financial support in this research project. SM is supported by the LabEx UnivEarthS (ANR10LABX0023 and ANR18IDEX0001), of the European Gravitational Observatory, and the Paris Center for Cosmological Physics.
References
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017a, ApJ, 848, L12 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017b, Nature, 551, 85 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017c, Phys. Rev. Lett., 119, 161101 [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2019, Phys. Rev. D, 100, 104036 [CrossRef] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020a, Liv. Rev. Rel., 23, 3 [CrossRef] [Google Scholar]
 Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2020b, ApJ, 892, L3 [Google Scholar]
 Andreoni, I., Kool, E. C., Sagues Carracedo, A., et al. 2020, ApJ, 904, 155 [CrossRef] [Google Scholar]
 Ascenzi, S., Oganesyan, G., Branchesi, M., & Ciolfi, R. 2021, J. Plasma Phys., 87, 845870102 [CrossRef] [Google Scholar]
 Balasubramanian, A., Corsi, A., Mooley, K. P., et al. 2021, ApJ, 914, L20 [CrossRef] [Google Scholar]
 Bellm, E. C. 2016, PASP, 128, 084501 [Google Scholar]
 Bellm, E. C., Kulkarni, S. R., Graham, M. J., et al. 2019, PASP, 131, 018002 [NASA ADS] [CrossRef] [Google Scholar]
 Beniamini, P., Petropoulou, M., Barniol Duran, R., & Giannios, D. 2019, MNRAS, 483, 840 [NASA ADS] [CrossRef] [Google Scholar]
 Beniamini, P., Granot, J., & Gill, R. 2020, MNRAS, 493, 3521 [CrossRef] [Google Scholar]
 Biwer, C. M., Capano, C. D., De, S., et al. 2019, PASP, 131, 024503 [Google Scholar]
 ChassandeMottin, E., Leyde, K., Mastrogiovanni, S., & Steer, D. A. 2019, Phys. Rev. D, 100, 083514 [CrossRef] [Google Scholar]
 Chen, H.Y. 2020, Phys. Rev. Lett., 125, 201301 [CrossRef] [Google Scholar]
 Chen, H.Y., Fishbach, M., & Holz, D. E. 2018, Nature, 562, 545 [Google Scholar]
 Chen, H.Y., Vitale, S., & Narayan, R. 2019, Phys. Rev. X, 9, 031028 [Google Scholar]
 Chen, H.Y., Holz, D. E., Miller, J., et al. 2021, Class. Quant. Grav., 38, 055010 [CrossRef] [Google Scholar]
 Coughlin, M. W., Dietrich, T., Margalit, B., & Metzger, B. D. 2019a, MNRAS, 489, L91 [CrossRef] [Google Scholar]
 Coughlin, M. W., Ahumada, T., Anand, S., et al. 2019b, ApJ, 885, L19 [NASA ADS] [CrossRef] [Google Scholar]
 Coughlin, M. W., Dietrich, T., Antier, S., et al. 2020a, MNRAS, 492, 863 [CrossRef] [Google Scholar]
 Coughlin, M. W., Antier, S., Dietrich, T., et al. 2020b, Nat. Commun., 11, 4129 [CrossRef] [Google Scholar]
 Coughlin, M. W., Dietrich, T., Heinzel, J., et al. 2020c, Phys. Rev. Res., 2, 022006 [CrossRef] [Google Scholar]
 Cutler, C., & Flanagan, É. E. 1994, Phys. Rev. D, 49, 2658 [NASA ADS] [CrossRef] [Google Scholar]
 Del Pozzo, W., Li, T. G. F., & Messenger, C. 2017, Phys. Rev. D, 95, 043502 [CrossRef] [Google Scholar]
 Dewdney, P. E., Hall, P. J., Schilizzi, R. T., & Lazio, T. J. L. W. 2009, IEEE Proc., 97, 1482 [Google Scholar]
 Dhawan, S., Bulla, M., Goobar, A., Sagués Carracedo, A., & Setzer, C. N. 2020, ApJ, 888, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Dobie, D., Kaplan, D. L., Hotokezaka, K., et al. 2020, MNRAS, 494, 2449 [CrossRef] [Google Scholar]
 Doctor, Z. 2020, ApJ, 892, L16 [CrossRef] [Google Scholar]
 Duque, R., Daigne, F., & Mochkovitch, R. 2019, A&A, 631, A39 [CrossRef] [EDP Sciences] [Google Scholar]
 Duque, R., Beniamini, P., Daigne, F., & Mochkovitch, R. 2020, A&A, 639, A15 [CrossRef] [EDP Sciences] [Google Scholar]
 Farrow, N., Zhu, X.J., & Thrane, E. 2019, ApJ, 876, 18 [CrossRef] [Google Scholar]
 Feeney, S. M., Mortlock, D. J., & Dalmasso, N. 2018, MNRAS, 476, 3861 [CrossRef] [Google Scholar]
 Fernández, J. J., Kobayashi, S., & Lamb, G. P. 2021, MNRAS, submitted [arXiv:2101.05138] [Google Scholar]
 Fishbach, M., Gray, R., Magaña Hernandez, I., et al. 2019, ApJ, 871, L13 [CrossRef] [Google Scholar]
 Fong, W., Berger, E., Margutti, R., & Zauderer, B. A. 2015, ApJ, 815, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Freedman, W. L. 2017, Nat. Astron., 1, 0169 [NASA ADS] [CrossRef] [Google Scholar]
 Ghirlanda, G., Salafia, O. S., Pescalli, A., et al. 2016, A&A, 594, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ghirlanda, G., Salafia, O. S., Paragi, Z., et al. 2019, Science, 363, 968 [NASA ADS] [CrossRef] [Google Scholar]
 Gill, R., & Granot, J. 2018, MNRAS, 478, 4128 [NASA ADS] [CrossRef] [Google Scholar]
 Gottlieb, O., Nakar, E., & Piran, T. 2019, MNRAS, 488, 2405 [NASA ADS] [CrossRef] [Google Scholar]
 Gray, R., Hernandez, I. M., Qi, H., et al. 2020, Phys. Rev. D, 101, 122001 [CrossRef] [Google Scholar]
 Guidorzi, C., Margutti, R., Brout, D., et al. 2017, ApJ, 851, L36 [CrossRef] [Google Scholar]
 Hajela, A., Margutti, R., Alexander, K. D., et al. 2019, ApJ, 886, L17 [NASA ADS] [CrossRef] [Google Scholar]
 Hajela, A., Margutti, R., Laskar, T., et al. 2020, GRB Coordinates Network, 27414, 1 [Google Scholar]
 Heinzel, J., Coughlin, M. W., Dietrich, T., et al. 2021, MNRAS, 502, 3057 [CrossRef] [Google Scholar]
 Holz, D. E., & Hughes, S. A. 2005, ApJ, 629, 15 [Google Scholar]
 Hotokezaka, K., Kiuchi, K., Shibata, M., Nakar, E., & Piran, T. 2018, ApJ, 867, 95 [NASA ADS] [CrossRef] [Google Scholar]
 Hotokezaka, K., Nakar, E., Gottlieb, O., et al. 2019, Nat. Astron., 3, 940 [NASA ADS] [CrossRef] [Google Scholar]
 Ivezic, Z., Axelrod, T., Brandt, W. N., et al. 2008, Serb. Astron. J., 176, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Kagra Collaboration (Akutsu, T., et al.) 2019, Nat. Astron., 3, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Karki, S., Tuyenbayev, D., Kandhasamy, S., et al. 2016, Rev. Sci. Instrum., 87, 114503 [CrossRef] [Google Scholar]
 Kashyap, R., Raman, G., & Ajith, P. 2019, ApJ, 886, L19 [CrossRef] [Google Scholar]
 Kasliwal, M. M., Nakar, E., Singer, L. P., et al. 2017, Science, 358, 1559 [NASA ADS] [CrossRef] [Google Scholar]
 Kasliwal, M. M., Anand, S., Ahumada, T., et al. 2020, ApJ, 905, 145 [CrossRef] [Google Scholar]
 Kathirgamaraju, A., Giannios, D., & Beniamini, P. 2019, MNRAS, 487, 3914 [CrossRef] [Google Scholar]
 Kawaguchi, K., Shibata, M., & Tanaka, M. 2020, ApJ, 889, 171 [CrossRef] [Google Scholar]
 Lamb, G. P., & Kobayashi, S. 2017, MNRAS, 472, 4953 [NASA ADS] [CrossRef] [Google Scholar]
 Lazzati, D., Perna, R., Morsony, B. J., et al. 2018, Phys. Rev. Lett., 120, 241103 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Maggiore, M., Van Den Broeck, C., Bartolo, N., et al. 2020, JCAP, 2020, 050 [CrossRef] [Google Scholar]
 Mandel, I., Farr, W. M., & Gair, J. R. 2019, MNRAS, 486, 1086 [NASA ADS] [CrossRef] [Google Scholar]
 Mastrogiovanni, S., Steer, D. A., & Barsuglia, M. 2020, Phys. Rev. D, 102, 044009 [CrossRef] [Google Scholar]
 Messenger, C., & Read, J. 2012, Phys. Rev. Lett., 108, 091101 [CrossRef] [Google Scholar]
 Mooley, K. P., Deller, A. T., Gottlieb, O., et al. 2018a, Nature, 561, 355 [Google Scholar]
 Mooley, K. P., Frail, D. A., Dobie, D., et al. 2018b, ApJ, 868, L11 [Google Scholar]
 Mortlock, D. J., Feeney, S. M., Peiris, H. V., Williamson, A. R., & Nissanke, S. M. 2019, Phys. Rev. D, 100, 103523 [Google Scholar]
 Mukherjee, S., Wandelt, B. D., & Silk, J. 2020, MNRAS, 494, 1956 [Google Scholar]
 Mukherjee, S., Lavaux, G., Bouchet, F. R., et al. 2021a, A&A, 646, A65 [CrossRef] [EDP Sciences] [Google Scholar]
 Mukherjee, S., Wandelt, B. D., Nissanke, S. M., & Silvestri, A. 2021b, Phys. Rev. D, 103, 043520 [CrossRef] [Google Scholar]
 Nakar, E., & Piran, T. 2021, ApJ, 909, 114 [CrossRef] [Google Scholar]
 Nakar, E., Piran, T., & Granot, J. 2002, ApJ, 579, 699 [NASA ADS] [CrossRef] [Google Scholar]
 Nakar, E., Gottlieb, O., Piran, T., Kasliwal, M. M., & Hallinan, G. 2018, ApJ, 867, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Nissanke, S., Holz, D. E., Dalal, N., et al. 2013, ArXiv eprints [arXiv:1307.2638] [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Poisson, E., & Will, C. M. 1995, Phys. Rev. D, 52, 848 [Google Scholar]
 Punturo, M., Lück, H., & Beker, M. 2014, in Advanced Interferometers and the Search for Gravitational Waves, ed. M. Bassan, Astrophys. Space Sci. Lib., 404, 333 [CrossRef] [Google Scholar]
 Reitze, D., Adhikari, R. X., Ballmer, S., et al. 2019, BAAS, 51, 35 [NASA ADS] [Google Scholar]
 Resmi, L., Schulze, S., IshwaraChandra, C. H., et al. 2018, ApJ, 867, 57 [NASA ADS] [CrossRef] [Google Scholar]
 Riess, A. G., Casertano, S., Yuan, W., Macri, L. M., & Scolnic, D. 2019, ApJ, 876, 85 [Google Scholar]
 Saleem, M., Pai, A., Misra, K., Resmi, L., & Arun, K. G. 2018, MNRAS, 475, 699 [NASA ADS] [CrossRef] [Google Scholar]
 Schutz, B. F. 1986, Nature, 323, 310 [Google Scholar]
 SoaresSantos, M., Palmese, A., Hartley, W., et al. 2019, ApJ, 876, L7 [CrossRef] [Google Scholar]
 The LIGO Scientific Collaboration, the Virgo Collaboration, Abbott, R., et al. 2020, ArXiv eprints [arXiv:2010.14533] [Google Scholar]
 Troja, E., van Eerten, H., Ryan, G., et al. 2019, MNRAS, 489, 1919 [NASA ADS] [Google Scholar]
 Villar, V. A., Guillochon, J., Berger, E., et al. 2017, ApJ, 851, L21 [NASA ADS] [CrossRef] [Google Scholar]
 Wanderman, D., & Piran, T. 2015, MNRAS, 448, 3026 [Google Scholar]
 Wollaeger, R. T., Korobkin, O., Fontes, C. J., et al. 2018, MNRAS, 478, 3298 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Statistical framework for the inferring of H_{0} with multimessenger data
According to Bayes’ theorem under selection effects, the posterior distribution on H_{0} obtained given a multimessenger set d of GW and electromagnetic data is (Mandel et al. 2019):
Here, p(H_{0}) is prior information before the measurement and p(dH_{0}) is the likelihood of the data. The function β(H_{0}) is generally referred to as selection function or Malmquist bias that accounts for the estimation bias arising from the observation of a distribution in space of astrophysical objects with a fluxlimited survey. The selection function β(H_{0}) corrects for this bias and it reads:
where p_{det}(z, cos ι) is the probability of making the measurement on a system located at redshift z and inclination ι and p_{pop}(z, cos ιH_{0}) is the distribution of the overall population in redshift and inclination.
As explicit in Eq. (A.1), the H_{0} posterior will only be impacted by selection effects if the selection function is not flat, that is, has a significant variation with H_{0}. Our detection probabilities are best expressed in terms of d_{L}:
with L = L_{1}, L_{2}, L_{3} covering our three detection scenario levels.
We consider the formation of binaries uniform in comoving volume (Mortlock et al. 2019; Mastrogiovanni et al. 2020). Therefore, at low redshift, we have .
As explained in Sect. 3, due to the small redshift of the sources we consider, we can neglect the effects of redshift on the detected chirp mass and last stable circular orbit for the GW part, and the electromagnetic wavelengths for counterpart searches. Therefore, the detection probabilities p^{L} we consider have no explicit dependence on the source redshift. In these circumstances, it is clear from Eq. (A.3) that β(H_{0}) will be independent of H_{0} if the cosmology is assumed linear: dz = H_{0}/c dd_{L}. Indeed, in this case, the integrand in β(H_{0}) will be:
leaving no dependence on H_{0}.
Therefore, assuming (i) linear cosmology, (ii) uniformincomovingvolume system formation rate, and (iii) no explicit redshift dependence in the detection process, the selection effects on the multimessenger measurement of H_{0} are null. This fact was already underlined in Mandel et al. (2019).
If, however, a general cosmology is considered, such a simplification does not occur and selection effects can appear. In Fig. A.1, we plot the selection function for our various observing scenarios assuming a flat Universe with current dark matter density Ω_{m} = 0.308 (Planck Collaboration VI 2020).
Fig. A.1. Selection function β(H_{0}) normalized at an arbitrary nominal value of H_{0} = 40 km s^{−1} Mpc^{−1}, for the different observing scenarios we considered. 
All Tables
Description and numerical values of various constants used in the electromagnetic emission and detection models.
Average fraction of GW events observed with different electromagnetic counterpart levels.
Average observation time (years) needed to detect one BNS merger with different EM counterpart levels.
All Figures
Fig. 1. GW detection probability as a function of the BNS luminosity distance and inclination for the detector network HLV using a S/N threshold of 14. The horizontal axis is scaled to the BNS 0.2% response distance , which is Mpc for O2; Mpc for O3; Mpc for O4. 

In the text 
Fig. 2. Probabilities of detection for the electromagnetic counterparts considered in the study: (top), (bottomleft) and (bottomright) under the two different hypotheses for the jet’s kinetic energy distribution denoted G16 and WP15. See Sect. 3.2.2 for details on these hypotheses. 

In the text 
Fig. 3. Histograms for the distributions of inclination angle ι, luminosity distance d_{L} and optimalfilter GW signaltonoise ratio (see Eq. (1)) for all simulated BNS mergers (black line), those detected with GW only (red), and those detected with Level 1 (green line) or Level 3 (blue line) for an O2type run. 

In the text 
Fig. 4. Same as Fig. 3, for an O4type run. 

In the text 
Fig. 5. Posterior on H_{0} obtained for ten BNS events observed with different level of electromagnetic information on ι in an O2type run. Left: Level 1 (no electromagnetic information on ι), center: Level 2 (electromagneticbased ι precision of 12 deg), right: Level 3 (electromagneticbased ι precision of 4 deg). 

In the text 
Fig. 6. Convergence of H_{0} with number of detected events. Bottom: precision on the estimation of H_{0} with 68.3% confidence intervals as a function of the number of events detected at each EM level. Here, all events are assumed to be at a given electromagnetic information level, as denoted in different colors. We also indicate the tension in the Hubble constant and the precision on H_{0} obtained with GW170817 and counterparts. Top: expected number of years of continuous observation required to detect events in various multimessenger scenarios. 

In the text 
Fig. 7. Effective singleevent standard deviation Θ in different observing scenarios, assuming all events are at a given electromagnetic information level, in an O2type GW run. 

In the text 
Fig. 8. Same as Fig. 6 (bottom panel), with the bottom axis counting the total number of GW events, regardless of the nature of the EM counterpart. Therefore, this considers the realistic detection rate of events with different electromagnetic counterparts. 

In the text 
Fig. 9. Same as Fig. 7, but assuming a realistic rate of electromagnetic counterpart detection. 

In the text 
Fig. A.1. Selection function β(H_{0}) normalized at an arbitrary nominal value of H_{0} = 40 km s^{−1} Mpc^{−1}, for the different observing scenarios we considered. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.