tablenum SIXtablenum

# Accurate standard siren cosmology with joint gravitational–wave and
$\gamma$–ray burst observations

Michele Mancarella  0000-0002-0675-508X<EMAIL_ADDRESS>Aix-
Marseille Univ., Université de Toulon, CNRS, CPT, Marseille, France
Dipartimento di Fisica “G. Occhialini”, Universitá degli Studi di Milano-
Bicocca, Piazza della Scienza 3, 20126 Milano, Italy  INFN, Sezione di Milano-
Bicocca, Piazza della Scienza 3, 20126 Milano, Italy  Francesco Iacovelli
0000-0002-4875-5862 Département de Physique Théorique & Gravitational Wave
Science Center (GWSC), Université de Genève, 24 quai Ernest Ansermet, 1211
Genève 4, Switzerland  Stefano Foffa  0000-0002-4530-3051 Département de
Physique Théorique & Gravitational Wave Science Center (GWSC), Université de
Genève, 24 quai Ernest Ansermet, 1211 Genève 4, Switzerland  Niccolò Muttoni
0000-0002-4214-2344 Département de Physique Théorique & Gravitational Wave
Science Center (GWSC), Université de Genève, 24 quai Ernest Ansermet, 1211
Genève 4, Switzerland  Michele Maggiore  0000-0001-7348-047X Département de
Physique Théorique & Gravitational Wave Science Center (GWSC), Université de
Genève, 24 quai Ernest Ansermet, 1211 Genève 4, Switzerland

(August 28, 2024)

###### Abstract

Joint gravitational-wave and $\gamma$–ray bursts (GRB) observations are among
the best prospects for standard siren cosmology. However, the strong selection
effect for the coincident GRB detection, which is possible only for sources
with small inclination angles, induces a systematic uncertainty that is
currently not accounted for. We show that this severe source of bias can be
removed by inferring the a–priori unknown electromagnetic detection
probability directly from multimessenger data. This leads at the same time to
an unbiased measurement of the Hubble constant, to constrain the properties of
GRB emission, and to accurately measure the viewing angle of each source. Our
inference scheme is applicable to real data already in the small–statistics
regime, a scenario that might become reality in the near future. Additionally,
we introduce a novel likelihood approximant for GW events which treats the
dependence on distance and inclination as exact.

_Introduction._ Gravitational waves are rapidly establishing as a new pillar
of concordance cosmology [1]. A Hubble diagram can be constructed with
luminosity distance measurements from gravitational–wave (GW) events in
conjunction with redshift determinations [2]. The detection of a few to tens
“bright sirens”, following GW170817 [3, 4, 5], can lead to an independent
percent–level measurement of the Hubble constant $H_{0}$ [6, 7, 8], required
to set the strong inconsistency between local [9] and CMB [10] measurements.
With such precision, systematics will start playing a prominent role and
prevent solving the Hubble tension unless properly accounted for – in
particular, calibration uncertainty [11, 12], peculiar velocities corrections
[13, 14, 15], and the EM selection effects [16].

The latter is the dominant source of bias when the counterpart is a short
$\gamma$–ray burst (GRB) [16], which is observed only if the angular momentum
of the binary is aligned with the line of sight, i.e. at small inclination.
Since inclination and luminosity distance are inherently correlated in the GW
waveform, with the correlation becoming stronger for small viewing angles,
neglecting this effect leads to a strong bias on $H_{0}$ [16]. Knowing a
priori the GRB emission mechanism, we could impose a prior on the inclination,
but the limited knowledge of these phenomena makes this impossible.
Independent measurements of the inclination from the observation of jet
motions with Very Long Baseline Interferometry (VLBI) can help breaking the
distance–inclination correlation [17, 18, 19, 20], but these are not always
guaranteed [21] and are themselves subject to modeling systematics [16].
Recently, proposals for mitigating the bias in absence of direct information
on the inclination and on the GRB emission were put forward. These rely on
estimating the EM selection effect from a larger sample of BNS without
counterparts [22, 23].

In this work, we instead propose to remove the bias by determining the EM
selection function _solely from the joint sample of GW+GRB observations_. This
scheme can be applied even with just two multimessenger events, in which case
we find that neglecting the EM selection bias can lead to a
$\mathcal{O}(10\%)$ shift of the posterior peak for $H_{0}$. This could happen
already during currently planned runs of the LVK collaboration [24]. As for a
larger number of events, the bias can reach the $\mathchar
21016\relax\,3\sigma$ ($\mathchar 21016\relax\,2\sigma$) level for
$\mathcal{O}(50)$($\mathcal{O}(10)$) detections, making the correction
crucial. Our implementation can be directly applied to real data, demonstrated
with a proof of principle application to GW170817.

_Statistical formulation._ We model the problem with a hierarchical Bayesian
approach [25, 26, 27]. We consider a set of $N_{\rm obs}$ events with
single–event detector–frame parameters
$\\{\boldsymbol{\theta}_{i}\\}_{i=1}^{N_{\rm obs}}$ and data ${\cal
D}=\\{{\cal D}_{\rm GW}^{i},{\cal D}_{\rm EM}^{i}\\}_{i=1}^{N_{\rm obs}}$. The
EM data consist only in a redshift measurement for each event, following the
identification of the host galaxy from the GRB. We instead remain agnostic to
the data behind the detection of the GRB and only use the fact that such
detection happened. For clarity, we introduce a set of boolean variables
$\\{{\rm det}^{i}_{\rm GW},{\rm det}^{i}_{\rm EM}\\}_{i=1}^{N_{\rm obs}}$
where ${\rm\det}^{i}=1$ if the event is detected and 0 otherwise. We assume
that the GW catalog is obtained with a selection cut which is a known,
deterministic function of the observed data only, as customary in hierarchical
Bayesian analyses. On the other hand, the EM detections are also subject to a
strong selection effect, but their selection function is not known a priori,
since we do not know the mechanism behind GRB emission. Moreover, the
possibility of measuring one source parameter – the redshift – depends on
another one – the inclination – which is not directly measured. This
introduces a dependence of the selection function on the (also unknown) single
event parameters $\boldsymbol{\theta}$. More specifically, we will model the
EM detection probability, denoted by $P({\rm det}^{i}_{\rm
EM}|\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\rm EM})$, as a function of
the inclination angle $\iota$ and on the luminosity distance $d_{L}$ (for the
usual scaling of the flux with $d_{L}^{-2}$) with universal hyperparameters
$\boldsymbol{\lambda}_{\rm EM}$ describing the GRB jet structure and the GRB
detection threshold.111We assume that the viewing angle of the GRB is the same
as (or a good proxy of) $\iota$ in the GW signal, and that it has a universal
structure (see e.g. [28, 29]). We do not know the parameters
$\boldsymbol{\lambda}_{\rm EM}$ a priori, but we can infer them from the GW
and redshift data and, crucially, from the information that the counterpart
has been detected, i.e. ${\rm det}^{i}_{\rm EM}=1$ for each event in the
catalog.

The joint posterior on
$\\{\boldsymbol{\lambda},\\{\boldsymbol{\theta}_{i}\\}\\}$ can be written as

$\begin{split}p\big{(}&\boldsymbol{\lambda},\\{\boldsymbol{\theta}_{i}\\}|\\{{\cal
D}_{\rm GW}^{i}\\},\\{{\cal D}_{\rm EM}^{i}\\},\\{{\rm det}_{\rm
GW}^{i}\\},\\{{\rm det}_{\rm EM}^{i}\\}\big{)}\propto\\\
&\frac{\pi(\boldsymbol{\lambda})}{P({\rm det}|\boldsymbol{\lambda})^{N_{\rm
obs}}}\times\prod_{i=1}^{N_{\rm obs}}\mathcal{L}\big{(}{\cal D}_{\rm
GW}^{i}|\boldsymbol{\theta}_{i}\big{)}\,\mathcal{L}\big{(}{\cal D}_{\rm
EM}^{i}|\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\rm c}\big{)}\\\ &\times
p_{\rm pop}(\boldsymbol{\theta}_{i}|\boldsymbol{\lambda})\,P({\rm
det}^{i}_{\rm EM}|\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\rm
EM})\,.\end{split}$ (1)

A derivation and detailed discussion of this posterior is given in App.
A.222See also [30] for a similar discussion in the context of GW data only. In
Eq. (1), $p_{\rm pop}(\boldsymbol{\theta}_{i}|\boldsymbol{\lambda})$ is the
probability of a source to have parameters $\boldsymbol{\theta}_{i}$ given
hyperparameters $\boldsymbol{\lambda}$, $\boldsymbol{\lambda}_{\rm
c}=\\{H_{0},\Omega_{\rm m,0}\\}$ are the cosmological parameters,
$\mathcal{L}\big{(}{\cal D}_{\rm GW}^{i}|\boldsymbol{\theta}_{i}\big{)}$ and
$\mathcal{L}\big{(}{\cal D}_{\rm
EM}^{i}|\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\rm c}\big{)}$ are the
GW and EM likelihoods, $\pi(\boldsymbol{\lambda})$ is the prior on
hyperparameters, and $P(\rm det|\boldsymbol{\lambda})$ denotes the fraction of
detectable events [26]. This term contains both the GW and EM selection
effects, with the latter given by $P({\rm det}^{i}_{\rm
EM}|\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\rm EM})$.

The key difference between Eq. (1) and the standard hierarchical Bayesian
posterior in GW population studies is the presence of the term $P({\rm
det}^{i}_{\rm EM}|\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\rm EM})$ at
the numerator. A similar contribution is not present for GW data because
detectability is a property of data only, which allows us to write
$\mathcal{L}({\rm det}_{\rm GW}^{i},{\cal D}_{\rm
GW}^{i}|\boldsymbol{\theta}_{i})=P({\rm det}_{\rm GW}^{i}|{\cal D}_{\rm
GW}^{i})\mathcal{L}({\cal D}_{\rm
GW}^{i}|\boldsymbol{\theta}_{i})=\mathcal{L}({\cal D}_{\rm
GW}^{i}|\boldsymbol{\theta}_{i})$, where the first equality follows from the
composition of probabilities and the second from the fact that for detected
events $P({\rm det}_{\rm GW}^{i}|{\cal D}_{\rm GW}^{i})=1$. A similar
simplification does not hold for the EM likelihood due to the discussed
dependence of the detection probability on $\boldsymbol{\theta}_{i}$ and
$\boldsymbol{\lambda}_{\rm EM}$. The interpretation of the detection
probability $P({\rm det}^{i}_{\rm
EM}|\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\rm EM})$ in Eq. (1) is that
this term will reshape the GW likelihood essentially acting as a modification
of the prior [31, 30] carrying information about the fact that the detection
of the GRB emission preferentially selects small viewing angles. This breaks
the distance–inclination correlation of each individual GW event likelihood
leading to an unbiased estimate of $H_{0}$ (as well as of other cosmological
parameters) and to the reconstruction of the maximum viewing angle.

