Cosmology and Astrophysics with Standard Sirens and Galaxy Catalogs in View of Future Gravitational Wave Observations

The core of the method is the construction of a population prior on the GW source redshift from a galaxy catalog. This is given by the second term of Equation (8) that we factorize as follows:

$$\begin{eqnarray}&&p(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}},{{\boldsymbol{\lambda }}}_{{\rm{z}}})=\displaystyle \frac{{p}_{\mathrm{gal}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}}){p}_{\mathrm{rate}}(z\,| \,{{\boldsymbol{\lambda }}}_{{\rm{z}}})}{\int {p}_{\mathrm{gal}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}}){p}_{\mathrm{rate}}(z\,| \,{{\boldsymbol{\lambda }}}_{{\rm{z}}}){dz}\,d\hat{{\rm{\Omega }}}},\end{eqnarray} \tag{ 9 }$$

where p_gal is the probability that there is a galaxy at ($z,\hat{{\rm{\Omega }}}$) and p_rate the probability of a galaxy at redshift z to host a GW event. This takes into account that the probability for a galaxy to host a merger can have a nontrivial redshift dependence. We assume that this can be parameterized as

$$\begin{eqnarray}&&{p}_{\mathrm{rate}}(z\,| \,{{\boldsymbol{\lambda }}}_{{\rm{z}}})\propto \displaystyle \frac{\psi (z;{{\boldsymbol{\lambda }}}_{{\rm{z}}})}{(1+z)}\end{eqnarray} \tag{ 10 }$$

where ψ(z; λ _z) is the merger rate evolution of compact objects with redshift and the term (1 + z)⁻¹ takes into account the conversion between source and detector time. Note that the normalization integral in Equation (9), as well as any rate normalization factor in Equation (10), simplifies between numerator and denominator in Equation (1) and does not need to be computed explicitly.

As a first approximation, we assume to have a complete galaxy catalog, i.e., this contains all the potential host galaxies. In this case, we denote as ${p}_{\mathrm{cat}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}})$ the probability distribution built from galaxies in the catalog. We can write ${p}_{\mathrm{gal}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}})={p}_{\mathrm{cat}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}})$ and compute the latter probability as a sum over the contribution of each galaxy. Given a set of ${{\boldsymbol{d}}}^{\mathrm{EM}}=\{{{\boldsymbol{d}}}_{g}^{\mathrm{EM}}\}$ (g = 1,...,N_gal) EM observations of galaxies we have

$$\begin{eqnarray}&&{p}_{\mathrm{cat}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}})=\displaystyle \frac{{\sum }_{g}{w}_{g}\,p(z\,| \,{{\boldsymbol{d}}}_{g}^{\mathrm{EM}},{{\boldsymbol{\lambda }}}_{{\rm{c}}})\delta (\hat{{\rm{\Omega }}}-{\hat{{\rm{\Omega }}}}_{g})}{{\sum }_{g}{w}_{g}}\end{eqnarray} \tag{ 11 }$$

where w_g weights the probability of each galaxy to host a GW event (e.g., by the galaxy luminosity), $p(z\,| \,{{\boldsymbol{d}}}_{g}^{\mathrm{EM}},{{\boldsymbol{\lambda }}}_{{\rm{c}}})$ is the galaxy's redshift distribution that we want to use as a prior, and $\delta (\hat{{\rm{\Omega }}}-{\hat{{\rm{\Omega }}}}_{g})$ is a Dirac delta distribution of each galaxy's sky localization that can be treated as errorless.

The galaxy catalog contains redshift measurements ${\tilde{z}}_{g}$ and associated uncertainties ${\tilde{\sigma }}_{z,g}$ for each observed galaxy, i.e., ${{\boldsymbol{d}}}_{g}^{\mathrm{EM}}\,=\{{\tilde{z}}_{g},{\tilde{\sigma }}_{z,g}\}$. From these quantities, we construct the likelihoods, which we assume to be Gaussian, $p({\tilde{z}}_{g}\,| \,z)={ \mathcal N }(z;{\tilde{z}}_{g},{\tilde{\sigma }}_{z,g}^{2})$. This is a probability distribution over the observed values ${\tilde{z}}_{g}$. To get $p(z\,| \,{{\boldsymbol{d}}}_{g}^{\mathrm{EM}},{{\boldsymbol{\lambda }}}_{{\rm{c}}})$ we need to multiply it by a prior on the redshift distribution, which in the absence of other information is naturally chosen as uniform in comoving volume (Gair et al. 2023). Using Bayes' theorem, we get

$$\begin{eqnarray}&&p(z\,| \,{{\boldsymbol{d}}}_{g}^{\mathrm{EM}},{{\boldsymbol{\lambda }}}_{{\rm{c}}})=\displaystyle \frac{{ \mathcal N }(z;{\tilde{z}}_{g},{\tilde{\sigma }}_{z,g}^{2})\tfrac{{{dV}}_{c}}{{dz}}(z;{{\boldsymbol{\lambda }}}_{{\rm{c}}})}{\int { \mathcal N }(z;{\tilde{z}}_{g},{\tilde{\sigma }}_{z,g}^{2})\tfrac{{{dV}}_{c}}{{dz}}(z;{{\boldsymbol{\lambda }}}_{{\rm{c}}}){dz}},\end{eqnarray} \tag{ 12 }$$

where dV_c /dz is the differential comoving volume element in a flat universe. With this definition, Equation (11) is normalized so that if $p({\tilde{z}}_{g}\,| \,z)=\delta (z-{\tilde{z}}_{g})$ and in case of uniform weights, we get ${n}_{\mathrm{cat}}(z)=(1/{V}_{c}){\sum }_{g}\delta (z-{\tilde{z}}_{g})$, which is consistent with Equation (3.39) of Finke et al. (2021), i.e., the comoving density of galaxies n_cat(z) is estimated by counting the objects in the catalog and dividing by the total volume. If, instead, the likelihood is completely uninformative, $p({\tilde{z}}_{g}\,| \,z)=1$, we correctly get a comoving density constant in redshift. Note that the assumption of a uniform-in-comoving-volume prior introduces a dependence on the cosmology in the redshift prior, since dV_c /dz is a function of λ _c. Once normalizing correctly $p(z\,| \,{{\boldsymbol{d}}}_{g}^{\mathrm{EM}},{{\boldsymbol{\lambda }}}_{{\rm{c}}})$, however, the dependence on H₀ cancels, as this is an overall multiplicative factor in dV_c /dz, but a dependence on Ω_m,0 remains and should be accounted for. The same would be true if one wished to include more general cosmologies, e.g., a Dark Energy equation of state w_DE ≠ −1. Alternatively, one may adopt a flat prior, which eliminates the extra dependence on the cosmology, as done, e.g., in Gray et al. (2023).