_A new GW likelihood approximant._ The GW data are usually given, for each
event $i$, in terms of samples from the posterior probability
$p\big{(}\boldsymbol{\theta}_{i}|{\cal D}_{\rm
GW}^{i}\big{)}\propto\mathcal{L}\big{(}{\cal D}_{\rm
GW}^{i}|\boldsymbol{\theta}_{i}\big{)}\pi({\boldsymbol{\theta}_{i}})$ where
$\pi({\boldsymbol{\theta}_{i}})$ is the prior used to analyze the data.
However, full Bayesian simulations are quite expensive, especially for BNS,
and are often circumvented by resorting to the Fisher Information Matrix (FIM)
approximation, i.e. a Gaussian approximation of the likelihood around its peak
[32]. This however fails precisely in the corner of the parameter space which
is relevant for GRBs, i.e. at small or vanishing inclination angle [33, 34],
regardless of possible regularizations or of the use of priors [32, 35, 33,
34].

Fig. 1: Posterior probability in the distance–inclination plane for the
highest SNR source in the simulations used for the results of this work.
Orange contours are obtained using our new likelihood approximant which
encodes the exact dependence of the likelihood on $d_{L}$ and $\iota$. Blue
contours show for comparison the Fisher matrix approximation, which is
inaccurate for sources at low inclination. Green contours are obtained from a
Bayesian parameter estimation with parallel bilby.

For this study, we develop a new likelihood approximant that treats the
likelihood in the subspace $\\{d_{L},\iota\\}$ as exact, while resorting to
the FIM in the rest of the parameter space. This extends existing approaches,
e.g. [36, 37, 38, 39]. The result is given in terms of a simple combination of
FIM elements, computable with open–source packages [40, 35, 41]. We also
introduce a second extension of the FIM approximation, in order to properly
encode the effects of noise fluctuations [31]. This is based on an expansion
around the point maximizing the likelihood in presence of noise, which we
obtain numerically from an explicit realization of the noise in the frequency
domain. The procedure is described in detail App. B and the result for the
likelihood can be found in Eq. (19).

Fig. 1 shows posterior contours in the distance–inclination plane for the
highest SNR source in our simulation, obtained with our extended likelihood
approximant (orange), compared to the FIM (blue), and, as proof of the
validity of the extended approximant, to a full Bayesian simulation with
parallel bilby [42, 43] (green). We see that the FIM completely mismodels the
shape of the likelihood, while the agreement between our approximant and the
full likelihood is excellent. The biased reconstruction of the marginal
distance posterior in Fig. 1 is a consequence of the small inclination of the
source and is therefore intrinsic to any multimessenger event with a GRB, and
the origin of the systematic uncertainty on $H_{0}$.333The bias is further
exacerbated by the use of a prior on inclination flat in $\cos\iota$,
following the assumptions that inclinations have a uniform distribution on the
sphere.

This simulation scheme is indispensable for studies of third–generation (3G)
detectors, for which an efficient and scalable Bayesian inference sampling
tool is still missing (but see [44, 45, 46, 47] for recent progress), while
detection rates are expected to be large [48, 33].

_Case study._ We now illustrate the method on a simple yet realistic example
for the current generation of GW and EM experiments, namely, a
LIGO–Virgo–KAGRA–LIGO India (LVKI) network [49, 50, 51, 52, 53] and Fermi–GBM
[54].

We simulate a population of BNS with a Madau–Dickinson profile for the merger
rate [55], a Gaussian mass distribution, and inclinations uniformly
distributed on the sphere. We assume an LVKI network at design sensitivity
with a $100\%$ duty cycle and threshold on the observed SNR $\rho_{\rm
obs}\geq 12$ [53]. For the EM detection model, we assume that the GRB jet has
a Gaussian–shaped profile with amplitude $A_{0}$ and half–opening $\iota_{\rm
c}$, compatible with the multimessenger analysis of GW170817 [56] (see Eq.
(25)). We consider an event detected if the flux exceeds the Fermi–GBM–like
sensitivity of $F_{\rm
th}=$2\text{\times}{10}^{-7}\text{\,}\mathrm{erg}\text{\,}{\mathrm{cm}}^{-2}\text{\,}{\mathrm{s}}^{-1}$$
[54]. This model predicts a maximum viewing angle $\Theta_{\rm max}={\rm
min}\\{\iota_{\rm max},\,180^{\circ}-\iota_{\rm max}\\}$ of order $\Theta_{\rm
max}\sim 11^{\circ}$ ($0.2\text{\,}\mathrm{rad}$) in the redshift range of
interest for our sample (see Fig. 4 below). Our simulated dataset consists of
a set of samples from the posterior $p(\boldsymbol{\theta}_{i}|{\cal D}_{\rm
GW}^{i})\propto\mathcal{L}\big{(}{\cal D}_{\rm
GW}^{i}|\boldsymbol{\theta}_{i}\big{)}/\pi(\boldsymbol{\theta}_{i})$ obtained
with the likelihood approximant introduced before, and a redshift measurement
with uncertainty $\sigma_{z}=10^{-3}$ for each GW event. We sample Eq. (1) in
the high–dimensional space of
$\boldsymbol{\lambda},\\{\boldsymbol{\theta}_{i}\\}$ using PyMC [57] and the
Hamiltonian Monte Carlo–based scheme NUTS [58, 59]. The hyperparameters
$\boldsymbol{\lambda}$ include cosmology, mass, and redshift distributions,
together with the parameters $\boldsymbol{\lambda}_{\rm
EM}=\\{A_{0},\,\iota_{\rm c},F_{\rm th}\\}$ describing the EM detection
probability. Tab. 1 in App. C contains a summary, with the same appendix
containing a detailed description of the population model, of the mock catalog
and of the inference scheme.

Fig. 2:  Marginal posterior probability on $H_{0}$, without (red) and with
(green) the inclusion of the EM detection probability. Different line styles
correspond to a different number of events analyzed, as in the legend. The
inset shows the bias on the measurement when neglecting the EM selection
probability.

Fig. 2 shows the effect of neglecting (red lines) and inferring (green lines)
the EM detection probability on the posterior probability on $H_{0}$,
marginalized over all the remaining hyperparameters and individual event
parameters, for 2, 5, 10 and 50 events. When neglecting the EM detection
probability a non–negligible ($\mathchar 21016\relax\,7-10\%$) shift of the
peak is already visible for 2 events, remaining stable as the number of
detections increases. The simultaneous fit of the EM selection function can
remove the bias as far as $\geq 2$ events are included in the analysis. Of
course, for low statistics, the bias is partially compensated by the large
statistical uncertainty. The inset in Fig. 2 quantifies the bias as the ratio
between the shift of the posterior peak from the injected value and the half
$68.3\%$ Highest Density Interval, which we denote by $\sigma$ being this
quantity the standard deviation in the Gaussian case. This reaches the
$\mathchar 21016\relax\,2\sigma$ level between 10 and 20 events and $\mathchar
21016\relax\,3\sigma$ for 50 events. We conclude that while for low detection
rates the shift of the peak could be somewhat compatible with a noise
fluctuation, only the inclusion of the EM selection effect can ensure a robust
measurement. For a larger number of detections, as might be expected in a few
years of O4–O5 observations [60] and should be expected for 3G detectors [48],
the bias is statistically significant and must be corrected.

Fig. 3:  Joint constraint on $H_{0}$ and the half–width of the jet $\iota_{\rm
c}$. Dashed (solid) lines indicate the $95\%$ ($68\%$) C.L..

Fig. 3 shows the correlation between $H_{0}$ and the half–width of the jet
$\iota_{\rm c}$, when the EM detection probability is inferred from the data.
In particular, for $\gtrsim 10$ events, large values of $\iota_{\rm c}$
($\gtrsim 5^{\circ}$) are excluded. This shows that adopting the strategy of
putting a prior on the inclination of single events without inferring it from
data can lead to a biased estimate of $H_{0}$.

Furthermore, our inference scheme allows us to reconstruct the posterior
probabilities of individual source parameters for each event, properly
reweighted by the population prior and the EM detection probability.
Remarkably, an accurate measurements of the viewing angle $\Theta_{i}={\rm
min}\\{\iota_{i},\,180^{\circ}-\iota_{i}\\}$ for each event is obtained just
from the deconvolution of the latter, even in absence of a direct measurement
of this angle in the EM (see [28] for a similar approach). These unbiased
measurements are shown in Fig. 4 as a function of redshift for the case of 10
events as green dots,444We find that 50 events at O5 sensitivity give only a
marginal improvement due to the small horizon of the GW detectors. and recover
the injected values (black stars) within $\mathchar 21016\relax\,1\sigma$ with
good accuracy. The inferred values are below the reconstructed maximum viewing
angle $\Theta_{\rm max}$ at any redshift (green band). In contrast, the
measurement from GW data only (blue points) is largely inaccurate, as expected
(see Fig. 1).

Finally, to provide an idea of the information that we can extract with this
technique from the only event with counterpart detected so far, and to show
that our pipeline can be straightforwardly applied to actual detections, the
violin point in yellow in Fig. 4 shows the result obtained with the analysis
of GW170817 and its counterpart, assuming a Gaussian shaped profile for the
jet.555We use the “high spin” posterior samples available on GWOSC [61], the
redshift given in [62] and compute the detection probability with an injection
set generated with O2 sensitivity and an optimal SNR cut $\rho_{\rm opt}\geq
10$ [63]. We find that the posterior probability on $H_{0}$ is perfectly
compatible with the result of [62] even when including the EM detection
probability, as expected with a sample of 1 event only and in line with the
results of our simulations.

Fig. 4: Reconstruction of the maximum inclination angle of the GRB as a
function of redshift from 10 simulated multimessenger events. The fiducial
model is shown as a black, dashed line. Dots with errorbars indicate the
measurements without (blue) and with (green) the inclusion of the EM detection
probability. Black stars denote the injected values. The yellow violin plot
shows, for comparison, the result for the viewing angle of GW170817 obtained
with the method proposed in this work.

_Conclusions and discussion._ The bias arising from EM selection effects can
be a major limitation for standard siren cosmology already in the short term.
In this work, we provided a new method to properly deconvolve it. We can
summarize the main novelties as follows:

  * •

_Unbiased inference of $H_{0}$ in the short term._ It is possible to
accurately measure $H_{0}$ by inferring the EM detection probability from
multimessenger data only, and without the need of modeling explicitly the GRB
detector, starting from as few as two sources, which could become a concrete
scenario in the near future. In this case, the inference scheme introduced
here would be directly applicable, demonstrated by applying our pipeline to
GW170817.

  * •

_Determination of the GRB emission profile and individual viewing angles._ The
proposed inference scheme leads to a reconstruction of the GRB emission
profile as well as to an accurate and precise measurement of the inclination
of each event. The same can be used in conjunction with a prior on $H_{0}$
(from either CMB or SNe) to determine the GRB structure only from GW data and
redshift measurements (see [64, 65, 29]).

These results open other possibilities:

  * •