In the real world, observed galaxy catalogs are subject to selection effects that result in having an incomplete set of potential hosts. To address this issue, the p_gal term in Equation (9) has to be modified including the probability of missing galaxies, ${p}_{\mathrm{miss}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}})$, as follows (Chen et al. 2018; Finke et al. 2021):

$$\begin{eqnarray}&&\begin{array}{l}{p}_{\mathrm{gal}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}})={f}_{{ \mathcal R }}\,{p}_{\mathrm{cat}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}})\\ \quad +\ (1-{f}_{{ \mathcal R }}){p}_{\mathrm{miss}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}}).\end{array}\end{eqnarray} \tag{ 13 }$$

The probability ${p}_{\mathrm{miss}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}})$ is constructed with two pieces of information. The first is the number of missing galaxies as a function of redshift and sky position, encoded in the completeness fraction of the catalog ${P}_{\mathrm{comp}}(z,\hat{{\rm{\Omega }}})$. This can be estimated with some assumption on the total comoving density of galaxies or on the luminosity distribution of a complete sample that is assumed to follow a Schechter function. We refer to Finke et al. (2021) for a rigorous definition of this quantity and for the discussion of different ways to compute it. In CHIMERA, we will use in particular the so-called mask completeness of Finke et al. (2021), which we will describe better in Section 2.4. From this, one can compute the weight ${f}_{{ \mathcal R }}$ in Equation (13) as

$$\begin{eqnarray}&&{f}_{{ \mathcal R }}\equiv \displaystyle \frac{1}{{V}_{c}({{\boldsymbol{\lambda }}}_{{\rm{c}}})}\int {P}_{\mathrm{comp}}(z,\hat{{\rm{\Omega }}})\,{{dV}}_{c},\end{eqnarray} \tag{ 14 }$$

where the integral extends to a region sufficiently large to encompass all the GW events under consideration and V_c ( λ _c) denotes the corresponding comoving volume.

The second bit of information amounts to specifying how the missing galaxies are distributed. Two extreme possibilities are a homogeneous completion, where the distribution is assumed to be uniform in comoving volume and sky position, and a multiplicative completion, where missing galaxies are assumed to trace the distribution of those present in the catalog (Finke et al. 2021). The two correspond to

$$\begin{eqnarray}&&{p}_{\mathrm{miss}}^{\mathrm{HOM}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}})=\displaystyle \frac{1-{P}_{\mathrm{comp}}(z,\hat{{\rm{\Omega }}})}{(1-{f}_{{ \mathcal R }}){V}_{c}({{\boldsymbol{\lambda }}}_{{\rm{c}}})}\,\displaystyle \frac{{{dV}}_{c}}{{dz}}(z;{{\boldsymbol{\lambda }}}_{{\rm{c}}}),\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&{p}_{\mathrm{miss}}^{\mathrm{MULT}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}})=\displaystyle \frac{{f}_{{ \mathcal R }}}{1-{f}_{{ \mathcal R }}}\,\displaystyle \frac{1-{P}_{\mathrm{comp}}(z,\hat{{\rm{\Omega }}})}{{P}_{\mathrm{comp}}(z,\hat{{\rm{\Omega }}})}\,{p}_{\mathrm{cat}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}}).\end{eqnarray} \tag{ 16 }$$

In practice, as suggested in Finke et al. (2021) and done in the code DarkSirensStat, it is convenient to interpolate between a multiplicative completion where the catalog is fairly complete, and a homogeneous completion where the catalog is largely incomplete. Finally, an even better option, suggested in Finke et al. (2021) and recently studied in detail in Dalang & Baker (2024), would be to complete the catalog using the typical correlation length of galaxies.

By putting together Equations (1)–(14) we can write the likelihood in its final form. In the following, we explicitly restrict to the subset of source parameters $\bar{{\boldsymbol{\theta }}}=\{z,\hat{{\rm{\Omega }}},{m}_{1},{m}_{2}\}$ and ${\bar{{\boldsymbol{\theta }}}}^{\det }=\{{d}_{L},\hat{{\rm{\Omega }}},{m}_{1}^{\det },{m}_{2}^{\det }\}$. This relies on the assumption that the remaining parameters of the GW waveform (see Equation (27) below) have a population distribution that coincides with the prior used in the analysis of individual events. We have

$$\begin{eqnarray}&&\begin{array}{l}p({{\boldsymbol{d}}}^{\mathrm{GW}}\,| \,{\boldsymbol{\lambda }})\propto \displaystyle \frac{1}{\xi {\left({\boldsymbol{\lambda }}\right)}^{{N}_{\mathrm{ev}}}}\displaystyle \prod _{i=1}^{{N}_{\mathrm{ev}}}\int {dz}\,d\hat{{\rm{\Omega }}}\\ \times \,{{ \mathcal K }}_{\mathrm{gw},i}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}},{{\boldsymbol{\lambda }}}_{{\rm{m}}})\,{p}_{\mathrm{gal}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}})\displaystyle \frac{\psi (z;{{\boldsymbol{\lambda }}}_{{\rm{z}}})}{1+z},\end{array}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}&&\begin{array}{l}{{ \mathcal K }}_{\mathrm{gw},i}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}},{{\boldsymbol{\lambda }}}_{{\rm{m}}})\equiv \displaystyle \int {{dm}}_{1}{{dm}}_{2}\\ \times \,\displaystyle \frac{p(z,{m}_{1},{m}_{2},\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{d}}}_{i}^{\mathrm{GW}},{{\boldsymbol{\lambda }}}_{{\rm{c}}})}{\pi ({d}_{L})\pi ({m}_{1}^{\det })\pi ({m}_{2}^{\det })}\,\displaystyle \frac{1}{\tfrac{{{dd}}_{L}}{{dz}}(z;{{\boldsymbol{\lambda }}}_{{\rm{c}}}){\left(1+z\right)}^{2}}\\ \,\times \,p({m}_{1},{m}_{2}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{m}}}),\end{array}\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&\begin{array}{l}\xi ({\boldsymbol{\lambda }})=\displaystyle \int d{\boldsymbol{\theta }}\,{P}_{\det }({\boldsymbol{\theta }},{{\boldsymbol{\lambda }}}_{{\rm{c}}})\\ \times \,p({m}_{1},{m}_{2}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{m}}})\,{p}_{\mathrm{gal}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}})\displaystyle \frac{\psi (z;{{\boldsymbol{\lambda }}}_{{\rm{z}}})}{1+z}.\end{array}\end{eqnarray} \tag{ 19 }$$

Note that we have omitted the overall normalization of the redshift prior, as this is independent of $z,\hat{{\rm{\Omega }}}$ and is present also in the selection function, leading to a simplification between numerator and denominator in Equation (17). A different normalization of this quantity just amounts to a different interpretation of the selection function, which does not correspond to the fraction of detectable events, as a consequence of the fact that p_pop loses its correct normalization to unity. We find the form of Equations (17)–(19) useful since it makes clear what quantities actually need to be computed in the numerical evaluation. On the other hand, we note that the normalization of the term p_gal has to be computed carefully, in order to correctly balance the contributions of the catalog and missing terms in Equation (13).

We now discuss the implementation of the method in CHIMERA. The code accurately computes Equation (17) across different regimes, i.e., when the EM redshift information in highly constraining (bright sirens), more loosely constraining (dark sirens), or not present at all (spectral sirens). We also prioritize computational efficiency to meet the demands of next-generation GW observatories and upcoming galaxy surveys.

In this section, we summarize the key aspects of the code, while a more detailed description is provided in Appendix A. We will use the real event GW170817 as an illustrative example, considering it a dark siren in conjunction with the GLADE+ galaxy catalog. This is useful because GW170817 is the best-localized GW event so far and its host galaxy, NGC 4993, has been uniquely identified (Abbott et al. 2017a), providing an ideal benchmark for dark siren techniques. Indeed, this is often used as a test case and allows comparison to other pipelines (Mastrogiovanni et al. 2023,  Gray et al. 2023). In Appendix B, we provide details about the data and cuts used for this example.

Figure 1 shows a compact representation of the workflow of CHIMERA, in the case of GW170817. In short, the galaxy catalog term is precomputed as a function of redshift (red lines) on a pixel-by-pixel basis inside the GW localization region (given by the gray squares). The GW kernel in Equation (18) is given by a kernel density estimation (KDE) of the GW posterior samples reweighted for the population prior on source-frame masses (blue area).

Zoom In Zoom Out Reset image size

Pixelization and precomputation of the galaxy term. The quantity ${p}_{\mathrm{gal}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}})$ can be cosmology dependent, but only when adopting a cosmology-dependent prior on the galaxies' redshifts, and in any case never on H₀, as explained in Section 2. If adopting a flat prior on the galaxies' redshifts, or a narrow prior on the other cosmological parameters (e.g., a Planck prior on the matter density), the quantity p_gal can be computed once for all, given a specific galaxy catalog. For each event in the GW catalog, CHIMERA selects galaxies enclosed in the GW localization volume, with a redshift range computed taking into account all of the prior range for H₀. Then, the sky area of the GW event is divided into equal-area pixels in the pixelization scheme of Healpix (Górski et al. 2005; Zonca et al. 2019). CHIMERA includes an adaptive pixelization procedure that allows to set the pixel size for each GW event depending on its localization area (see Appendix A for details). This step corresponds to the left panel of Figure 2, where we show the galaxies enclosed in the 90% localization volume of GW170817 as red crosses, overlapping the pixels. The quantity ${p}_{\mathrm{gal}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}})$ is then approximated as a redshift-dependent function of the form given in Equations (11) and (12) for each pixel, by summing the contribution of all the galaxies enclosed within it. The weights in Equation (11) can be chosen to be either uniform or corresponding to the galaxies' luminosity in the catalog. The redshift distributions in each pixel are shown in red in the central panel of Figure 2 for the example of GW170817. This method allows for the compression of the information coming from all galaxies in each pixel in a single function of redshift. In the case of GW170817, we find 88 galaxies in the 90% localization volume, while using 20 pixels, which corresponds to a factor of ∼5 improvement in the computational cost. This procedure is even more effective in the case of events with thousands of galaxies in the localization region. We note that, in general, the GW sky map changes on much larger angular scales than the galaxy distribution. This means that the compression introduced here (which basically amounts to approximating the galaxies' angular positions with the center of each pixel) does not lead to a loss of information as far as the pixel size is chosen appropriately.

Zoom In Zoom Out Reset image size

Catalog completeness. The completeness fraction ${P}_{\mathrm{comp}}(z,\hat{{\rm{\Omega }}})$ is computed following the mask method presented in Finke et al. (2021). First, similar pixels are grouped together in masks by applying an agglomerative clustering algorithm. The chosen feature for this purpose is the number of galaxies within each pixel. The number of masks is arbitrarily chosen depending on the survey properties. For each mask, the completeness fraction is computed by comparing the luminosity function at each redshift to a reference Schechter function. We refer to Finke et al. (2021) for a more detailed description of this procedure. In particular, this allows us to easily cut regions particularly under-covered by the galaxy survey, for example, the region of the sky covered by the Milky Way. It also makes straightforward use of galaxy catalogs with partial sky coverage. In the example of GW170817, it is complete.

3D GW kernel. We now turn to the term ${{ \mathcal K }}_{\mathrm{gw}}$ in Equation (18). For each GW event, we have a set of posterior samples of the detector frame quantities ${\bar{{\boldsymbol{\theta }}}}^{\det }\,=\{{d}_{L},\hat{{\rm{\Omega }}},{m}_{1}^{\det },{m}_{2}^{\det }\}$ approximating the posterior probability $p({m}_{1}^{\det },{m}_{2}^{\det },{d}_{L},\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{d}}}_{i}^{\mathrm{GW}})$. CHIMERA assumes that the prior on detector frame masses used for the individual event analysis is flat, $\pi ({m}_{1}^{\det })=\pi ({m}_{2}^{\det })=1$, while the prior on the luminosity distance is $\pi ({d}_{L})\propto {d}_{L}^{2}$. Then, for each value of λ , the samples are converted to source-frame quantities by inverting Equations (3)–(5). Each sample, labeled by j, is assigned a weight equal to

$$\begin{eqnarray}&&{w}_{{ \mathcal K },j}({\boldsymbol{\lambda }})=\displaystyle \frac{p({m}_{1,j},{m}_{2,j}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{m}}})}{{d}_{L}^{2}\,\tfrac{{{dd}}_{L}}{{dz}}({z}_{j};{{\boldsymbol{\lambda }}}_{{\rm{c}}}){\left(1+{z}_{j}\right)}^{2}},\end{eqnarray} \tag{ 20 }$$

and ${{ \mathcal K }}_{\mathrm{gw}}$ is obtained by a 3D weighted KDE in the space (z, R.A., decl.) with weights ${w}_{{ \mathcal K },j}({\boldsymbol{\lambda }})$. Using the reweighted samples in this subspace corresponds to performing the integral in Equation (18) with the correct population prior p(m₁, m₂ ∣ λ _m) on the source-frame masses. This will encode a dependence on the cosmological parameters, which enters in the conversion of the posterior samples from the detector frame to the source frame, resulting in a reshaping of the kernel ${{ \mathcal K }}_{\mathrm{gw}}$. A nonflat source-frame mass spectrum will add redshift information at the statistical level through this term, even in the absence of a galaxy catalog, in which case we recover the spectral siren case. For GW170817, the GW kernel is plotted in blue in Figure 2(central panel), assuming a value of H₀ = 70 km s⁻¹ Mpc⁻¹ and a flat mass distribution between 1 and 3 M_☉