_More accurate modeling of the EM emission._ More general models of the EM
selection function can be straightforwardly included, such as a limited
redshift range for the EM detections, a dependence of the EM emission on the
masses, a selection of GW sources based on their sky localization, and the
Field of View of the EM instrument. We plan to address these aspects in future
work.

  * •

_Model–independent approaches._ In this work, for simplicity, the functional
dependence of the flux from inclination used for inference was the same used
to generate the EM data. Importantly, this assumption can be relaxed and is
not essential to the method. One can adopt a more agnostic approach and model
the EM detection probability as a flexible function (e.g. a polynomial or a
Gaussian process) to be less prone to modeling systematics. On the other hand,
using physically motivated models can be of interest for determining if
different emission mechanisms provide a better fit to the data, e.g. by
comparing the respective Bayes factors [64, 65].

  * •

_Independent inclination measurements._ We assumed that independent
inclination measurements (e.g. from VLBI) are not present to stress the
generality of the method. Nonetheless, such measurements can be
straightforwardly included by adding the corresponding likelihood in Eq. (1).
In this case, however, additional sources of bias might be present, as shown
in [16].

  * •

_Population properties and perspectives for 3G detectors._ GW and EM
observatories of next generation are under active development [66, 67, 68].
Besides $H_{0}$, our inference scheme includes the reconstruction of the mass
and redshift distributions of the BNS population, as well as the full
distance–redshift relation, including the matter energy density and the Dark
Energy Equation of State. It can also be straightforwardly extended to test
the effect of modified GW propagation [69, 70]. With 3G sensitivities, we will
be able to precisely determine all those. We will provide forecasts in future
work.

With more multimessenger events in the near future and major new experiments
planned, we hope that this work will contribute to advancing towards accurate
and precise standard siren cosmology both in the short and long term.

_Acknowledgments._ We thank Simone Mastrogiovanni for LVK internal review; Om
Salafia, Matthew Mould, Davide Gerosa, Hsin-Yu Chen for useful comments on the
draft; Marica Branchesi and Samuele Ronchini for useful exchanges in the
initial stage of the work; Massimo Dotti, Costantino Pacilio, Arianna Renzini,
Pippa Cole, and Alberto Colombo for discussions. gwfast is publicly available
at github.com/cosmostatgw/gwfast. PyMC is availabe at pymc.io. M.Manc. is
supported by European Union’s H2020 ERC Starting Grant No. 945155–GWmining and
Cariplo Foundation Grant No. 2021-0555. The work of M.Manc. received support
from the French government under the France 2030 investment plan, as part of
the Initiative d’Excellence d’Aix-Marseille Université – A*MIDEX
AMX-22-CEI-02. F.I, S.F., M.Mag and N.M. are supported by the Swiss National
Science Foundation, grant 200020$\\_$191957, and by the SwissMap National
Center for Competence in Research. This work made use of the clusters
Yggdrasil and Baobab at the University of Geneva. This document has LIGO DCC
number LIGO-P2400167.