Integration. The double integral in Equation (17) is performed by integrating numerically dz in each pixel first. Then, the integral in $d\hat{{\rm{\Omega }}}$ can be approximated as a discrete sum of the results in each pixel, multiplied by the pixel area. The latter, for each given GW event, is equal for all pixels and inversely proportional to the number of Healpix pixels of the map. So, in the end, the angular integral in Equation (17) is obtained as an average over pixels.

Selection effects. The selection bias term ξ( λ ) from Equation (2) is computed by using the standard reweighted Monte Carlo (MC) method (Tiwari 2018; Farr 2019). In this case, it is crucial to use detector frame quantities, since the detectability in a GW experiment is a function of those only, i.e., ${P}_{\det }={P}_{\det }({{\boldsymbol{\theta }}}^{\det })$. This allows us to compute this function only once. The integral in Equation (19) can then be computed with the change of variables ${\boldsymbol{\theta }}\mapsto {{\boldsymbol{\theta }}}^{\det }({\boldsymbol{\theta }},{{\boldsymbol{\lambda }}}_{{\rm{c}}})$, where the Jacobian factors are given in Equation (7). CHIMERA takes as input a set of simulated events, each with detector frame parameters ${{\boldsymbol{\theta }}}_{i}^{\det }$, drawn from a reference distribution with probability ${p}_{\mathrm{draw}}({{\boldsymbol{\theta }}}_{i}^{\det })$. These must have been previously injected in the same detection pipeline used to obtain the GW catalog, storing those that pass the detection threshold.

Then, the selection function in Equation (19) is computed with MC integration as

$$\begin{eqnarray}\begin{array}{rcl}\xi ({\boldsymbol{\lambda }}) & = & \displaystyle \frac{1}{{N}_{\mathrm{inj}}}\displaystyle \sum _{i=1}^{{N}_{\det }}\displaystyle \frac{1}{{p}_{\mathrm{draw}}({{\boldsymbol{\theta }}}_{i}^{\det })}\displaystyle \frac{1}{\tfrac{{{dd}}_{L}}{{dz}}({z}_{i},{{\boldsymbol{\lambda }}}_{{\rm{c}}}){\left(1+{z}_{i}\right)}^{2}}\\ & & \times \ p({m}_{1,i},{m}_{2,i}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{m}}})\,{p}_{\mathrm{gal}}({z}_{i},{\hat{{\rm{\Omega }}}}_{i}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}})\displaystyle \frac{\psi ({z}_{i};{{\boldsymbol{\lambda }}}_{{\rm{z}}})}{1+{z}_{i}},\end{array}\end{eqnarray} \tag{ 21 }$$

where N_inj is the total number of injections, N_det the number of detected ones, and the source-frame quantities z_i , m_1,i , and m_2,i are understood to be computed from the (fixed) detector frame ones via inversion of Equations (3)–(5), and are thus functions of the cosmological parameters λ _c. Finally, when computing the selection bias term, CHIMERA employs a galaxy catalog interpolant p_gal constructed on the full sky.

Accuracy of the MC integrals. Finally, when computing Equation (21), we ensure that the effective number of independent draws after reweighting, N_eff, for a given population is sufficiently high (${N}_{\mathrm{eff}}\gt 5{N}_{\det }$, see Farr 2019). This criterion is also applied to the GW kernel weights to ensure numerical stability.

In conclusion, the right panel of Figure 2 shows the results for the posterior probability on H₀ (solid line). We obtain a value of ${H}_{0}={73}_{-22}^{+58}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$ in the dark siren case, in good agreement with Fishbach et al. (2019). As can be seen from Figure 2 (central panel), the primary contribution to the posterior arises from galaxies at approximately z ∼ 0.01. This value combined with the GW measurement of d_L ∼ 40 Mpc, implies H₀ ∼ 70 km s⁻¹ Mpc⁻¹. The galaxy groups at z ∼ 0.02 (0.04) provide a shallower contribution to H₀ ∼ 150 (300)km s⁻¹ Mpc⁻¹ due to the presence of selection effects that disfavor high values of H₀ $[\xi ({H}_{0})\sim {H}_{0}^{3}$]. In this case, the assumption of a population model for GW170817 provides no significant improvements to H₀.

For comparison and to test the performance of the code in the bright siren regime, we repeat the analysis assuming the identification of the host galaxy with z = 0.0100 ± 0.0005 (Abbott et al. 2017a). In this case, we assign a weight w = 1 to NGC 4993 and w = 0 to all the other galaxies (see Equation (11)). The resulting posterior on H₀ is shown in Figure 2(c) (dashed line). In the bright siren case, we obtain a value of ${H}_{0}={69}_{-8}^{+15}\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{\mathrm{Mpc}}^{-1}$, in very good agreement with Abbott et al. (2023b).

We generate our mock galaxy catalog (hereafter, parent sample) from the MICE Grand Challenge light-cone simulation (v2), which populates one octant of the sky (close to 5157 deg²) designed to mimic a complete DES-like survey up to an observed magnitude of i < 24 at redshift z < 1.4 (Carretero et al. 2015; Crocce et al. 2015; Fosalba et al. 2015a, 2015b; Hoffmann et al. 2015). MICE assumes a flat ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m,0 = 0.25, and Ω_Λ,0 = 0.75.

While we ideally require a complete catalog with high number density, as a simplifying assumption in this paper we consider only galaxies with stellar masses $\mathrm{log}{M}_{\star }/{M}_{\odot }\gt 10.5$. This cut is consistent with the idea that the binary merger rate is traced by stellar mass, as also adopted in current standard sirens analysis via absolute magnitude cuts and luminosity weighting (Fishbach et al. 2019; Finke et al. 2021; Gray et al. 2022; Abbott et al. 2023b; Mastrogiovanni et al. 2023; Gray et al. 2023). A similar cut in mass is also considered (Muttoni et al. 2023) in the context of simulations for the Einstein Telescope.

We subsample the MICEv2 catalog to reproduce the density for the cut described above extracting the galaxies to get a uniform in comoving volume distribution. In the end, we obtain a parent sample of about 1.6 million massive galaxies.