## References

  * Moresco  _et al._ [2022] M. Moresco, L. Amati, L. Amendola, _et al._ , Living Rev. Rel. 25, 6 (2022), arXiv:2201.07241 [astro-ph.CO] .
  * Schutz [1986] B. F. Schutz, Nature 323, 310 (1986).
  * Abbott  _et al._ [2017a] B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017a), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2017b] B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ (LIGO Scientific, Virgo, Fermi GBM, INTEGRAL, IceCube, AstroSat Cadmium Zinc Telluride Imager Team, IPN, Insight-Hxmt, ANTARES, Swift, AGILE Team, 1M2H Team, Dark Energy Camera GW-EM, DES, DLT40, GRAWITA, Fermi-LAT, ATCA, ASKAP, Las Cumbres Observatory Group, OzGrav, DWF (Deeper Wider Faster Program), AST3, CAASTRO, VINROUGE, MASTER, J-GEM, GROWTH, JAGWAR, CaltechNRAO, TTU-NRAO, NuSTAR, Pan-STARRS, MAXI Team, TZAC Consortium, KU, Nordic Optical Telescope, ePESSTO, GROND, Texas Tech University, SALT Group, TOROS, BOOTES, MWA, CALET, IKI-GW Follow-up, H.E.S.S., LOFAR, LWA, HAWC, Pierre Auger, ALMA, Euro VLBI Team, Pi of Sky, Chandra Team at McGill University, DFN, ATLAS Telescopes, High Time Resolution Universe Survey, RIMAS, RATIR, SKA South Africa/MeerKAT), Astrophys. J. Lett. 848, L12 (2017b), arXiv:1710.05833 [astro-ph.HE] .
  * Abbott  _et al._ [2019] B. P. Abbott, R. Abbott, T. Abbott, _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 011001 (2019), arXiv:1805.11579 [gr-qc] .
  * Nissanke  _et al._ [2013] S. Nissanke, D. E. Holz, N. Dalal, S. A. Hughes, J. L. Sievers, and C. M. Hirata, (2013), arXiv:1307.2638 [astro-ph.CO] .
  * Chen  _et al._ [2018] H.-Y. Chen, M. Fishbach, and D. E. Holz, Nature 562, 545 (2018), arXiv:1712.06531 [astro-ph.CO] .
  * Feeney  _et al._ [2019] S. M. Feeney, H. V. Peiris, A. R. Williamson, S. M. Nissanke, D. J. Mortlock, J. Alsing, and D. Scolnic, Phys. Rev. Lett. 122, 061105 (2019), arXiv:1802.03404 [astro-ph.CO] .
  * Riess  _et al._ [2022] A. G. Riess, W. Yuan, L. M. Macri, _et al._ , Astrophys. J. Lett. 934, L7 (2022), arXiv:2112.04510 [astro-ph.CO] .
  * Aghanim  _et al._ [2020] N. Aghanim, Y. Akrami, F. Arroja, _et al._ (Planck), Astron. Astrophys. 641, A1 (2020), arXiv:1807.06205 [astro-ph.CO] .
  * Mozzon  _et al._ [2022] S. Mozzon, G. Ashton, L. K. Nuttall, and A. R. Williamson, Phys. Rev. D 106, 043504 (2022), arXiv:2110.11731 [astro-ph.CO] .
  * Huang  _et al._ [2022] Y. Huang, H.-Y. Chen, C.-J. Haster, L. Sun, S. Vitale, and J. Kissel, (2022), arXiv:2204.03614 [gr-qc] .
  * Nicolaou  _et al._ [2020] C. Nicolaou, O. Lahav, P. Lemos, W. Hartley, and J. Braden, Mon. Not. Roy. Astron. Soc. 495, 90 (2020), arXiv:1909.09609 [astro-ph.CO] .
  * Howlett and Davis [2020] C. Howlett and T. M. Davis, Mon. Not. Roy. Astron. Soc. 492, 3803 (2020), arXiv:1909.00587 [astro-ph.CO] .
  * Mukherjee  _et al._ [2021] S. Mukherjee, G. Lavaux, F. R. Bouchet, J. Jasche, B. D. Wandelt, S. M. Nissanke, F. Leclercq, and K. Hotokezaka, Astron. Astrophys. 646, A65 (2021), arXiv:1909.08627 [astro-ph.CO] .
  * Chen [2020] H.-Y. Chen, Phys. Rev. Lett. 125, 201301 (2020), arXiv:2006.02779 [astro-ph.HE] .
  * Mooley  _et al._ [2018] K. Mooley, A. Deller, O. Gottlieb, E. Nakar, G. Hallinan, S. Bourke, D. Frail, A. Horesh, A. Corsi, and K. Hotokezaka, Nature 561, 355 (2018), arXiv:1806.09693 [astro-ph.HE] .
  * Hotokezaka  _et al._ [2019] K. Hotokezaka, E. Nakar, O. Gottlieb, S. Nissanke, K. Masuda, G. Hallinan, K. P. Mooley, and A. Deller, Nature Astron. 3, 940 (2019), arXiv:1806.10596 [astro-ph.CO] .
  * Gianfagna  _et al._ [2024] G. Gianfagna, L. Piro, F. Pannarale, H. Van Eerten, F. Ricci, and G. Ryan, Mon. Not. Roy. Astron. Soc. 528, 2600 (2024), arXiv:2309.17073 [astro-ph.HE] .
  * Gianfagna  _et al._ [2023] G. Gianfagna, L. Piro, F. Pannarale, H. Van Eerten, F. Ricci, G. Ryan, and E. Troja, Mon. Not. Roy. Astron. Soc. 523, 4771 (2023), arXiv:2212.01104 [astro-ph.HE] .
  * Mastrogiovanni  _et al._ [2021] S. Mastrogiovanni, R. Duque, E. Chassande-Mottin, F. Daigne, and R. Mochkovitch, Astron. Astrophys. 652, A1 (2021), arXiv:2012.12836 [astro-ph.HE] .
  * Chen  _et al._ [2023] H.-Y. Chen, C. Talbot, and E. A. Chase, (2023), arXiv:2307.10402 [astro-ph.CO] .
  * Gagnon-Hartman  _et al._ [2023] S. Gagnon-Hartman, J. Ruan, and D. Haggard, Mon. Not. Roy. Astron. Soc. 520, 1 (2023), arXiv:2301.05241 [astro-ph.CO] .
  * Colombo  _et al._ [2022a] A. Colombo, O. S. Salafia, F. Gabrielli, G. Ghirlanda, B. Giacomazzo, A. Perego, and M. Colpi, Astrophys. J. 937, 79 (2022a), arXiv:2204.07592 [astro-ph.HE] .
  * Loredo [2004] T. J. Loredo, AIP Conf. Proc. 735, 195 (2004), arXiv:astro-ph/0409387 .
  * Mandel  _et al._ [2019] I. Mandel, W. M. Farr, and J. R. Gair, Mon. Not. Roy. Astron. Soc. 486, 1086 (2019), arXiv:1809.02063 [physics.data-an] .
  * Vitale  _et al._ [2022] S. Vitale, D. Gerosa, W. M. Farr, and S. R. Taylor, in  _Handbook of Gravitational Wave Astronomy_ (Springer, Singapore, 2022) p. 45, arXiv:2007.05579 [astro-ph.IM] .
  * Farah  _et al._ [2020] A. Farah, R. Essick, Z. Doctor, M. Fishbach, and D. E. Holz, Astrophys. J. 895, 108 (2020), arXiv:1912.04906 [astro-ph.HE] .
  * Salafia  _et al._ [2023] O. S. Salafia, M. E. Ravasio, G. Ghirlanda, and I. Mandel, Astron. Astrophys. 680, A45 (2023), arXiv:2306.15488 [astro-ph.HE] .
  * Mould  _et al._ [2024] M. Mould, C. J. Moore, and D. Gerosa, Phys. Rev. D 109, 063013 (2024), arXiv:2311.12117 [gr-qc] .
  * Essick and Fishbach [2024] R. Essick and M. Fishbach, Astrophys. J. 962, 169 (2024), arXiv:2310.02017 [gr-qc] .
  * Vallisneri [2008] M. Vallisneri, Phys. Rev. D 77, 042001 (2008), arXiv:gr-qc/0703086 .
  * Iacovelli  _et al._ [2022a] F. Iacovelli, M. Mancarella, S. Foffa, and M. Maggiore, Astrophys. J. 941, 208 (2022a), arXiv:2207.02771 [gr-qc] .
  * Dupletsa  _et al._ [2024] U. Dupletsa, J. Harms, K. K. Y. Ng, J. Tissino, F. Santoliquido, and A. Cozzumbo, (2024), arXiv:2404.16103 [gr-qc] .
  * Dupletsa  _et al._ [2023] U. Dupletsa, J. Harms, B. Banerjee, M. Branchesi, B. Goncharov, A. Maselli, A. C. S. Oliveira, S. Ronchini, and J. Tissino, Astron. Comput. 42, 100671 (2023), arXiv:2205.02499 [gr-qc] .
  * Chen  _et al._ [2019] H.-Y. Chen, S. Vitale, and R. Narayan, Phys. Rev. X 9, 031028 (2019), arXiv:1807.05226 [astro-ph.HE] .
  * Chassande-Mottin  _et al._ [2019] E. Chassande-Mottin, K. Leyde, S. Mastrogiovanni, and D. A. Steer, Phys. Rev. D 100, 083514 (2019), arXiv:1906.02670 [astro-ph.CO] .
  * de Souza and Sturani [2023] J. M. S. de Souza and R. Sturani, Phys. Rev. D 108, 043027 (2023), arXiv:2302.07749 [gr-qc] .
  * Li  _et al._ [2024] M. Li, J. Yu, and Z. Pan, (2024), arXiv:2403.01846 [gr-qc] .
  * Borhanian [2021] S. Borhanian, Class. Quant. Grav. 38, 175014 (2021), arXiv:2010.15202 [gr-qc] .
  * Iacovelli  _et al._ [2022b] F. Iacovelli, M. Mancarella, S. Foffa, and M. Maggiore, Astrophys. J. Supp. 263, 2 (2022b), arXiv:2207.06910 [astro-ph.IM] .
  * Ashton  _et al._ [2019] G. Ashton, M. Hübner, P. D. Lasky, C. Talbot, _et al._ , Astrophys. J. Suppl. 241, 27 (2019), arXiv:1811.02042 [astro-ph.IM] .
  * Smith  _et al._ [2020] R. J. E. Smith, G. Ashton, A. Vajpeyi, and C. Talbot, Mon. Not. Roy. Astron. Soc. 498, 4492 (2020), arXiv:1909.11873 [gr-qc] .
  * Nitz and Dal Canton [2021] A. H. Nitz and T. Dal Canton, Astrophys. J. Lett. 917, L27 (2021), arXiv:2106.15259 [astro-ph.HE] .
  * Smith  _et al._ [2021] R. Smith, S. Borhanian, B. Sathyaprakash, _et al._ , Phys. Rev. Lett. 127, 081102 (2021), arXiv:2103.12274 [gr-qc] .
  * Wong  _et al._ [2023] K. W. K. Wong, M. Isi, and T. D. P. Edwards, Astrophys. J. 958, 129 (2023), arXiv:2302.05333 [astro-ph.IM] .
  * Wouters  _et al._ [2024] T. Wouters, P. T. H. Pang, T. Dietrich, and C. Van Den Broeck, (2024), arXiv:2404.11397 [astro-ph.IM] .
  * Ronchini  _et al._ [2022] S. Ronchini, M. Branchesi, G. Oganesyan, B. Banerjee, U. Dupletsa, G. Ghirlanda, J. Harms, M. Mapelli, and F. Santoliquido, Astron. Astrophys. 665, A97 (2022), arXiv:2204.01746 [astro-ph.HE] .
  * Aasi  _et al._ [2015] J. Aasi, B. P. Abbott, R. Abbott, _et al._ (LIGO Scientific), Class. Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc] .
  * Acernese  _et al._ [2015] F. Acernese, M. Agathos, K. Agatsuma, _et al._ (Virgo), Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Aso  _et al._ [2013] Y. Aso, Y. Michimura, K. Somiya, M. Ando, O. Miyakawa, T. Sekiguchi, D. Tatsumi, and H. Yamamoto (KAGRA), Phys. Rev. D 88, 043007 (2013), arXiv:1306.6747 [gr-qc] .
  * LIG [2011] LIGO-India, Proposal of the Consortium for Indian Initiative in Gravitational-wave Observations, https://dcc.ligo.org/LIGO-M1100296/public (2011).
  * Abbott  _et al._ [2020] B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ (KAGRA, LIGO Scientific, Virgo), Living Rev. Rel. 23, 3 (2020), arXiv:1304.0670 [gr-qc] .
  * Meegan  _et al._ [2009] C. Meegan, G. Lichti, P. N. Bhat, _et al._ , Astrophys. J. 702, 791 (2009), arXiv:0908.0450 [astro-ph.IM] .
  * Madau and Dickinson [2014] P. Madau and M. Dickinson, Ann. Rev. Astron. Astrophys. 52, 415 (2014), arXiv:1403.0007 [astro-ph.CO] .
  * Troja  _et al._ [2018] E. Troja, L. Piro, G. Ryan, H. van Eerten, R. Ricci, M. Wieringa, S. Lotti, T. Sakamoto, and S. B. Cenko, Mon. Not. Roy. Astron. Soc. 478, L18 (2018), arXiv:1801.06516 [astro-ph.HE] .
  * Abril-Pla  _et al._ [2023] O. Abril-Pla, V. Andreani, C. Carroll, L. Dong, C. J. Fonnesbeck, M. Kochurov, R. Kumar, J. Lao, C. C. Luhmann, O. A. Martin, M. Osthege, R. Vieira, T. V. Wiecki, and R. Zinkov, PeerJ Computer Science 9, e1516 (2023).
  * Brooks  _et al._ [2011] S. Brooks, A. Gelman, G. Jones, and X.-L. Meng,  _Handbook of Markov Chain Monte Carlo_ (Chapman and Hall/CRC, 2011).
  * Hoffman and Gelman [2014] M. D. Hoffman and A. Gelman, Journal of Machine Learning Research 15, 1593 (2014), arXiv:1111.4246 [stat.CO] .
  * Colombo  _et al._ [2022b] A. Colombo, O. S. Salafia, F. Gabrielli, G. Ghirlanda, B. Giacomazzo, A. Perego, and M. Colpi, Astrophys. J. 937, 79 (2022b), arXiv:2204.07592 [astro-ph.HE] .
  * Abbott  _et al._ [2021] R. Abbott, T. D. Abbott, S. Abraham, _et al._ (LIGO Scientific, Virgo), SoftwareX 13, 100658 (2021), arXiv:1912.11716 [gr-qc] .
  * Abbott  _et al._ [2017c] B. P. Abbott, R. Abbott, T. D. Abbott, _et al._ (LIGO Scientific, Virgo, 1M2H, Dark Energy Camera GW-E, DES, DLT40, Las Cumbres Observatory, VINROUGE, MASTER), Nature 551, 85 (2017c), arXiv:1710.05835 [astro-ph.CO] .
  * Abbott  _et al._ [2023] R. Abbott, T. D. Abbott, F. Acernese, _et al._ (KAGRA, Virgo, LIGO Scientific), Phys. Rev. X 13, 011048 (2023), arXiv:2111.03634 [astro-ph.HE] .
  * Hayes  _et al._ [2020] F. Hayes, I. S. Heng, J. Veitch, and D. Williams, Astrophys. J. 891, 124 (2020), arXiv:1911.04190 [astro-ph.HE] .
  * Hayes  _et al._ [2023] F. Hayes, I. S. Heng, G. Lamb, E.-T. Lin, J. Veitch, and M. J. Williams, Astrophys. J. 954, 92 (2023), arXiv:2305.06275 [astro-ph.HE] .
  * Maggiore  _et al._ [2020] M. Maggiore, C. Van Den Broeck, N. Bartolo, _et al._ , JCAP 2020 (03), 050, arXiv:1912.02622 [astro-ph.CO] .
  * Evans  _et al._ [2021] M. Evans, R. X. Adhikari, C. Afle, _et al._ , (2021), arXiv:2109.09882 [astro-ph.IM] .
  * Branchesi  _et al._ [2023] M. Branchesi, M. Maggiore, D. Alonso, _et al._ , JCAP 2023 (07), 068, arXiv:2303.15923 [gr-qc] .
  * Belgacem  _et al._ [2018a] E. Belgacem, Y. Dirian, S. Foffa, and M. Maggiore, Phys. Rev. D 97, 104066 (2018a), arXiv:1712.08108 [astro-ph.CO] .
  * Belgacem  _et al._ [2018b] E. Belgacem, Y. Dirian, S. Foffa, and M. Maggiore, Phys. Rev. D 98, 023510 (2018b), arXiv:1805.08731 [gr-qc] .
  * Finke  _et al._ [2021] A. Finke, S. Foffa, F. Iacovelli, M. Maggiore, and M. Mancarella, JCAP 2021 (08), 026, arXiv:2101.12660 [astro-ph.CO] .
  * Maggiore [2007] M. Maggiore,  _Gravitational Waves. Vol. 1: Theory and Experiments_, Oxford Master Series in Physics (Oxford University Press, 2007).
  * Kole  _et al._ [2023] M. Kole, F. Iacovelli, M. Mancarella, and N. Produit, Astron. Astrophys. 669, A77 (2023), arXiv:2211.12403 [astro-ph.HE] .
  * Farrow  _et al._ [2019] N. Farrow, X.-J. Zhu, and E. Thrane, Astrophys. J. 876, 18 (2019), arXiv:1902.03300 [astro-ph.HE] .
  * Yu  _et al._ [2021] J. Yu, H. Song, S. Ai, H. Gao, F. Wang, Y. Wang, Y. Lu, W. Fang, and W. Zhao, Astrophys. J. 916, 54 (2021), arXiv:2104.12374 [astro-ph.HE] .
  * Zhang  _et al._ [2004] B. Zhang, X.-Y. Dai, N. M. Lloyd-Ronning, and P. Meszaros, Astrophys. J. Lett. 601, L119 (2004), arXiv:astro-ph/0311190 .
  * Howell  _et al._ [2019] E. J. Howell, K. Ackley, A. Rowlinson, and D. Coward, Monthly Notices of the Royal Astronomical Society 485, 1435 (2019), arXiv:1811.09168 [astro-ph.HE] .
  * Belgacem  _et al._ [2019] E. Belgacem, Y. Dirian, S. Foffa, E. J. Howell, M. Maggiore, and T. Regimbau, JCAP 2019 (08), 015, arXiv:1907.01487 [astro-ph.CO] .
  * Karamanis  _et al._ [2021] M. Karamanis, F. Beutler, and J. A. Peacock, Mon. Not. Roy. Astron. Soc. 508, 3589 (2021), arXiv:2105.03468 [astro-ph.IM] .
  * Karamanis and Beutler [2021] M. Karamanis and F. Beutler, Stat. Comput. 31, 61 (2021), arXiv:2002.06212 [stat.ML] .
  * Tiwari [2018] V. Tiwari, Class. Quant. Grav. 35, 145009 (2018), arXiv:1712.00482 [astro-ph.HE] .
  * Farr  _et al._ [2019] W. M. Farr, M. Fishbach, J. Ye, and D. Holz, Astrophys. J. Lett. 883, L42 (2019), arXiv:1908.09084 [astro-ph.CO] .
  * Pedregosa  _et al._ [2011] F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Prettenhofer, R. Weiss, V. Dubourg, _et al._ , J. Mach. Learn. Res. 12, 2825 (2011), arXiv:1201.0490 [cs.LG] .