For the redshift uncertainties, we consider two cases. First of all, we explore the possibility of maximizing the galaxy catalog information by having a spectroscopic catalog. This could be done by expanding the currently available catalogs (GLADE+ Dálya et al. 2022) in the future by exploiting the information provided by the next large spectroscopic surveys. As an example, the ESA mission Euclid (Laureijs et al. 2011) will provide an all-sky map of spectroscopic redshift in the range of 0.9 < z < 1.8, with an accuracy of σ_z /(1 + z) ≲ 0.001, and the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration et al. 2016) is planned to observe ∼14,000 deg² covering the redshift range of 0.4 < z < 2.1. Second, we study how the information extracted changes when using photometric redshift, assuming an uncertainty σ_z /(1 + z) = 0.05. This is easily accessible with current ongoing surveys like DES (e.g., DES has reached σ_z ∼ 0.01 (Myles et al. 2021), and this limit can be pushed to σ_z ∼ 0.007 with improved techniques, e.g., Buchs et al. 2019) on a smaller area, and in future surveys like Euclid and Rubin Observatory (Ivezić et al. 2019) are planned to extend it to the entire sky to a depth of Euclid H magnitude H_E ∼ 24 (H_E ∼ 26 in the Deep Survey) and an expected uncertainty σ_z /(1 + z) ≲ 0.05 (Euclid Collaboration et al. 2020, 2022).

Therefore, in this paper, we consider two regimes of photometric (hereafter, "phot") and spectroscopic (hereafter, "spec") redshift uncertainties:

$$\begin{eqnarray}{\sigma }_{z}=\left\{\begin{array}{ll}0.001(1+z) & ({z}_{\mathrm{spec}})\\ 0.05(1+z) & ({z}_{\mathrm{phot}})\end{array}\right..\end{eqnarray} \tag{ 22 }$$

We generate mock GW events from the parent sample by fixing cosmological hyperparameters λ _c and astrophysical population hyperparameters λ _z and λ _m. We describe the redshift and mass distributions in turn.

3.2.1. Rate Evolution

The source-frame merger rate is parameterized as follows (see Madau & Dickinson 2014):

$$\begin{eqnarray}&&\psi (z;{{\boldsymbol{\lambda }}}_{{\rm{z}}})=\displaystyle \frac{{\left(1+z\right)}^{\gamma }}{1+{\left(\tfrac{1+z}{1+{z}_{{\rm{p}}}}\right)}^{\gamma +\kappa }},\end{eqnarray} \tag{ 23 }$$

where ψ(z) ∝ (1 + z)^γ at low z, then reaches its peak near z_p, and subsequently declines as ψ(z) ∝ (1 + z)^−κ . To build the catalog, we use γ = 2.7 consistent with the LVK GWTC-3 results (Abbott et al. 2023c). The limited detection range of current GW detectors restricts our ability to determine the merger rate at higher redshifts, therefore we adopt z_p = 2, κ = 3, consistent with the idea that ψ(z; λ _z) follows the galaxy's star formation rate density with parameters from Madau & Dickinson (2014). The catalog of potential sources is then obtained by sampling the parent sample using a weight proportional to the detector frame merger rate ψ(z; λ _z)/(1 + z).

3.2.2. Mass Distribution

For the mass distribution, we adopt the phenomenological "PowerLaw+Peak" (PLP) model following the LVK GWTC-3 results (Abbott et al. 2023c). The mass term of Equation (8) is factorized as follows:

$$\begin{eqnarray}&&p({m}_{1},{m}_{2}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{m}}})=p({m}_{1}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{m}}})p({m}_{2}\,| \,{m}_{1},{{\boldsymbol{\lambda }}}_{{\rm{m}}}).\end{eqnarray} \tag{ 24 }$$

The probability of the primary BH mass is given by

$$\begin{eqnarray}&&p({m}_{1}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{m}}})\propto \left[(1-{\lambda }_{{\rm{g}}}){ \mathcal P }({m}_{1})+{\lambda }_{{\rm{g}}}\,{ \mathcal G }({m}_{1})\right]{ \mathcal S }({m}_{1}),\end{eqnarray} \tag{ 25 }$$

where ${ \mathcal P }(m)\propto {m}^{-\alpha }$ is a power law truncated in the domain m ∈ [m_low, m_high], ${ \mathcal G }(m)\propto { \mathcal N }({\mu }_{{\rm{g}}};{\sigma }_{{\rm{g}}}^{2})$ is a Gaussian component whose contribution relative to ${ \mathcal P }$ is regulated by λ_g, and ${ \mathcal S }(m)\in [0,1]$ is a smoothing piece-wise function with a tapering parameter δ_m fully described in Appendix B of Abbott et al. (2021c) . The secondary BH mass is modeled by a power law with an index β in the domain m ∈ [m_low, m₁].

3.2.3. Summary

With the above assumptions, the cosmological and astrophysical hyperparameters to be studied are

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{\lambda }}}_{{\rm{c}}} & = & \{{H}_{0},{{\rm{\Omega }}}_{{\rm{m}},0}\}\\ {{\boldsymbol{\lambda }}}_{{\rm{z}}} & = & \{\gamma ,k,{z}_{{\rm{p}}}\}\\ {{\boldsymbol{\lambda }}}_{{\rm{m}}} & = & \{\alpha ,\beta ,{\delta }_{m},{m}_{\mathrm{low}},{m}_{\mathrm{high}},{\mu }_{{\rm{g}}},{\sigma }_{{\rm{g}}},{\lambda }_{{\rm{g}}}\}.\end{array}\end{eqnarray} \tag{ 26 }$$

The fiducial values and the prior ranges chosen for this work are reported in Table 1. For the cosmological parameters, we will consider the value of the matter density to be fixed to its fiducial value.

Table 1. Summary of the Population Hyperparameters and Priors Adopted

Parameter	Description	Fiducial Value	Prior
	Cosmology (flat ΛCDM)
H0	Hubble constant [km s−1 Mpc−1]	70.0	${ \mathcal U }(10.0,200.0)$
Ωm,0	Matter energy density	0.25	Fixed
	Rate evolution (Madau-like)
γ	Slope at z < zp	2.7	${ \mathcal U }(0.0,12.0)$
κ	Slope at z > zp	3	${ \mathcal U }(0.0,6.0)$
zp	Peak redshift	2	${ \mathcal U }(0.0,4.0)$
	Mass distribution (PowerLaw+Peak)