##  Appendix A Derivation of the hierarchical posterior

We consider a set of events with single–event parameters
$\\{\boldsymbol{\theta}_{i}\\}_{i=1}^{N_{\rm obs}}$ and data $\\{{\cal
D}\\}_{i=1}^{N_{\rm obs}}$. We also introduce a set of boolean variables
$\\{{\rm det}^{i}\\}_{i=1}^{N_{\rm obs}}$ where ${\rm\det}^{i}=1$ if the event
is detected and 0 otherwise. As common in GW population studies, the
hierarchical model is an inhomogeneous Poisson process [25, 26]. The expected
number of sources is defined as $N_{\rm exp}(\boldsymbol{\lambda})=N\,P(\rm
det|\boldsymbol{\lambda})$ with

$P({\rm
det}|\boldsymbol{\lambda})=\int\differential\boldsymbol{\theta}\differential{\cal
D}\,p_{\rm pop}(\boldsymbol{\theta}|\boldsymbol{\lambda})\mathcal{L}({\cal
D}|\boldsymbol{\theta})\,P({\rm det}|{\cal
D},\boldsymbol{\theta},\boldsymbol{\lambda})\,.$ (2)

In the literature, $P(\rm det|\boldsymbol{\lambda})$ is often denoted as
$\alpha({\boldsymbol{\lambda}})$,666The notation
$\beta({\boldsymbol{\lambda}})$ or $\sigma({\boldsymbol{\lambda}})$ is also
adopted. and represents the fraction of detectable events given
hyperparameters ${\boldsymbol{\lambda}}$, while usually the integral
$\int\differential{\cal D}\mathcal{L}({\cal D}|\boldsymbol{\theta})\,P({\rm
det}|{\cal
D},\boldsymbol{\theta},\boldsymbol{\lambda})=\int_{\raisebox{2.0pt}{$\scriptstyle
f({\cal D})\geq 0$}}\differential{\cal D}\mathcal{L}({\cal
D}|\boldsymbol{\theta})$ is denoted as $P_{\rm det}(\boldsymbol{\theta})$ [or
more correctly $P({\rm det}|\boldsymbol{\theta})$], under the assumption that
$P({\rm det}_{i}|{\cal D},\boldsymbol{\theta},\boldsymbol{\lambda})=P({\rm
det}|{\cal D})$ and that $P({\rm det}|{\cal D})$ is an indicator function
equal to 1 if $f({\cal D})\geq 0$ and 0 otherwise, see below. Concretely, in
GW population studies one often assumes $f({\cal D})=\rho_{\rm obs}-\rho_{\rm
th}$ where $\rho_{\rm obs}$ is the observed SNR in the detector and $\rho_{\rm
th}$ a given threshold that defines a detection.

Including explicitly the variables $\\{{\rm det}^{i}\\}_{i=1}^{N_{\rm obs}}$,
the hierarchical likelihood can be written as [26, 31]

$\begin{split}p\big{(}\\{\boldsymbol{\theta}_{i}\\},\\{{\cal D}_{i}\\},\\{{\rm
det}_{i}\\}|\boldsymbol{\lambda}\big{)}=&\,N^{N_{\rm obs}}e^{-NP(\rm
det|\boldsymbol{\lambda})}\\\ &\,\times\prod_{i=1}^{N_{\rm obs}}p({\rm
det}_{i},{\cal
D}_{i},\boldsymbol{\theta}_{i}|\boldsymbol{\lambda})\,,\end{split}$ (3)

where in full generality

$\begin{split}p({\rm det}_{i},{\cal
D}_{i},\boldsymbol{\theta}_{i}|\boldsymbol{\lambda})&=\mathcal{L}({\rm
det}_{i},{\cal D}_{i}|\boldsymbol{\theta}_{i})p_{\rm
pop}(\boldsymbol{\theta}_{i}|\boldsymbol{\lambda})\\\ &=P({\rm det}_{i}|{\cal
D}_{i},\boldsymbol{\theta}_{i},\boldsymbol{\lambda})\mathcal{L}({\cal
D}_{i}|\boldsymbol{\theta}_{i})p_{\rm
pop}(\boldsymbol{\theta}_{i}|\boldsymbol{\lambda})\,.\end{split}$ (4)

At this point, the usual assumption in GW astronomy is that detectability is a
function of the data only, and not of the (true, unknown) source parameters
$\boldsymbol{\theta}_{i}$, nor of the hyperparameters $\boldsymbol{\lambda}$.
This implies $P({\rm det}_{i}|{\cal
D}_{i},\boldsymbol{\theta}_{i},\boldsymbol{\lambda})=P({\rm det}_{i}|{\cal
D}_{i})$. Finally, assuming that detection is a deterministic function of the
data, one has $P({\rm det}_{i}|{\cal D}_{i})=1\,\forall i$ since each event in
the catalog has obviously been detected.777It can be shown that the likelihood
is independent of this term even in the more general case in which the
detection is a probabilistic function of the data provided that it is
independent of the true source parameters $\boldsymbol{\theta}_{i}$ [31]. This
is what we assume for GW data also in this work – an event is detected if the
_observed_ SNR exceeds a given threshold.

However, we might be in a situation where detectability depends on some
(unknown) source parameters that are not measured by the experiment. In this
case, we have no data available for all the parameters determining
detectability, but we know that a detection happened, and we can include this
information. In particular, we assume that detectability has an a–priori
unknown dependence on some of the (equally unknown) source parameters
$\boldsymbol{\theta}_{i}$ through a subset of hyperparameters
$\tilde{\boldsymbol{\lambda}}$ that we wish to determine from data. In this
case, $P({\rm det}_{i}|{\cal
D}_{i},\boldsymbol{\theta}_{i},\tilde{\boldsymbol{\lambda}})$ is not equal to
one for all $\boldsymbol{\theta}_{i}$. This term will appear in the likelihood
and the parameters $\tilde{\boldsymbol{\lambda}}$ can be determined from the
data. This is the case for the EM data in this work, as we now detail.

Turning to the specific case of interest for this paper, we assume that the
likelihood is composed by GW ad EM data with independent noise,
$\mathcal{L}({\cal D}|\boldsymbol{\theta},z)=\mathcal{L}\big{(}{\cal D}_{\rm
GW}|\boldsymbol{\theta}\big{)}\,\mathcal{L}\big{(}z_{\rm EM}^{\rm
obs}|z\big{)}$. We denote by $\\{{\cal D}_{\rm GW}^{i},{\cal D}_{\rm
EM}^{i}\\}_{i=1}^{N_{\rm obs}}$ the two datasets, with ${\cal D}_{\rm
EM}=\\{z_{\rm EM}^{\rm obs}\\}$. The parameters $\boldsymbol{\theta}$ refer to
detector–frame parameters of the GW waveform while $z$ denotes the true
redshift. The GW likelihood is discussed in App. B. As anticipated, the GW
measurement is modeled as a deterministic function of the data.

Turning to the EM sector, crucially, in this work we assume that we are not
able to obtain a measurement of all the variables determining the EM
detections. Explicitly, for the problem under consideration in this paper, the
EM data consist only in a redshift measurement, following from the
identification of the host galaxy from the GRB emission.888Notice that, even
if we assume an EM likelihood for the redshift only, crucially the counterpart
detection provides also a very accurate measurement of the sky position of the
source and of the GPS time of the event. This information in our setting is
already incorporated in the GW likelihood, which we compute keeping these
parameters fixed. We expect that the detection of the source can happen only
under certain conditions on its inclination (namely, the inclination has to be
small enough in order to observe a relativistic jet), but we remain agnostic
about the data from the GRB observatory (see [29] for an explicit modeling),
and assume that we are just given an observed redshift $z_{\rm EM}^{\rm obs}$
(for which we will assume a Gaussian likelihood with standard deviation
$\sigma_{z}$).

We then assume that the EM detection has a dependence on the source parameters
and the hyperparameters $\boldsymbol{\lambda}_{\rm EM}$ describing the GRB jet
structure. Hence, $P({\rm det}_{i}|{\cal
D}_{i},\boldsymbol{\theta}_{i},\tilde{\boldsymbol{\lambda}})=P({\rm det}_{\rm
GW}|{\cal D}_{\rm GW})P({\rm det}_{\rm
EM}|\boldsymbol{\theta},\boldsymbol{\lambda}_{\rm EM})=P({\rm det}_{\rm
EM}|\boldsymbol{\theta},\boldsymbol{\lambda}_{\rm EM})$.

Finally, the population prior is

$p_{\rm
pop}(\boldsymbol{\theta}|\boldsymbol{\lambda})\propto{\frac{\psi(z,\boldsymbol{\lambda}_{\rm
z})}{1+z}\frac{\differential V_{c}}{\differential
z}\;\frac{p(m_{1},m_{2}|\boldsymbol{\lambda}_{\rm
m})\,}{\frac{\partial{d_{L}}}{\partial z}(z,\boldsymbol{\lambda}_{\rm
c})\,(1+z)^{2}}}\Bigg{|}_{\raisebox{6.0pt}{$\scriptstyle\begin{subarray}{c}m_{i}=\frac{m_{i}^{\rm
det}}{1+z(d_{L},\boldsymbol{\lambda}_{\rm c})}\\\
z=z(d_{L},\boldsymbol{\lambda}_{\rm c})\end{subarray}$}}\,,$ (5)

where $p(m_{1},m_{2}|\boldsymbol{\lambda}_{\rm m})$ is the source–frame mass
distribution with hyperparameters $\boldsymbol{\lambda}_{\rm m}$,
$\psi(z,\boldsymbol{\lambda}_{\rm z})$ is the source–frame merger rate
distribution with hyperparameters $\boldsymbol{\lambda}_{\rm z}$,
$\differential V_{c}/\differential z$ is the differential comoving volume, the
factor $1/(1+z)$ takes into account time dilation from source to detector
frame, and the factor $(1+z)^{2}\partial{d_{L}}/\partial
z(z,\boldsymbol{\lambda}_{\rm c})$ at the denominator comes from conversion of
the prior from detector– to source–frame variables. Explicitly, the
distance–redshift relation in a flat $\Lambda$CDM Universe is given by

$d_{L}(z,\boldsymbol{\lambda}_{\rm
c})=\frac{c(1+z)}{H_{0}}\int_{0}^{z}\frac{\differential
z^{\prime}}{\sqrt{\Omega_{\rm m,0}(1+z^{\prime})^{3}+1-\Omega_{\rm m,0}}}\,.$
(6)

The cosmological hyperparameters are denoted by $\boldsymbol{\lambda}_{\rm
c}=\\{H_{0},\Omega_{\rm m,0}\\}$. The full set of hyperparameters is defined
by $\boldsymbol{\lambda}=\\{\boldsymbol{\lambda}_{\rm
c},\boldsymbol{\lambda}_{\rm m},\boldsymbol{\lambda}_{\rm
z},\boldsymbol{\lambda}_{\rm EM}\\}$.

We will work in detector–frame; hence, in Eq. (5) the source–frame quantities
are computed from the detector–frame ones with a dependence on the cosmology.
In particular, the redshift can always be re–expressed as a function of the
luminosity distance and of the cosmological parameters
$\boldsymbol{\lambda}_{\rm c}$ [71]; we will assume to have performed this
change of variables and denote the EM likelihood as $\mathcal{L}\big{(}{\cal
D}_{\rm EM}^{i}|\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\rm
c}\big{)}=\mathcal{L}\big{(}z_{{\rm EM},i}^{\rm
obs}|z(d_{L,i},\boldsymbol{\lambda}_{\rm c})\big{)}$.

Putting everything together in Eq. (3), adding a prior
$\pi(\boldsymbol{\lambda})$ on the hyperparameters, and marginalizing over the
overall number of expected events $N$ with a prior $\propto 1/N$ [26], we
arrive at the hierarchical likelihood in Eq. (1).

Interestingly, we can also make an explicit link to the modeling of the
likelihood of the GRB experiment [29].999We are grateful to Om Salafia for
discussions on this point. In this case, one would observe a flux with the
observation described by a likelihood of the form $\mathcal{L}(F^{\rm
obs}|F)$. Each term in the product at the numerator of Eq. (1) would then be
of the form

$\begin{split}&\mathcal{L}\big{(}{\cal D}_{\rm
GW}^{i}|\boldsymbol{\theta}_{i}\big{)}\,\mathcal{L}\big{(}z_{{\rm EM},i}^{\rm
obs}|z(d_{L,i},\boldsymbol{\lambda}_{\rm c})\big{)}\,\mathcal{L}(F^{\rm
obs}_{i}|F_{i})\,\\\ &\ \times p_{\rm
pop}(\boldsymbol{\theta}_{i}|\boldsymbol{\lambda})\,p_{\rm
pop}(F|\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\rm jet})\,,\end{split}$
(7)

where $\boldsymbol{\lambda}_{\rm jet}$ are the parameters describing the GRB
jet structure (e.g. $\\{A_{0},\iota_{\rm c}\\}$ in this paper). Being agnostic
on the GRB data is equivalent to marginalize the likelihood over the observed
flux $F^{\rm obs}_{i}$, which gives for each event a contribution

$\int\differential F^{\rm obs}_{i}\,\mathcal{L}(F^{\rm obs}_{i}|F_{i})\,p_{\rm
pop}(F_{i}|\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\rm jet})\,.$ (8)

The population prior on the flux is given by the relation between the latter
and the source parameters determining its shape, $p_{\rm
pop}(F|\boldsymbol{\theta}_{i},\boldsymbol{\lambda})=\delta\big{(}F-F(\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\rm
jet})\big{)}$, so $F$ can be integrated out of the likelihood (see e.g. Eq.
(25) below for the model of
$F(\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\rm jet})$ used in this
work). Moreover, the integral in Eq. (8) can be restricted to the domain of
observable EM data (i.e. those with $F^{\rm obs}>F_{\rm th}$ for a given
detection threshold $F_{\rm th}$) by writing $\mathcal{L}(F^{\rm
obs}_{i}|F_{i})=\mathcal{L}({\rm det}_{\rm EM,i},F^{\rm
obs}_{i}|F_{i})+\mathcal{L}(\neg{\rm det}_{\rm EM,i},F^{\rm obs}_{i}|F_{i})$
where $\neg{\rm det}_{\rm EM,i}$ denotes a non–detection. By definition, the
second term is zero while the first is $\mathcal{L}({\rm det}_{\rm
EM,i},F^{\rm obs}_{i}|F_{i})=P({\rm det}_{\rm EM,i}|F^{\rm
obs}_{i})\mathcal{L}(F^{\rm obs}_{i}|F_{i})=\mathcal{L}(F^{\rm
obs}_{i}|F_{i})$ because for detected events $P({\rm det}_{\rm EM,i}|F^{\rm
obs}_{i})=1$. So, Eq. (8) becomes

$\int_{F^{\rm obs}_{i}\geq F_{\rm th}}\differential F^{\rm
obs}_{i}\,\mathcal{L}(F^{\rm
obs}_{i}|F(\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\rm jet}))\equiv
P({\rm det}^{i}_{\rm EM}|\boldsymbol{\theta}_{i},\boldsymbol{\lambda}_{\rm
EM})\,,$ (9)

which leads again to Eq. (1), and we defined $\boldsymbol{\lambda}_{\rm
EM}=\\{A_{0},\iota_{\rm c},F_{\rm th}\\}$.

##  Appendix B Extended likelihood approximant

###  B.1 Effects of detector noise

We start by discussing the inclusion of the effects of detector noise in the
simulation. The GW likelihood is

$-2\,\log\mathcal{L}({\cal D}_{\rm GW}\;|\;{\bf\it\theta})\propto\left({\cal
D}_{\rm GW}-h({\bf\it\theta})\,|\,{\cal D}_{\rm
GW}-h({\bf\it\theta})\right)\,,$ (10)

with ${\cal D}_{\rm GW}=h_{0}+n$ denoting the signal,
$h_{0}=h({\bf\it\theta}_{0})$ the waveform evaluated at the true parameters
${\bf\it\theta}_{0}$, and $n$ the noise. The notation $(\cdot|\cdot)$
indicates the inner product in the frequency domain [72],
$(a|b)=4{\text{Re}}\int_{0}^{\infty}\tilde{a}^{*}(f)\tilde{b}(f)/S_{n}(f)$,
with $S_{n}(f)$ the noise one–sided power spectral density (PSD), defined by
$\langle\tilde{n}^{*}(f)\tilde{n}(f^{\prime})\rangle=\delta(f-f^{\prime})S_{n}(f)/2$.
The usual FIM approximation consists in expanding the likelihood at linear
order around the point $h_{0}$, i.e. $h({\bf\it\theta})\approx
h_{0}+h_{i}\delta{{\theta}}^{i}$ (with
$h_{i}\equiv{\partial_{{\theta}_{i}}h}|_{{\bf\it\theta}_{0}}$ and
$\delta\theta^{i}=\theta^{i}-\theta_{0}^{i}$), obtaining 101010Here and in the
following, we drop parameter–independent normalization factors.

$-2\log\mathcal{L}({\cal D}_{\rm
GW}\,|\,{\bf\it\theta})\propto\delta{{\theta}}^{i}\delta{{\theta}}^{j}\left(h_{i}\,|\,h_{j}\right)-2\delta{{\theta}}^{i}\left(n|\,h_{i}\right)\,.$
(11)

The covariance ${(h_{i}\,|\,h_{j})}^{-1}$ obtained this way does not encode
the effect of noise fluctuations, as it is computed at the point $h_{0}$. This
is not fully consistent with the selection cut applied to obtain the GW
catalog, which, as discussed previously, should depend only on the observed
data, including noise fluctuations [31].111111In practice the impact of
ignoring this effect on the results of simulations depends on the specific
problem, sensitivity, and catalog size, and should be evaluated case–by–case
[31].

We include the effect of detector noise as follows. We restrict for clarity to
the case of the FIM, neglecting higher–order corrections, but the argument
does not depend on this assumption. For each simulated event with true
parameters ${\bf\it\theta}_{0}$, and for each GW detector in the network under
consideration, we generate an explicit realization of the noise from the PSD
in the frequency domain.

For events that pass the selection cut on the observed SNR $\rho_{\rm obs}$,
computed on the specific noise realization, we then find the point
${\bf\it\hat{\theta}}$ that minimizes the negative log–likelihood in Eq. (10)
via numerical minimization.121212More specifically, after drawing a specific
noise realization over a linear frequency grid, we perform a minimization
helped by the fact that we know ${\bf\it\theta}_{0}$, so the initialization
helps the minimizer to converge quickly, and we incorporate physical prior
ranges for all parameters. Then, we consider the expansion of the likelihood
around the ML point ${\bf\it\hat{\theta}}$, i.e. $h({\bf\it\theta})\approx
h({\bf\it\hat{\theta}})+\hat{h}_{i}\delta\hat{\theta}^{i}$ with
$\delta\hat{\theta}^{i}=\theta^{i}-\hat{\theta}^{i}$. The notation
$\hat{h}_{i}$ indicates the partial derivative evaluated at the point
${\bf\it\hat{\theta}}$:
$\hat{h}_{i}\equiv{\partial_{{\theta}_{i}}h}|_{{\bf\it\hat{\theta}}}$. We also
denote $h({\bf\it\hat{\theta}})=\hat{h}$.

At first order we get

$\begin{split}-2\log\mathcal{L}({\cal D}_{\rm
GW}\,|\,{\bf\it\theta})\propto&-\delta\hat{{\theta}}^{i}\delta\hat{{\theta}}^{j}\left(\hat{h}_{i}\,|\,\hat{h}_{j}\right)\\\
&+2\delta\hat{{\theta}}^{i}\left(h_{0}+n-\hat{h}\,|\,\hat{h}_{i}\right)\,.\end{split}$
(12)

The first term is the usual expansion with a Fisher matrix computed around the
ML point, while the second differs from Eq. (11) because we are expanding the
likelihood around ${\bf\it\hat{\theta}}$ and not around
${{\bf\it\theta}}_{0}$. By definition, the point ${\bf\it\hat{\theta}}$ is
such that (neglecting overall normalization factors that are not relevant)
$-2\log\mathcal{L}({\cal D}_{\rm GW}\,|\,{\bf\it\hat{\theta}})\approx 0$, so
we must have ${\cal D}_{\rm GW}=h_{0}+n\approx
h({\bf\it\hat{\theta}})=\hat{h}$, i.e the second term vanishes. In conclusion,
after finding the ML point, in the FIM approximation we can just approximate
the likelihood as a multivariate Gaussian centered around
${\bf\it\hat{\theta}}$ and with covariance
$(\hat{h}_{i}\,|\,\hat{h}_{j})^{-1}$.

###  B.2 Explicit form of the likelihood approximant

Here we present some expressions of the extended likelihood approximant, which
is exact in some subsets of parameters $\boldsymbol{\beta}$; while
$\boldsymbol{\beta}=\\{d_{L},\iota\\}$ is the focus of the present work, it is
also worth considering other cases.