α	(Primary) slope of the power law	3.4	${ \mathcal U }(1.5,12.0)$
β	(Secondary) slope of the power law	1.1	${ \mathcal U }(-4.0,12.0)$
δm	(Primary) smoothing parameter [M⊙]	4.8	${ \mathcal U }(0.01,10.0)$
mlow	Lower value [M⊙]	5.1	${ \mathcal U }(2.0,50.0)$
mhigh	Upper value [M⊙]	87.0	${ \mathcal U }(50.0,200.0)$
μg	(Primary): mean of the Gaussian component [M⊙]	34.0	${ \mathcal U }(2.0,50.0)$
σg	(Primary): standard deviation of the Gaussian component [M⊙]	3.6	${ \mathcal U }(0.4,10.0)$
λg	(Primary): fraction of the Gaussian component	0.039	${ \mathcal U }(0.01,0.99)$

Note. ${ \mathcal U }(\cdot )$ denotes a uniform distribution.

Download table as:  ASCIITypeset image

For the simulation of GW detections, we use the simulation pipeline GWFAST ¹⁴ (Iacovelli et al. 2022a, 2022b).

We assume quasi-circular non-precessing BBH systems. Their waveform is characterized by the detector frame parameters

$$\begin{eqnarray}&&{{\boldsymbol{\theta }}}^{\det }=\left\{{{ \mathcal M }}_{c},\eta ,{\chi }_{1,z},{\chi }_{2,z},{d}_{L},\theta ,\phi ,\iota ,\psi ,{t}_{c},{{\rm{\Phi }}}_{c}\right\},\end{eqnarray} \tag{ 27 }$$

where ${{ \mathcal M }}_{c}$ is the detector frame chirp mass, η is the symmetric mass ratio, χ_1/2,z are the adimensional spin parameters along the direction of the orbital angular momentum, d_L is the luminosity distance, θ = π/2 − decl. and ϕ = R.A. are the sky position angles, ι refers to the inclination angle of the binary's orbital angular momentum with respect to the line of sight, ψ is the polarization angle, t_c is the coalescence time, and Φ_c is the phase at coalescence.

First of all, we generate a population of GW events following the prescriptions given in Section 3.2, so that each source is characterized by a set of source-frame parameters θ . For the parameters in Equation (27) that are not explicitly mentioned in Section 3.2, we assume the following distributions: the sky position angles are uniform on the octant covered by the parent sample, the inclination angles have a uniform distribution in $\cos \iota$ in the range [0, π], the polarization angle and coalescence phase have a uniform distribution in [0, π] and [0, 2π], respectively, and the time of coalescence, given as Greenwich mean sidereal time and expressed in units of fraction of a day, is uniform in the interval [0, 1]. For the spin parameters instead, we adopt a flat distribution between [−1, 1] for the components aligned with the orbital angular momentum. The source-frame parameters are then converted to a detector frame using the fiducial cosmology (see Table 1).

For each source, we simulate GW emission and compute the network matched-filtered signal-to-noise ratio (S/N) using the IMRPhenomHM (London et al. 2018; Kalaghatgi et al. 2020) waveform approximant, which includes the contribution of subdominant modes to the signal, that are of great relevance to describe in particular the merger phase of BBH systems. We consider two network configurations. The first, denoted as O4-like, is composed of a network of the two LIGO interferometers at Hanford, Washington, and Livingston, Louisiana, USA, the Virgo interferometer in Cascina, Italy, and the KAGRA interferometer in Japan. We assume sensitivity curves representative of the O4 run of the LVK Collaboration, with public sensitivity curves from Abbott et al. (2016). For the second configuration, denoted as O5-like the network includes the two LIGO, Virgo, and KAGRA instruments, as well as a LIGO detector located in India. For O4, we use sensitivity curves representative of the O5 run of the LVK Collaboration ¹⁵ (Abbott et al. 2016). We assume a 100% duty cycle in all cases.

Then, we select samples of GW events as follows:

• 1.  
  O4-like: 100 events with a network S/N > 12;
• 2.  
  O5-like: 100 events with a network S/N > 25.

These cuts are designed to yield the 100 best events for each configuration over approximately 1 yr of observation. We determined these numbers simulating with GWFAST a 1 yr observing run for each of the two scenarios, with a population of BBHs with an overall merger rate calibrated on the latest LVK constraints (see, e.g., Iacovelli et al. 2022a). We note that the fact that the simulation is performed on one octant of the sky is irrelevant in our case, as we do not constrain direction-dependent hyperparameters.

We approximate the GW posterior probability with the Fisher information matrix (FIM) approximation, where the GW likelihood $p({{\boldsymbol{d}}}_{i}^{\mathrm{GW}}\,| \,{{\boldsymbol{\theta }}}_{i}^{\det })$ is assumed to be given by a multivariate Gaussian distribution. This approximation is valid for the high S/N events that we are considering in this work (see Iacovelli et al. 2022a for details). We compute the FIM with GWFAST. For each event, we then draw 5000 samples from a multivariate Gaussian distribution with covariance equal to the inverse of the FIM, further imposing priors that are chosen to be flat in detector frame masses with the condition ${m}_{2}^{\det }\lt {m}_{1}^{\det }$, $\propto {d}_{L}^{2}$ for the luminosity distance, and equal to the distributions used to sample the events for the remaining parameters. We then use the samples of ${m}_{1}^{\det },{m}_{2}^{\det },{d}_{L},$ R.A., decl. to compute the GW likelihood as described in Section 2.4.

The main properties of the galaxy and GW catalogs are summarized in Figure 3. The top panels show the GW skymaps of the events detected in the two configurations overlaid to the galaxy sky distribution. The central panels present the scatter plots of the sky localization area versus the error on d_L (with a color scale giving information on the S/N). Finally, the bottom panels show the redshift and mass distribution of the detected GW events, as well as the distribution of the number of galaxies found within their localization volumes.

Zoom In Zoom Out Reset image size

The O4-like events are at redshift z ≲ 0.9 with sky areas between a few and a few hundred square degrees. This typically results in more than ∼500 potential host galaxies for 90% of the events (Figure 3, bottom-right panel). A particularly lucky event is present with a small localization area (∼2 deg²) and high S/N (∼55, see the mid-left panel of Figure 3). While this event represents an outlier of the distribution that can be ascribed to a statistical fluctuation, it still represents a possibility with this configuration. For the O5-like events, the redshift distribution remains limited to z ≲ 1 as a consequence of the higher S/N cut, while the larger and more sensitive detector network substantially improves the localization capabilities. This typically results in more than ∼50 potential host galaxies for 90% of the events (Figure 3, bottom-right panel). In this configuration, there is a significant tail of events with just a few tens of galaxies, corresponding to sky localization regions of less than 1 deg²  with an S/N that can exceed 100. Overall, while the number of galaxies in the localization volume depends on the assumption on the galaxy catalog employed, it is interesting to observe the 10× reduction in the number of potential host galaxies between the O4 and O5 networks. This improvement, following from the much smaller localization volumes, is a key factor in obtaining improved dark siren constraints as discussed in Section 4.