To begin with, we consider the case $\boldsymbol{\beta}=\\{d_{L}\\}$. The
waveform parameters are then split as
$\boldsymbol{\theta}\equiv\\{d_{L},\bar{\boldsymbol{\theta}}\\}$ and the
waveform can be expanded around the ML point as

$h(d_{L}\,,\bar{\boldsymbol{\theta}})\simeq
h(d_{L}\,,\hat{\bar{\boldsymbol{\theta}}})+(\bar{\theta}^{i}-\hat{\bar{\theta}}^{i})\cdot\partial_{\bar{\theta}_{i}}h(d_{L}\,,\bar{\boldsymbol{\theta}})\big{|}_{\bar{\boldsymbol{\theta}}=\hat{\bar{\boldsymbol{\theta}}}}\,.$
(13)

The approximant obtained with this expansion will encode the effect of noise
fluctuations in its width. Replacing this expansion into Eq. (10) leads to

$\begin{split}-2\,\log{\cal L}({\cal D}_{\rm
GW}|d_{L}\,,&\bar{\boldsymbol{\theta}})\propto\left(\frac{\hat{d}_{L}}{d_{L}}\right)^{2}\Bigg{[}\left(\delta
d_{L}\right)^{2}\Gamma_{d_{L}d_{L}}\\\ &\ +2\delta\bar{\theta}_{i}\ \delta
d_{L}\Gamma_{id_{L}}+\delta\bar{\theta}_{i}\Gamma_{ij}\delta\bar{\theta}_{j}\Bigg{]}\,,\end{split}$
(14)

with

$\delta d_{L}\equiv
d_{L}-\hat{d}_{L}\,,\quad\delta\bar{\theta}_{i}=\bar{\theta}^{i}-\hat{\bar{\theta}}^{i}\,.$
(15)

The prefactor $(\hat{d}_{L}/d_{L})^{2}$ naturally suppresses the likelihood
for small values of $d_{L}$, thus making it unnecessary to enforce the prior
$d_{L}>0$. This version of the likelihood could be useful to treat nearby GW
sources, where the Taylor expansion in $d_{L}$ is less likely to be accurate.

We next consider the case $\boldsymbol{\beta}=\\{d_{L}\,,\iota\\}$, which has
been used to produce the results of this paper. We split the waveform
parameters as
$\boldsymbol{\theta}\equiv\\{d_{L},\iota,\bar{\boldsymbol{\theta}}\\}$, and
linearly expand the likelihood in Eq. (10) only in the subspace of the
parameters $\bar{\boldsymbol{\theta}}$. This time the waveform expansion can
be written as

$h(d_{L}\,,\iota\,,\bar{\boldsymbol{\theta}})\simeq
h(d_{L}\,,\iota\,,\hat{\bar{\boldsymbol{\theta}}})+(\bar{\theta}^{i}-\hat{\bar{\theta}}^{i})\cdot\partial_{\bar{\theta}_{i}}h(d_{L}\,,\iota\,,\bar{\boldsymbol{\theta}})\big{|}_{\bar{\boldsymbol{\theta}}=\hat{\bar{\boldsymbol{\theta}}}}\,.$
(16)

To express the final result, we define the polarization–decomposed FIM
elements as

$\Gamma^{ab}_{ij}\equiv\frac{(f^{a}_{,i}|f^{b}_{,j})}{{\hat{d}_{L}}^{2}}\,,\quad
a,b\in\\{+,\,\crossproduct\\}\,,$ (17)

with $f^{a,b}$ defined by

$\displaystyle
h(\boldsymbol{\theta})=\frac{c_{+}(\iota)}{d_{L}}f^{+}(\bar{\boldsymbol{\theta}})+i\frac{c_{\crossproduct}(\iota)}{d_{L}}f^{\crossproduct}(\bar{\boldsymbol{\theta}})\,,$
(18) $\displaystyle c_{+}(\iota)=\frac{1+\cos^{2}\iota}{2}\,,\quad
c_{\crossproduct}(\iota)=\sin\iota\,.$

where a subscript ,i denotes the partial derivative with respect to
$\bar{\theta}_{i}$. The explicit expression for $f^{a,b}$, which is not
relevant for the discussion here, can be obtained by comparing the above
equations to the full frequency–domain waveform, e.g. Eq. (22) of [33]. We
find that the result can be written as

$\begin{split}&-2\,\log{\cal L}({\cal D}_{\rm
GW}|d_{L}\,,\iota\,,\bar{\boldsymbol{\theta}})\propto\left(\frac{\hat{d}_{L}}{d_{L}}\right)^{2}\Bigg{[}\tilde{\delta}d^{a}_{L}\Gamma^{ab}_{d_{L}d_{L}}\tilde{\delta}d^{b}_{L}\\\
&\quad+(\tilde{\delta}\bar{\theta}^{a}_{i}\
\tilde{\delta}d^{b}_{L}+\tilde{\delta}\bar{\theta}^{b}_{i}\
\tilde{\delta}d^{a}_{L})\Gamma^{ab}_{id_{L}}+\tilde{\delta}\bar{\theta}^{a}_{i}\Gamma_{ij}^{ab}\tilde{\delta}\bar{\theta}^{b}_{j}\Bigg{]}\,,\end{split}$
(19)

with

$\begin{split}&\Gamma^{ab}_{id_{L}}\equiv\frac{(f^{a}_{,i}|f^{b})}{{\hat{d}_{L}}^{3}}\,,\quad\Gamma^{ab}_{d_{L}d_{L}}\equiv\frac{(f^{a}|f^{b})}{{\hat{d}_{L}}^{4}}\,,\\\
&\tilde{\delta}\bar{\theta}_{i}^{+\,,\crossproduct}\equiv\delta\bar{\theta}_{i}\cdot
c_{+\,,\crossproduct}(\iota)\,,\\\
&\tilde{\delta}d_{L}^{+\,,\crossproduct}\equiv\delta d_{L}\cdot
c_{+\,,\crossproduct}(\hat{\iota})-\hat{d}_{L}[c_{+\,,\crossproduct}(\iota)-c_{+\,,\crossproduct}(\hat{\iota})]\,.\end{split}$
(20)

The polarization–decomposed FIM elements
$\Gamma^{ab}_{d_{L}d_{L}}\,,\Gamma^{ab}_{d_{L}i}\,,\Gamma^{ab}_{ij}$ can be
derived from a combination of the usual FIM elements, evaluated at three
conveniently chosen values of $\iota$. Using in particular
$\iota=\\{0\,,\pi/2\,,\pi\\}$, we have
$\Gamma_{IJ}^{++}=4\Gamma_{IJ}(\iota=\pi/2)$,
$\Gamma_{IJ}^{\times\times}=-\left\\{-\left[\Gamma_{IJ}(\iota=0)+\Gamma_{IJ}(\iota=\pi)\right]/2+4\Gamma_{IJ}(\iota=\pi/2)\right\\}$,
$\Gamma_{IJ}^{+\times}+\Gamma_{IJ}^{\times+}=\left[\Gamma_{IJ}(\iota=0)-\Gamma_{IJ}(\iota=\pi)\right]/2$,
with $I(J)$ labeling either $d_{L}$ or $i(j)$ and $i,j$ running over
$\bar{\boldsymbol{\theta}}$. While in this paper we use this strategy, we
point out that a more efficient approach would be to compute the FIM for each
polarization directly, in which case the computational cost would be identical
to that of a single FIM. Using the same argument as in App. B.1, we find that
additional corrections to Eq. (19) at first order in
$\delta\hat{{\theta}}^{i}$ cancel, so evaluating Eq. (19) around the maximum
likelihood point fully encodes the effects of noise fluctuations.

As a further extension, we also report the result for the case
$\boldsymbol{\beta}=\\{d_{L}\,,\iota\,,\psi\\}$, being $\psi$ the polarization
angle. This is relevant for forecasting joint GW+GRB polarization measurements
[73]. In this case, the following definitions are needed

$\Gamma^{a\alpha,b\beta}_{ij}\equiv(f^{a\alpha}_{,i}|f^{b\beta}_{,j})\,,\quad
a,b\in\\{+,\,\crossproduct\\}\,,\quad\alpha,\beta\in\\{c,\,s\\}\,,$ (21)

with

$\displaystyle
h(\boldsymbol{\theta})=\frac{c_{+}(\iota)}{d_{L}}\left[\cos{2\psi}f^{+c}(\bar{\boldsymbol{\theta}})+\sin{2\psi}f^{+s}(\bar{\boldsymbol{\theta}})\right]$
(22) $\displaystyle\qquad\
+i\frac{c_{\crossproduct}(\iota)}{d_{L}}\left[\cos{2\psi}f^{\crossproduct
c}(\bar{\boldsymbol{\theta}})+\sin{2\psi}f^{\crossproduct
s}(\bar{\boldsymbol{\theta}})\right]\,.$

The likelihood Eq. (10) in this approximation can be formally written in a way
that recalls the previous case

$\begin{split}&-2\,\log{\cal L}({\cal D}_{\rm
GW}|d_{L}\,,\iota\,,\,\psi,\,\bar{\boldsymbol{\theta}})\propto\left(\frac{\hat{d}_{L}}{d_{L}}\right)^{2}\Bigg{[}\tilde{\delta}d^{a\alpha}_{L}\Gamma^{a\alpha,b\beta}_{d_{L}d_{L}}\tilde{\delta}d^{b\beta}_{L}\\\
&+(\tilde{\delta}\bar{\theta}^{a\alpha}_{i}\
\tilde{\delta}d^{b\beta}_{L}+\tilde{\delta}\bar{\theta}^{b\beta}_{i}\
\tilde{\delta}d^{a\alpha}_{L})\Gamma^{a\alpha,b\beta}_{id_{L}}+\tilde{\delta}\bar{\theta}^{a\alpha}_{i}\Gamma_{ij}^{a\alpha,b\beta}\tilde{\delta}\bar{\theta}^{b\beta}_{j}\Bigg{]}\,,\end{split}$
(23)

with

$\begin{split}&\tilde{\delta}\bar{\theta}_{i}^{+c}\equiv\delta\bar{\theta}_{i}\cdot
c_{+}(\iota)\cos{2\psi}\\\ &\tilde{\delta}d_{L}^{+c}\equiv\delta d_{L}\cdot
c_{+}(\hat{\iota})\cos{2\hat{\psi}}-\hat{d}_{L}[c_{+}(\iota)\cos{2\psi}-c_{+}(\hat{\iota})\cos{2\hat{\psi}}]\,,\end{split}$
(24)

and all the other values of $\tilde{\delta}\bar{\theta}_{i}^{a\alpha}$ and
$\tilde{\delta}d_{L}^{a\alpha}$ for the indexes $a=\\{+,\crossproduct\\}$ and
$\alpha=\\{c,s\\}$ are obtained by replacing $+$ with $\times$ and/or $\cos
2\psi$ with $\sin 2\psi$ in Eq. (24).

Again, the decomposed FIM matrices can be either computed directly with an
adapted code or derived from the usual FIM, evaluated at specific values of
$\iota$ and $\psi$ (a practical set of values is $\iota=\\{0\,,\pi/2\,,\pi\\}$
and $\psi=\\{0\,,\pi/8\,,\pi/4\\}$).

Finally, we note that this formalism can be extended to all parameters
entering the signal amplitude, in particular the sky position, which is
crucial for obtaining detailed predictions on localization capabilities.

##  Appendix C Population and detector models

Parameter | Description | Fiducial Value
---|---|---
| Cosmology (flat $\bm{\Lambda}$CDM) – $\boldsymbol{\lambda}_{\rm c}$ |
$H_{0}$ | Hubble constant $[$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{Mpc}}^{-1}$]$ | 70
$\Omega_{\rm m,0}$ | Matter energy density | 0.3
| Rate evolution (Madau–Dickinson) – $\boldsymbol{\lambda}_{\rm z}$ |
$\gamma$ | Slope at $z<z_{\rm p}$ | 1.5
$\kappa$ | Slope at $z>z_{\rm p}$ | 5
$z_{\rm p}$ | Peak redshift | 2
| Mass distribution (Gaussian) – $\boldsymbol{\lambda}_{\rm m}$ |
$\mu_{\rm M}$ | Mean $[$\mathrm{M_{\odot}}$]$ | 1.3
$\sigma_{\rm M}$ | Standard deviation $[$\mathrm{M_{\odot}}$]$ | 0.09
| EM detection probability – $\boldsymbol{\lambda}_{\rm EM}$ |
$A_{0}$ | Amplitude of the flux $[$\mathrm{erg}\text{\,}{\mathrm{s}}^{-1}$]$ |  $2.7\text{\times}{10}^{51}$
$\iota_{\rm c}$ | Half–opening of the Gaussian jet profile $[\rm deg]$ | 3.27
$F_{\rm th}$ | GRB detector flux threshold $[$\mathrm{erg}\text{\,}{\mathrm{cm}}^{-2}\text{\,}{\mathrm{s}}^{-1}$]$ |  $2.7\text{\times}{10}^{-07}$
| Single event parameters – $\\{\boldsymbol{\theta}_{\rm i}\\}_{i=1}^{N_{\rm obs}}$ |
$m_{1}^{\rm det}$ | Detector–frame primary mass $[$\mathrm{M_{\odot}}$]$ | –
$m_{2}^{\rm det}$ | Detector–frame secondary mass $[$\mathrm{M_{\odot}}$]$ | –
$d_{L}$ | Luminosity distance $[\rm Gpc]$ | –
$\iota$ | Inclination $[\rm deg]$ | –
Tab. 1: Summary of the population hyperparameters $\boldsymbol{\lambda}$ with
their fiducial values, and of the individual event parameters
$\\{\boldsymbol{\theta}_{\rm i}\\}_{i=1}^{N_{\rm obs}}$ used in the
hierarchical model.

###  C.1 Populations and detectors

For the GW population we assume $p(m_{1},m_{2}|\boldsymbol{\lambda}_{\rm
m})=\mathcal{N}(m_{1}|\mu_{M},\sigma_{M})\,\mathcal{N}(m_{2}|\mu_{M},\sigma_{M})$
with $\mathcal{N}(x|\mu_{M},\sigma_{M})$ a Gaussian distribution with fiducial
mean and standard deviation $\\{\mu_{M,{\rm
fid}}=$1.3\text{\,}\mathrm{M_{\odot}}$,\,\sigma_{M,{\rm
fid}}=$0.09\text{\,}\mathrm{M_{\odot}}$\\}$ [74]. The merger rate evolution
with redshift in Eq. (5) is modeled as a Madau–Dickinson (MD) profile [55],
$\psi(z;\boldsymbol{\lambda}_{\rm
z})=(1+z)^{\gamma}/\\{1+[(1+z)/(1+z_{\mathrm{p}})]^{\gamma+\kappa}\\}$. We use
fiducial parameters $\boldsymbol{\lambda}_{\rm z}=\\{\gamma_{\rm
fid}=1.5,\,\kappa_{\rm fid}=5,\,z_{{\rm p},{\rm fid}}=2\\}$. The fiducial
values we use are obtained by convolving a MD profile with parameters
determined by the Star Formation Rate [55] with a time–delay distribution
chosen to be log–flat. We refer to App. A of [33] for details. Inclinations
are drawn from a distribution uniform on the sphere, i.e. flat in $\cos\iota$.
Since this is the same prior used for the analysis of individual GW events,
there is no need to include explicitly this term in the hierarchical
likelihood since it will simplify between numerator and denominator [26].
Finally, cosmology is assumed to be described by a flat $\Lambda$CDM model
with
$H_{0}=$70\text{\,}\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{Mpc}}^{-1}$$,
$\Omega_{\rm m,0}=0.3$.

We consider an observing scenario representative of the best case for the
current generation of GW detectors, with a LIGO–Virgo–KAGRA–LIGO India network
composed of two LIGO detectors in the US, the Virgo and KAGRA detectors in
Italy and Japan, respectively, and a LIGO detector located in India, with A$+$
sensitivity.131313Specifically, we use the AplusDesign public sensitivity
curve for the three LIGO detectors, avirgo_O5low_NEW for Virgo, and
kagra_80Mpc for KAGRA [53], available at
https://dcc.ligo.org/LIGO-T2000012/public. We assume a $100\%$ duty cycle and
a SNR threshold of $\rho_{\rm obs}\geq 12$, where the “observed” SNR is
simulated encoding the effects of detector noise.

For the EM detection model, we assume that the jet has a Gaussian–shaped
profile, $E(\iota)=E_{0}\exp{-\left(\iota^{2}/\iota_{\rm c}^{2}\right)/2}$,
and we model the detection probability with a threshold on the observed
$\gamma$–ray flux of each event.141414We assume the uncertainty in the
measurement of the flux to be negligible. This is defined as $F_{\rm
GRB}=\eta\,E(\iota)/(4\pi d_{L}^{2}T_{90})$, with $\eta=0.1$ being the
radiative efficiency [75, 76], and $T_{90}$ the duration of the period during
which $90\%$ of the burst’s energy is emitted, which we assume to be
$2\text{\,}\mathrm{s}$. Therefore, we define $A_{0}=\eta E_{0}/T_{90}$, and we
have explicitly

$F_{\rm GRB}=\frac{A_{0}}{4\pi
d_{L}^{2}}\,e^{-\nicefrac{{\iota^{2}}}{{2\iota_{c}^{2}}}}\,.$ (25)

We assume as fiducial values $E_{0,{\rm
fid}}=$10^{52.73}\text{\,}\mathrm{erg}$,\;\iota_{\rm c,{\rm
fid}}=$3.27\text{\,}\mathrm{\SIUnitSymbolDegree}$$
($0.057\text{\,}\mathrm{rad}$), obtained by a multimessenger analysis of
GW170817 [56]. We will consider an event as detected if the flux exceeds the
Fermi–GBM–like sensitivity of $F_{\rm
th}=$2\text{\times}{10}^{-7}\text{\,}\mathrm{erg}\text{\,}{\mathrm{cm}}^{-2}\text{\,}{\mathrm{s}}^{-1}$$
[54].151515This model corresponds to the one used in [77, 78], with the slight
simplification that they introduce a distribution on the peak luminosity
$E_{0}$, while here we consider it fixed. We have checked that the scaling of
$\Theta_{\rm max}$ with redshift is consistent among the two. In conclusion,
we define $\boldsymbol{\lambda}_{\rm EM}=\\{A_{0},\iota_{\rm c},F_{\rm
th}\\}$, and we have

$P({\rm det}_{\rm EM}|\iota,d_{L},A_{0},\iota_{\rm
c})=\Bigg{\\{}\begin{array}[]{c l}1&{\rm if}\ F_{\rm
GRB}(\iota,d_{L},A_{0},\iota_{\rm c})\geq F_{\rm th}\\\ 0&{\rm
otherwise}\end{array}\,.$ (26)

Tab. 1 summarizes the population hyperparameters and individual event
parameters used in the paper.

###  C.2 Mock dataset

We simulate the GW likelihood for each joint detection by using the likelihood
approximant described in App. B. We use the public software gwfast [41, 33]
for the evaluation of the FIM, and we extract samples from the posterior given
by Eq. (19) multiplied by a prior (denoted by
$\pi(\boldsymbol{\theta}_{i})\,\forall i$) flat in $\cos(\iota)$ with
$\iota\in[0,\pi]$, flat in detector–frame masses, and flat in the other
variables. This reproduces actual data–analysis choices. We use the
affine–invariant MCMC sampler zeus [79, 80]. We assume that the sky position
is known and fixed to the host galaxy, and the GPS time of detection is
known.161616To limit the computational cost, we also fix the phase at
coalescence, polarization angle, and tidal deformabilities. We checked that
this does not lead to underestimation of the uncertainty on the masses,
distance, and inclination. On the contrary, fixing the spins or neglecting
them would lead to substantial underestimation of the uncertainty on the
detector–frame mass.

For the EM data, we assume a Gaussian likelihood with standard deviation
$\sigma_{z}=10^{-3}$. We simulate a redshift measurement by drawing, for each
event that passes the selection cut in the GW and EM, an “observed” value from
this likelihood.

In summary, this procedure yields a set of samples from the posterior
$p(\boldsymbol{\theta}_{i}|{\cal D}_{\rm
GW}^{i})\propto\mathcal{L}\big{(}{\cal D}_{\rm
GW}^{i}|\boldsymbol{\theta}_{i}\big{)}/\pi(\boldsymbol{\theta}_{i})$ for each
GW event and a redshift measurement, as in data analysis of real events.

Finally, to estimate the selection effect $P({\rm det}|\boldsymbol{\lambda})$
in Eq. (1) we use reweighted Monte Carlo integration [81, 82]. We generate a
large injection set and apply the same selection cut used for constructing the
catalog for GW observations, while crucially no EM selection is applied to the
injections, as the EM selection probability will be inferred during inference.

###  C.3 Hierarchical inference

Starting from the mock data, we sample Eq. (1) in the high–dimensional space
of $\boldsymbol{\lambda},\\{\boldsymbol{\theta}_{i}\\}$ using PyMC [57] and
the Hamiltonian Monte Carlo [58] based scheme NUTS [59]. This choice is
dictated by the necessity of accurately exploring the posterior distribution
of individual GW events in the corner of the parameter space at small
inclination which is preferentially selected by the EM detection probability.
Using standard approaches that rely on Monte Carlo integrals over existing
posterior samples [26] would either result in numerically unreliable results
or require an unreasonably large amount of individual posterior samples. In
our implementation, prior to inference we instead fit a continuous interpolant
of the individual GW likelihoods with a Gaussian Mixture Model over posterior
samples [83], from which we can continuously resample at inference time.