Ultimately, for the computation of the selection bias, we generate injection sets for both the O4-like and the O5-like scenarios with GWFAST, adopting the same selection cuts as for the GW catalogs. The injections cover the same sky area as the catalogs and a distance range that arrives up to the detector horizon for the given S/N cuts. The injection set is made of N_inj = 2 × 10⁷ and 4 × 10⁷ events, resulting in about 1.5 × 10⁶ and 1 × 10⁶ detected events in the O4- and O5-like scenario, respectively. These are then used to estimate the selection bias as in Equation (21).

We start by discussing the results for the best-case scenario, consisting of a complete galaxy catalog with spectroscopic redshift measurements. The marginal posterior distributions for the selection of parameters {H₀, α, μ_g, σ_g, γ} are shown in Figure 6.

Zoom In Zoom Out Reset image size

We find that in 1 yr of observations in the O4-like configuration, the LVK data are able to constrain H₀ with 7% uncertainty (at the 1σ level) from BBH in a combined cosmological and astrophysical population inference. This is a remarkable improvement with respect to the current constraints, as the analysis of the 42 BBHs at S/N > 12 observed so far yields a 46% measurement of H₀ (Mastrogiovanni et al. 2023). To arbitrate the Hubble tension, percent-level measurements are required. If we assume that the uncertainty on H₀ scales as $1/\sqrt{N}$, in the O4 configuration it is not possible to reach 1% in the planned schedule of about 2 yr.

In contrast, with the O5-like configuration it is possible to reach 1% uncertainty in about 1 yr of operation. We stress that this is a best-case scenario, relying on having a complete catalog of potential hosts. In general, the actual completeness can vary among different galaxy surveys and galaxy types. For example, with Euclid it would change between the photometric or the spectroscopic survey mode and the north or south direction, requiring a more detailed assessment in a future study. Moreover, even if we marginalize the astrophysical parameters, the results still rely on an assumption of the functional form of the BBH mass distribution. This dependence should be investigated in more detail in the future. On the other hand, this analysis is based only on the BBH population, and further improvements can be obtained including NSBH and BNS events and their potential EM counterparts.

We now move to the population hyperparameters. Currently, in population studies the cosmology is typically fixed (e.g., Abbott et al. 2023c). Here we study how well the fiducial models are recovered when taking into account potential correlations with cosmological hyperparameters. In Figure 7, we show the reconstructed primary mass distribution. For both the O4 and O5 scenarios, the Gaussian peak at 34 M_⊙ is clearly visible and its mean value μ_g is recovered with a precision of 7% and 3%, respectively. The second best-constrained mass parameter is the slope α of the primary BBH mass distribution that is recovered with fractional uncertainties of 15% and 13%, respectively. Overall, for the mass function parameters these are small improvements with respect to the full cosmological and astrophysical inference with the current GWTC-3 catalog, which yields ${\sigma }_{{\mu }_{{\rm{g}}}}/{\mu }_{{\rm{g}}}\sim 13 \%$ and σ_α /α ∼ 11% (see Mastrogiovanni et al. 2023). The same is true for the population analysis at fixed cosmology, which gives ${\sigma }_{{\mu }_{{\rm{g}}}}/{\mu }_{{\rm{g}}}\sim 9 \%$ and σ_α /α ∼ 11% (Abbott et al. 2023c). Similarly, the constraints on the rate parameter γ remain essentially unchanged with respect to the current knowledge, even when considering the transition from O4 to O5. This can be explained by the higher S/N threshold adopted in our O5 catalog, resulting in the same number of GW events that map a redshift range comparable to that of O4.

Zoom In Zoom Out Reset image size

In conclusion, we recall that our results are based on the full astrophysical and cosmological analysis of the best 100 GW events detectable in 1 yr for each configuration. Population studies typically benefit from the inclusion of all confident events detected and are carried out by fixing the cosmological parameters (e.g., Abbott et al. 2023c). In this sense, our population constraints should not be taken as representative of the overall performance of O4 and O5.

While obtaining a complete spectroscopic catalog poses challenges and awaits future facilities, ongoing surveys such as DES and Euclid are already building extensive photometric galaxy catalogs. The full z_phot and z_spec configurations for both the O4- and O5-like configurations are compared in Figure 4.

With photometric redshifts, the constraint on H₀ for O4 is notably less accurate, with a measurement uncertainty that is three times greater (${\sigma }_{{H}_{0}}/{H}_{0}\approx 18 \%$) compared to the spectroscopic approach. In the case of O5, this factor increases to 9 (${\sigma }_{{H}_{0}}/{H}_{0}\approx 9 \%$). Interestingly, this shows that from 100 O5 events at S/N > 25 it is not possible to achieve percent-level precision on H₀ using a photometric catalog, even under the assumption of completeness. Of course, one may consider lowering the S/N threshold to include more events; however, such events would also have much larger localization volumes, and thus a large number of potential hosts. This would likely limit their additional constraining power. Some information can be still retrieved for the mass distribution, whose reconstruction benefits from a larger sample; however, in the next section we will show that the constraints obtained in the absence of a galaxy catalog with our sample are only marginally worse than what can be obtained from a sample with a much lower S/N cut. This seems to suggest that it would be difficult to reach percent-level accuracy.

Another interesting result concerns the comparison of the full (z_phot) O5-like and full (z_spec) O4-like configurations (see Figure 5). We find that considering a galaxy catalog with spectroscopic redshift uncertainties in the O4-like scenario, we are able to achieve a more precise constraint on H₀ compared to having a larger LVK detector network at O5 design sensitivity with a photometric galaxy catalog. This occurs despite the factor 10 improvement in the GW localization volume (see Figure 3). Overall, these findings underline the importance of mapping the GW localization volume—at least for well-localized events—with dedicated spectroscopic surveys.

In the absence of the galaxy catalog information (spectral sirens case), the cosmological constraints are determined by the capability of the source-frame mass function to break the mass-redshift degeneracy.

In our spectral siren analysis, the Hubble parameter is recovered at ${\sigma }_{{H}_{0}}/{H}_{0}\approx 43 \%$ in O4-like and ${\sigma }_{{H}_{0}}/{H}_{0}\approx 32 \%$ in O5-like configuration. Our results are in good agreement with those of Mancarella et al. (2022) and Leyde et al. (2022). In particular, Leyde et al. (2022) find 38% (24%) uncertainty on H₀ using S/N > 12 spectral siren events for O4 (O5) obtaining a total of 87 (423) events. When comparing with our results for O5, we must consider that we applied a higher S/N > 25 threshold, resulting in a factor of 4 difference in the number of detected events. It is interesting to note that this large difference results only in a marginal improvement in the constraining power, suggesting that the constraints are mainly driven by the events with higher S/N.

In general, we conclude that the analysis of 100 well-localized GW events with a complete galaxy catalog with spectroscopic measurements would provide better constraints with respect to a pure spectral siren analysis under the assumed mass functions. This result also holds true in comparison to the 5 yr results of O5 spectral sirens by Farr et al. (2019), providing a relative error on H(z) of ≈3% at a pivot redshift of z ≃ 0.8. We note however that the authors assumed a much more optimistic event rate and a mass function with sharper edges, in line with the knowledge of the BBH population at the time.

In the absence of a galaxy catalog, there is a strong correlation between H₀ and μ_g, which drives the constraint on H₀. In this section, we comment on this correlation and compare the three scenarios considered in the paper. Figure 8 shows the constraints in the H₀−μ_g plane in the O5-like configuration.

Zoom In Zoom Out Reset image size

A relatively strong anticorrelation is present between μ_g and H₀ in both the spectral and full (z_phot) analyses. In terms of the Spearman's rank correlation coefficient computed from the MCMC chains, these correspond to ρ ∼ −0.8 and ρ ∼ −0.6, respectively.

This is a consequence of the fact that, in the absence of a precise redshift measurement, a higher value of H₀ would move the events with a given luminosity distance at higher z, requiring smaller values of the source-frame masses to reproduce the observed data (see Equation (4)). This leads to a smaller value of μ_g. This anticorrelation is observed also in the latest GWTC-3 analyses (Abbott et al. 2023b). Having a complete spectroscopic galaxy catalog allows instead to break the H₀−μ_g degeneracy, constraining H₀ with higher precision. This is visible in Figure 8 where no significant correlation is present in the full (z_spec) case. Interestingly, in the case where only photometric redshifts are available, this result shows that the information coming from the mass distribution still plays a role in the measurement of H₀ and has therefore to be accurately modeled. In particular, Figure 8 shows that a wrong assumption on the value of the mass scale would lead to a bias on H₀ even in the presence of a complete photometric galaxy survey.

The core modules of CHIMERA are Likelihood.py and Bias.py. In particular, the file Likelihood.py contains classes to compute the product of the integrals in Equation (17), while Bias.py allows for the computation of the selection bias term ξ( λ ). The computation of the likelihood uses functions to analyze the GW (from GW.py) and EM (from EM.py) information. These modules contain specific methods to perform all the preliminary computations that do not change during the likelihood evaluation (a more extended discussion can be found in the next paragraph). Data are loaded with specific classes present in DataGW.py (e.g., DataGWMockand DataLVK) and DataEM.py(e.g., MockGalaxiesMICEv2and GLADEPlus), respectively. The algorithm to compute the Likelihood proceeds as follows:

• 1.  
  First, the class MockLike stores all the population models (model_*), GW data (data_GW*), galaxy data (data_GAL*), pixelization parameters (nside_list, npix_event, sky_conf), and integration parameters (z_int_H0_prior, z_int_sigma, z_int_res).
• 2.  
  The GW class is initialized and the pixelization and redshift grids are precomputed. To optimize the computation in the case of large galaxy catalog analysis, the code not only restricts the sky localization, but also the redshift integration grid. The first task is done starting from posterior distributions in R.A. and decl. The user can choose the approximate number of pixels desired for each event (npix_event) to be found within a confidence level ellipse (sky_conf), given a list of possible pixelizations (nside_list). CHIMERA optimizes the pixelization of each event (in particular, the nside parameter of HealPix) to obtain the number of pixels closest to npix_event. The second task must take into account that by varying the cosmology, also ${{ \mathcal K }}_{\mathrm{gw}}$ varies, so to avoid biases, the integration grid must encompass all the redshift ranges explored during the inference. This is obtained starting from the GW posteriors on d_L and defining the range inside z_int_sigma standard deviations at a resolution of z_int_res spanning all the range ¹⁷ of H₀ explored.
• 3.  
  The EM class is initialized, and the quantity ${p}_{\mathrm{cat}}(z,\hat{{\rm{\Omega }}}\,| \,{{\boldsymbol{\lambda }}}_{{\rm{c}}})$ (Equations (11) and 12)) is precomputed pixel by pixel on the redshift grids. At this point, based on the included catalog (e.g., GLADE+, MICEv2), it is possible to activate catalog-related tasks (e.g., luminosity cut, or associate user-defined redshift uncertainties). Similarly, if defined, the completeness is computed, and the pixelized completeness function P_compl(z) (see the discussion in the above paragraph) is stored in the class to be accessible during the inference.

The pixelization approach and the precomputation of p_gal are essential for next-gen galaxy surveys. The first allows us to reduce the dimension, combining the probability of 10–10³ galaxies in one single pixel ensuring granularity in the sky analysis. The second allows us to avoid the computation of p_gal at each step of the MCMC inference. While this last approximation is not a problem for H₀ inference since both the numerator and denominator of Equation (12) have a ${H}_{0}^{3}$ dependence, the impact of a possible bias on Ω_m should be carefully assessed when future analyses using events at higher z will be carried out. The algorithm to compute the selection Bias proceeds as follows:

• 1.  
  First of all, the class Bias stores all the population models (model_*), the directory of the GW injection catalog data (file_inj), and the S/N threshold to be applied. Optional arguments also include the catalog interpolant function. If not given, the bias is evaluated for a uniform in comoving volume galaxy distribution (i.e., ∝dV_C /dz).
• 2.  
  The injection catalog is loaded by applying the chosen S/N cut. It is important to ensure that this cut is equivalent to the one adopted when creating the catalog of GW events to analyze.

Thus, following Equation (1), the full likelihood is conveniently derived by first calculating the log-likelihood for all the events and then subtracting the log-bias term, which is computed only once and multiplied by the number of events. The computation of the likelihood and selection bias is performed by calling the .compute() methods or, in logarithmic form, the .compute_ln() methods. In the latter case, the models should also be given in logarithmic form.

Currently, CHIMERA includes the following population models:

• 1.  
  Mass distributions (mass.py): logpdf_TPL (Truncated Power Law), logpdf_BPL (Broken Power Law), logpdf_PLP (Power Law + Peak), logpdf_PL2P (Power Law + 2 Peaks).
• 2.  
  Rate evolutions (rate.py): logphi_PL (Power Law), logphi_MD (Madau-like).
• 3.  
  Spin distributions (spin.py): logpdf_G (Gaussian), logpdf_U (Uniform).
• 4.  
  Cosmological models (cosmo.py): fLCDM, fLCDM modified gravity

All the above functions accept parameters in the form of dictionaries, e.g., lambda_cosmo={''H0'':70.0, ''Om0'':0.3}.